Environmental drivers of spring primary production in Hudson Bay

Over the past decades, the Arctic Ocean has undergone a significant decline in the previously dominant, thick multiyear ice (MYI) cover, leading to predictions of an ice-free (with sea-ice area < 1 million km²) Arctic summer before 2050 (Notz et al., 2020). This loss in the sea-ice cover has decreased habitat availability for Arctic top predators while increasing light availability for primary producers in the ice bottom and water column. Observations of high relative contributions of bottom-ice and subice algal communities to total annual production in the central Arctic Ocean (e.g., Gosselin et al., 1997; Boetius et al., 2013; Fernández-Méndez et al., 2014, 2015; Leu et al., 2015), followed by large under-ice phytoplankton blooms in spring (e.g., Mundy et al., 2014; Assmy et al., 2017; Oziel et al., 2019), and an increasing occurrence of secondary fall blooms (Ardyna et al., 2014) outline a productive polar ecosystem. However, these features are often localized and show high interannual variability, making predictions of the future timing and magnitude of primary production (PP) and its impacts on higher trophic levels in the changing Arctic Ocean on a pan-Arctic scale difficult. This difficulty also highlights the need to assess current seasonal PP patterns in multiple regions of the Arctic and sub-Arctic seas.

Hudson Bay, the world’s largest inland sea, at the southern margin of the Canadian Arctic, has so far received little attention during the spring peak of microalgal growth, even though it holds 10% of the seasonal ice cover found in the Arctic Ocean and provides a habitat for large populations of migratory birds and marine mammals (Ferguson et al., 2010). Furthermore, the summer ice extent has declined at a rate of –10.4 ± 3.2% per decade between 1968 and 2009, increasing the open water season by 3.1 weeks (Tivy et al., 2011; Hochheim and Barber, 2014). The current concentration and thickness distribution of the sea-ice cover, which is generally present from December to July, is controlled mainly by air temperature and wind forcing (Gagnon and Gough, 2005; Hochheim et al., 2011; Andrews et al., 2018; Kirillov et al., 2020). Especially strong northwesterly winds regularly open up a polynya in the northwest (NW polynya), which enhances ice formation in winter (Bruneau et al., n.d.; Saucier et al., 2004; Landy et al., 2017), but also makes northwestern Hudson Bay the first area to become ice-free during spring (Andrews et al., 2018). Due to the dominant northwesterly wind direction, sea ice is generally advected eastward, causing the ice cover along the east coast to grow dynamically to a thickness greater than 2 m (Prinsenberg, 1986; Landy et al., 2017). The last remaining sea ice is typically found offshore of the Hudson Bay lowlands in the southern part of the Bay (Landy et al., 2017; Kirillov et al., 2020).

Hudson Bay is also expected to undergo rapid changes in the influx of freshwater with predicted increases in precipitation and freshwater discharge from the surrounding watershed in response to the projected warming climate (Clair et al., 1998; Brown, 2010; Stadnyk et al., 2019), with major implications for PP (Hopwood et al., 2020). The inland sea already receives a river discharge of 630–870 km³ yr^–1, which corresponds to 12% of the total pan-Arctic runoff (Saucier et al., 2004; Déry et al., 2011; St-Laurent et al., 2011). This large runoff forms a strong pycnocline dividing the warmer and fresher surface layer from the underlying colder and saltier water (Prinsenberg, 1986).

In summary, these environmental conditions (ice dynamics, water circulation and stratification, and freshwater fluxes) together cause large spatial variations in phytoplankton production and biomass in Hudson Bay. Coastal areas as well as the entrance into Hudson Bay from Foxe Basin, hereafter called the Narrows, are 2–3 times more productive than central Hudson Bay in late summer to early fall (Bursa, 1961; Anderson and Roff, 1980; Harvey et al., 1997; Ferland et al., 2011; Heikkilä et al., 2014). Due to the inaccessibility of central Hudson Bay during spring, there are no previous observations of ice algal and phytoplankton production during the spring bloom. Based on historical postbloom measurements, annual production of Hudson Bay has been estimated to range from 21.5 to 39 g C m^–2 (Roff and Legendre, 1986; Jones and Anderson, 1994; Ferland et al., 2011; Bélanger et al., 2013). However, these studies likely underestimate total production, as recent observations of satellite-derived surface chlorophyll a (Chl a) concentration by Barbedo et al. (2020) found the highest phytoplankton biomass to occur in the NW polynya during the spring season. Additionally, ice algal blooms with intermediate biomass of <40 mg Chl a m^–2 (Gosselin et al., 1986; Michel et al., 1993) to very high biomass of up to 170 mg Chl a m^–2 (Welch et al., 1991) have been observed under the stable landfast ice at several locations around the Bay between March and May. Furthermore, after the ice algal bloom has sloughed off the ice, under-ice blooms have been observed below the landfast ice in June with Chl a of 2.5 mg m^–3 (Legendre et al., 1981; Michel et al., 1993). Hence, by not accounting for the spring bloom, previous estimates likely significantly underestimate the total annual production of Hudson Bay, and further analysis of this period is required.

This study addresses this shortcoming by providing the first measurements of ice algal and pelagic PP in Hudson Bay during the late spring season when incoming solar irradiance is at its seasonal maximum. We used a combination of physical and biogeochemical parameters collected during the Hudson Bay System Study (BaySys) scientific cruise onboard the Canadian Coast Guard icebreaker CCGS Amundsen in June 2018. Our objectives were to (1) characterize the environmental parameters driving spring PP, (2) investigate the community structure and photoacclimation of microalgae in the different habitats, and (3) estimate the onset and magnitude of PP by the different algal communities in relation to the melting sea-ice cover in central Hudson Bay. Results for both sea-ice and pelagic PP are presented, which demonstrate a high spatiotemporal variability of microalgal biomass and production in late spring.

The present study was conducted from June 2 to July 1, 2018, as part of the BaySys project, which aimed to understand the relative contributions of river regulation, to generate hydroelectric power, and climate change to freshwater–marine coupling in Hudson Bay. At the time of the cruise, the seasonal ice cover was still in place and unregulated river discharge was near its seasonal maximum. Water samples were collected in three regions: (1) the Narrows near the confluence of Foxe Basin, Hudson Bay, and Hudson Strait; (2) western Hudson Bay including the NW polynya; and (3) ice-covered central Hudson Bay (Figure 1). Ice cores from mobile sea ice were collected from drifting ice floes in the Narrows and central Hudson Bay. Sampling in eastern Hudson Bay was not possible due to heavy ice conditions. Additional long-term Chl a fluorescence data were recorded at approximately 32 m (2016–2017) and 28 m (2017–2018) by a mooring (AN01), which was deployed in southwestern Hudson Bay (59° 58.156’ N, 91° 57.144’ W) in September 2016, redeployed in September 2017, and recovered during our June 2018 cruise. At sampling stations, water depths ranged from 104 to 321 m in the Narrows and 31 to 185 m in Hudson Bay. At open water stations, sampling was comprised of vertical profiles of physical and biological variables including light measurements, while ice stations further included remotely piloted airborne system (RPAS) surveys of the ice floe surface, sampling of ice cores, and melt pond water for physical and biological variables as well as optical measurements above and beneath the sampled ice floe.

Sampled mobile ice floes were selected based on both spatial representation and suitability for bringing the Amundsen alongside to safely disembark the research team by an ice cage or their suitability to land a helicopter. On the ice, three to four areas of different ice-surface types were first identified for optical measurements to ensure an undisturbed snow or ice surface. The radiometer setup consisted of a surface reference for measuring incident downwelling planar irradiance, E_d(0, λ), and an under-ice arm equipped with a similar radiometer to measure transmitted downwelling planar irradiance, E_d(z, λ). Both hyperspectral radiometers (RAMSES-ACC, TriOS GmbH, Rastede, Germany) were equipped with internal pressure and tilt sensors and measured irradiance spectra in the wavelength range of 320 to 950 nm at a resolution of 3.3 nm (190 channels). Surface albedo and light transmittance were determined for the different surface types of snow, melt ponds, and white ice (i.e., snow-free ice with a white surface scattering layer). Transmitted irradiance was recorded via a hyperspectral radiometer that was attached to a custom-built double-hinged aluminum pole (hereafter L-arm; Ehn et al., 2008; Matthes et al., 2019). The L-arm was deployed through a 10-inch auger hole and positioned the radiometer directly beneath the ice bottom 1.5 m south of the hole. Snow and/or shaved ice were placed back into the hole to minimize the influence of elevated light levels on under-ice measurements.

For additional optical data processing, the fractional area for each surface type was estimated from RPAS surveys using a DJI Phantom 4 (DJI, Shenzhen, China), equipped with a 12 MP optical camera. Details on postprocessing mosaic image generation and surface type classification can be found in Harasyn et al. (2020). The surveys covered an area of 0.12 km² producing classified images of melt pond, snow-covered, white ice, and open ocean classes with 2.5-cm pixel resolution. Fractional area of melt ponds is expressed as a value of melt pond area over the total sea-ice area (sea ice plus melt pond area).

Sea-ice sampling was performed with a 9-cm core barrel (Mark II, Kovacs Enterprises, Roseburg, OR) on snow-covered and white ice areas. No ice cores were collected in melt ponds. Sea-ice thickness and freeboard were measured with an ice thickness gauge at each drilled hole through the ice floe. At each sampling location, two ice cores were extracted for vertical salinity and temperature profiles following Eicken et al. (2014). For biological sampling, the 5-cm bottom section of three ice cores were collected, pooled immediately in a dark isothermal container, and melted in 0.2-µm filtered seawater (FSW) at a ratio of 3:1 (three parts FSW and one part ice core volume) in the dark over 24 h to reduce osmotic stress (Campbell et al., 2019). Two additional independent 5-cm bottom sections were collected for bulk ice nutrient analysis. One 5-cm section was melted slowly in the dark without dilution for the analysis of silicic acid (Si(OH)₄) concentration. The other bottom section was melted rapidly in a sterile bag, which was submersed in 40 °C water to determine the concentration of nitrate (NO₃), nitrite (NO₂), and phosphate (PO₄). If melt ponds were present, water for the analysis of the same biological parameters was collected with a submersible pump (Cyclone®). Nutrient concentration of pond water was not determined.

Additionally, weekly Canadian Ice Service (CIS) ice charts provided for June 2018 were used to determine total ice concentration by different stages of development in the Narrows and central Hudson Bay at the time of ice sampling. Ice charts delineate different ice regimes with polygons that present the partial concentration (in tenths) of up to three different stages of development according to the World Meteorological Organization egg code. Stages of development considered within this study are new and young ice (<30 cm), thin (30–70 cm), medium (70–120 cm), and thick (>120 cm) first-year ice (FYI). Daily fields of sea-ice concentration were retrieved from passive microwave data (10-km OSI-430 global sea-ice concentration; http://www.osi-saf.org) and used to determine the number of days between an ice concentration falling below 15% and the day that location was sampled (hereafter referred to as days of open water, DOW). Sea-ice concentration is generally underestimated from passive microwave data sets during the melt period (Kern et al., 2020), which will introduce a slight positive bias to our estimates of DOW. We further estimated the average open water period throughout Hudson Bay from 2008 to 2018 to calculate annual PP in the open water. The number of melt days, defined as time period of surface air temperatures (SAT) > 0 °C and ice concentrations > 15%, were estimated from SAT over the central Hudson Bay extracted from the European Centre for Medium-Range Weather Forecasts (ERA-Interim) reanalysis (Dee et al., 2011).

At each open water and ice station, vertical profiles of physical and biological parameters were collected with the ship’s CTD-rosette system (CTD: conductivity, temperature, and depth). Temperature, salinity, and photosynthetically active radiation (PAR, 400–700 nm) were measured with a CTD probe (SBE-911, Sea-Bird Scientific, Bellevue, WA) and a spherical (scalar) radiometer (QSP-2300, Biospherical Instruments Inc., San Diego, CA). A surface reference (QCR-2200, Biospherical Instruments Inc., San Diego, CA), measuring incoming scalar PAR, was mounted to the ship’s main mast. In situ Chl a fluorescence, measured with the fluorometer (SCF, Seapoint Sensors Inc., Exeter, NH) attached to the rosette, was calibrated against ex situ Chl a measured in discrete water samples (see below). Additional chlorophyll fluorescence data were recorded every 15 min by ECO-Triplets (Sea-Bird Scientific, Bellevue, WA), attached to the mooring AN01, which were installed at approximately 32 m (2016–2017) and 28 m (2017–2018) and were averaged over a 24-h period. This data set was utilized to investigate the timing of the spring bloom in 2017–2018. To properly assess the absolute magnitude of the seasonal increase in Chl a was not possible due to the lack of calibration of the fluorescence sensor before or after the deployment.

The mixed layer depth, Z_m, was determined by finding the depth of the maximum buoyancy frequency (Brunt–Väisälä frequency, N²) following Carvalho et al. (2017). Before the rosette deployment, the optical depths at 100% (i.e., sea surface), 30%, 15%, 5%, 1%, and 0.2% for the water sampling were determined deploying a profiling fluorometer optical system (PNF-300A, Biospherical Instruments Inc., San Diego, CA) at the bow of the ship following Ferland et al. (2011). Water samples for the determination of nutrients, algal pigments, particulate organic carbon (POC), PP, and taxonomic composition were then collected with 12-L Niskin bottles at each optical depth and station, with the exception of no PP estimates for stations 20, 23, and 28. Nutrient samples were collected every 10 m between water depths of 0 and 100 m and every 20 m below 100 m. Water samples were prefiltered with a 200-µm mesh to avoid the influence of large grazers (meso-zooplankton) and stored in the dark containers at air temperatures of 0 °C until laboratory analyses. The prefiltration step could have influenced a slight underestimation of PP due to the potential removal of diatom and Phaeocystis colonies larger than the 200-µm mesh.

Collected hyperspectral irradiance data from ice sampling were interpolated to 1-nm steps and integrated over 400–700 nm to calculate surface albedo (or reflectance, R) and transmittance (T) for PAR. R(PAR) was calculated as the average ratio of five consecutive downwelling, E_d(0⁺, PAR, μmol photons m^–2 s^–1), and upwelling, E_u(0⁺, PAR, μmol photons m^–2 s^–1), irradiance readings. T(PAR) was calculated as the ratio of E_d(z₁, PAR) and E_d(0⁺, PAR) measured simultaneously at the ice bottom and surface, respectively. Under-ice light data were previously corrected for the larger refractive index of water compared to air.

To provide more accurate estimates of PAR at the ice bottom, regional surface albedo, $R¯(PAR)$⁠, and regional transmittance, $T¯(PAR)$⁠, which considers the spatial heterogeneity of the surface, were calculated. Following Matthes et al. (2020), $R¯(PAR)$ and $T¯(PAR)$ were calculated for each ice station with known fractions of open water, A_W, snow-covered ice, A_S, white ice, A_WI, and melt pond-covered ice, A_MP, (ΣA_i = 1) as

where R and T are the measured coefficients for each surface type. For open water, R_W(PAR) was set to the value of surface reflection at 5% (Kirk, 2011).

In the water column, the depth of the euphotic zone, Z_eu, was set at 0.2% of incident surface PAR (Ferland et al., 2011). The diffuse vertical attenuation coefficient for scalar irradiance, K_d0(PAR, m^–1) in the euphotic zone was determined by the slope of the linear regression between the natural logarithm of the measured vertical scalar rosette PAR profiles and depth. For the estimation of PP vertical scalar PAR profiles, E_d0(z₂, PAR), from 1 to 100 m, were calculated by applying Beer–Lambert’s Law:

including K_d0(PAR) and the measured downwelling scalar PAR beneath the surface, E_d0(z₁, PAR). Beer–Lambert’s Law is a commonly used approximation of PAR attenuation, despite the spectral nature of the downwelling irradiance in water (Wei and Lee, 2013), that we considered valid for the purpose of this study. Due to the artificially created open water area for the rosette deployment at ice stations, under-ice vertical PAR profiles were derived as follows:

including incident downwelling planar PAR, E_d(0⁺, PAR), from the surface TriOS measurement at the ice surface, calculated regional $T¯(PAR)$⁠, K_d0(PAR) from the vertical rosette profiles, and the average cosine for downwelling irradiance, μ_d, of 0.7 to convert planar PAR into scalar PAR at the ice bottom (z₁) following Matthes et al. (2019).

Water samples for dissolved inorganic nutrients (Si(OH)₄, NO₃, NO₂, and PO₄) were collected into acid-washed 15-ml polyethylene tubes after filtration through a 25-mm Whatman GF/F filter inserted into a filter holder to remove large particles. Nutrient concentrations were measured immediately onboard with a continuous-flow AutoAnalyzer III (Bran and Luebbe GmbH, Norderstedt, Germany) using a routine colorimetric method adapted from Hansen and Koroleff (1999). Analytical detection limits were 0.05 and 0.02 μmol L^–1 for NO₃ and NO₂, respectively, and 0.05 and 0.1 μmol L^–1 for PO₄ and Si(OH)₄, respectively. Nutrient ratios were calculated for different water depths and collected ice-bottom sections at each station. The N:P and N:Si ratios are defined here as the molar ratio of NO₃ + NO₂ to PO₄ and Si(OH)₄, respectively. Contour plots of nutrient and Chl a fluorescence were drawn using the ODV 5.1.5 software (Schlitzer, 2018).

POC was analyzed from water samples filtered onto precombusted (450 °C for 5 h) 25-mm Whatman GF/F filters. Filter blanks for each sampling station were produced by filtering 500 mL of FSW through a Whatman GF/F filter. Filters were then wrapped in tinfoil and stored at –80 °C for later analysis of POC following acidification of filters to remove particulate inorganic carbon at the University of British Columbia following the protocol of Glaz et al. (2014).

Extracted Chl a was measured with a fluorometer (10AU Field Fluorometer, Turner Designs, Sunnyvale, CA) onboard while the identification and concentration of selected algal pigments were determined by reverse-phase high-performance liquid chromatography (HPLC) after the cruise. Onboard, samples were filtered onto 25-mm Whatman GF/F filters using a vacuum pump. For fluorometric analysis, filters were subsequently soaked in 10 mL of 90% acetone at 5 °C for 18–24 h to extract Chl a. Fluorescence was measured before and after acidifying the sample with 5% hydrochloric acid (HCl, 1 N; Parsons et al., 1984); Chl a was determined from these measurements using the equations of Holm-Hansen et al. (1965). For HPLC analysis, filters were stored in 2-mL cryovials, wrapped in tinfoil, and flash-frozen in liquid nitrogen. Samples were then stored at –80 °C until analysis following Kilias et al. (2013). Pigments were extracted in 1.5 mL 100% acetone at –20 °C, homogenized (Precellys, Bertin Intruments, Montigny-le-Bretonneux, France) with glass beads and centrifuged for 5 min at 12,500 rpm in a cooled centrifuge (0 °C). The supernatant was filtered through 0.2-µm polytetrafluoroethylene filters, and samples were stored in Eppendorf tubes at –80 °C prior to analysis. Subsamples of the pigment extracts were measured with reverse-phase HPLC with a VARIAN Microsorb-MV3 C8 column (4.6 mm × 100 mm), using HPLC-grade solvents (Merck), a Waters 1525 binary pump equipped with an autosampler (OPTIMAS^TM, SunChrom GmbH, Frankfurt/Main, Germany), a Waters 2996 PDA (photodiode array detector), and the EMPOWER software. Chlorophyll, derivate and carotenoid absorption peaks were detected at 440 nm, while phaeopigments were detected at 410 nm. Pigments and derivates were identified based on retention time and the spectral properties of external pigment standards. In this study, total chlorophyll a (TChl a) concentration corresponds to the sum of Chl a and chlorophyllide a. The ratios of photoprotective carotenoids (PPC; including diadinoxanthin, diatoxanthin, violaxanthin, antheraxanthin, zeaxanthin, lutein, and β,β-carotene) to photosynthetic carotenoids (PSC; including fucoxanthin, peridinin, neoxanthin, alloxanthin, 19’-butanoyl-oxy-fucoxanthin, and 19’-hexanoyl-oxy-fucoxanthin) were also calculated following the pigment clustering of Kauko et al. (2019).

The taxonomic structure of the main protist groups for all water stations, collected ice bottom, and melt pond water samples was calculated from marker pigment ratios using the CHEMTAX Software V1.95 (Mackey et al., 1996; S. Wright, 2008). Initial pigment ratios were constrained as suggested by Higgins et al. (2011) based on microscopic examination of representative samples during the cruise and with published input matrices for ice algae (Alou-Font et al., 2013) and Arctic phytoplankton (Coupel et al., 2015; Fragoso et al., 2017) applied. Following Coupel et al. (2015), phytoplankton samples were divided into high-light surface samples (0–15 m) and low-light deep samples (16–50 m) to account for variations in pigment ratios due to light acclimation of the present phytoplankton groups. Melt pond and bottom-ice algal samples were grouped together to increase the number of samples for a successful CHEMTAX run. In the used CHEMTAX version, the initial matrices were optimized by generating 60 variants of the input ratio using the random function F = 1 + S × (R – 0.5) with a scaling factor S = 0.7 and R as a random number between 0 and 1 generated using the RAND function in Microsoft Excel as described in S. W. Wright et al. (2009). The best 10% of output matrices (n = 6) were averaged and used as new input matrix for a successive run of 60 variants of the new input matrix with S = 0.4 to reduce the standard deviation of results as recommended by Latasa (2007). The results of these six best output matrices were used to calculate the averages of the relative abundance estimates of the main protist groups. The final ratio matrices for bottom-ice algae and melt ponds (Table S2) and for phytoplankton (Tables S3 and S4) are displayed in the supplemental material.

Additionally, identification and enumeration of ice-bottom communities and phytoplankton at the subsurface chlorophyll maximum (SCM) was performed on 250-mL subsamples from melted bottom-ice scrapes and water samples. For the analysis of ice-bottom communities, the bottommost 1 cm of three ice cores was scraped off with a pocketknife into a container with FSW. Subsamples were preserved in acidic Lugol’s solution (Parsons et al., 1984) and stored in the dark at 4 °C until analysis. Cells were identified with a light microscope (Zeiss Axiovert 10 and Leica DMIL LED) following the inverted microscope method (Lund et al., 1958). Cell identification was performed to the lowest rank possible (groups, genus, or species; >2 µm) and referring primarily to Poulin and Cardinal (1982a, 1982b), Medlin and Priddle (1990), Tomas (1997), and von Quillfeldt (2001). Cell abundance was corrected for FSW dilution of ice-bottom samples.

Net primary production (NPP) of ice algal (from bottom-ice scrapes), melt pond, and phytoplankton communities were determined using the ¹⁴C assimilation method and applying photosynthesis–irradiance (P-E) relationships. Water samples in 1000-mL opaque Nalgene bottles were inoculated with initial NaH¹⁴CO₃ concentrations between 0.2 and 1.0 µCi mL^–1 depending on the strength of the Chl a fluorescence signal during the rosette cast and the length of the incubation. Out of each sampling bottle, subsamples of 50 mL were transferred to 12 clear culture flasks and one opaque flask, which were placed in a custom-made incubation chamber adapted after Babin et al. (1994). In the incubator, bottles were arranged in a row with the first bottle closest to the light source (7/9/15 W EIKO LED light bulb) and the dark bottle the furthest to provide a light gradient from 860 to 0 μmol photons m^–2 s^–1. They were incubated at –1.6 °C for 2–4 h. Three vials were also filled with 20 µL of the sample, 50 µL of ethanolamine, and 500 µL of MilliQ water to measure the initial activity and to determine the exact concentration of ¹⁴C in the samples. At the end of the incubation, samples were filtered onto 0.2-µm Millipore filters, and the filters were transferred into 20-mL scintillation vials to be spiked with 300 µL of 3.16% HCl. Vials were placed open on an orbital shaker for 2 h to evaporate the remaining inorganic ¹⁴C on the filter under a fume hood. Afterward, vials were filled with 10 mL of EcoLume Scintillation Cocktail (MP Biomedicals, Santa Ana, CA). The particulate radioactive carbon uptake was counted after the cruise at Université Laval using a Tri-Carb 2910 TR scintillation counter (PerkinElmer, Waltham, MA). The carbon uptake values in the opaque flask were subtracted from the corresponding clear flask’s carbon uptake values.

Samples for dissolved inorganic carbon (DIC), which was needed in the calculation of the amount of labeled carbon incorporation into the cell, were taken directly from the Niskin bottles and melt ponds into 250-mL or 500-mL borosilicate glass bottles with ground-glass stoppers and secured with electrical tape. All DIC samples were poisoned with 100 µL of a saturated HgCl₂ solution to halt biological activity and were stored in the dark at room temperature until being processed ashore. DIC was measured with a Single-Operator Multiparameter Metabolic Analyzer. The DIC concentration in the collected 5-cm ice-bottom core sections was not measured. Instead, DIC was calculated using the measured salinity of the core section and the equation presented in Parsons et al. (1984). Calculated carbon fixation rates (P^B, mg C mg^–1 Chl a h^–1) were normalized to measured Chl a and photosynthesis–irradiance relationships (P-E curves) were fitted by minimizing the sum of differences between the measured carbon uptake and the model proposed by Platt et al. (1980).

where $PsB$ (mg C mg^–1 Chl a h^–1) is the maximum carbon fixation rate if there is no photoinhibition, β^B, mg C mg^–1 Chl a h^–1 (μmol photons m^–2 s^–1)^–1; α^B, mg C mg^–1 Chl a h^–1 (μmol photons m^–2 s^–1)^–1, is the photosynthetic efficiency, defined as the initial slope of the P-E curve; and E (μmol photons m^–2 s^–1) is the irradiance measured in the incubation chamber. Only P-E curves with R² = .9 were included in the further analysis. Maximum carbon fixation rate $PmaxB$ was calculated as

The photoacclimation parameter, E_k, was calculated as $PmaxB$ divided by α^B. Production rates (mg C m^–3 h^–1) for each station were calculated by multiplying the P-E parameters of six optical depths with the vertical profiles of E_d0(z₂, PAR) for each hour of the day (24 h), which were generated from the performed light measurements and the change of the sun’s position over the day. Hourly production rates were then integrated over Z_eu and over the day using trapezoidal integration to calculate daily production rates (mg C m^–2 d^–1). Although short incubation times of 2–4 h were used to measure production, which does not account for respiration during nighttime and recycling of ¹⁴C fixed by photosynthesis, we consider our results to be only slightly different from NPP due to the integration of production over the euphotic zone. Prior studies have shown that the ¹⁴C-method with short incubation times provides good estimates of NPP at low growth rates, which was likely the case in the light-limited lower euphotic zone (Pei and Laws, 2013, and citations therein). Furthermore, short-term incubations minimize the potential for algae to acclimate to the constant light conditions in the incubator (Lewis and Smith, 1983). TChl a and nutrients were also integrated over Z_eu using trapezoidal integration. Mean integrated nutrient concentrations in the euphotic layer were obtained by dividing depth-integrated values by the integration depth.

Total annual PP of microalgal communities was estimated from historic field measurements and results from this study. Due to lack of direct PP measurements in early spring, Chl a of bottom-ice algae and under-ice phytoplankton were extracted from the literature, and production was calculated as net accumulation over the sampling period with a POC: Chl a ratio of 54 (Irwin, 1990). The incubation times of direct PP measurements differ between this study and historic estimations (24-h on-deck incubations in Ferland et al., 2011; 10-h on-deck incubations in Lapoussière et al., 2013). In the next calculation step of annual PP, seasonal production for early spring was calculated by multiplying the daily average of total production by ice algae and phytoplankton by 92 days. Late spring production during sea-ice melt was calculated by multiplying the daily average of total production by phytoplankton, bottom-ice algae, and Melosira arctica, measured in this study, with 34 melt days (i.e., where SATs were above 0 °C and ice concentration was >15%). Seasonal production in the ice-free water in summer and fall was calculated by multiplying the daily average of phytoplankton production with the average of 146 open water days between 2008 and 2018. PP during winter (December to February) is assumed to be negligible and was not included in the annual estimate.

A principal component analysis (PCA) was carried out on collected physical data to identify clusters of regions within Hudson Bay with similar physical environmental parameters. Included parameters were K_d0(PAR), Z_eu, Z_m, mean temperature, T_m, and salinity, S_m, in the mixed layer, integrated nutrient concentrations (Si(OH)₄, NO₃ + NO₂, PO₄) in the euphotic zone, ice concentration gained from the CIS ice charts for June 2018 and the DOW prior to sampling. The PCA was performed with the stats package in the R 5.5.1 software. Significant differences between the P-E curve parameters of phytoplankton communities in the different environments (open water and under-ice) and at different depths (surface: 0–15 m, deep: 16–50 m) were investigated using a two-way analysis of variance (ANOVA) in the R software packages “car” and “dplyr”. A log transformation was performed to achieve normal distribution of the data set and the homogeneity of variances for the use of parametric tests. Differences in P-E parameters of ice-associated communities were not statistically tested due to the low number of measured P-E curves. A two-way ANOVA was also used to investigate differences between the ratios of POC:Chl a and PPC:PSC of phytoplankton in the different environments and depths. Tukey's range test was performed to investigate the interaction between the groups further if significant results were identified during the ANOVA. Differences between the ratios of POC:Chl a and PPC:PSC of ice-associated communities (bottom-ice algae and melt pond) were tested with a student’s t test.

Based on the geographical location of sampling stations and the presence/absence of an ice cover, all sampling stations were grouped into three regions with distinct environmental conditions: (1) the partially ice-covered Narrows; (2) open water in western Hudson Bay, including the NW polynya; and (3) ice-covered central Hudson Bay (Figure 1). Open water stations close to the coast and in the NW polynya in western Hudson Bay were characterized by a depleted nutrient concentration in the upper euphotic zone and a warmer and deeper surface mixed layer (Table S1). Inshore stations also had the lowest surface salinities due to their proximity to river estuaries with the PCA analysis outlier, station 46, being located in front of the large and turbid estuaries of the Nelson and Hayes rivers (Figure 1). Stations in central Hudson Bay showed a higher nutrient concentration in the euphotic zone as well as a colder and shallower mixed layer compared to western Hudson Bay. An increased light attenuation, K_d0(PAR), was also observed at stations 32, 34, and 40 that were located east of the estuary.

To highlight the varying ice conditions in Hudson Bay, ice stations were separated into three subregions: (1) the Narrows, stations sampled in early June; (2) north-central Hudson Bay, stations sampled in mid-June; and (3) south-central Hudson Bay, stations sampled in late June (red polygons, Figure 2). The Narrows had an ice concentration of 67%, mainly composed of thick FYI that had no visible signs of surface melt (Figure 2b). Ice in the Narrows had a mean thickness of 114 ± 29 cm (hereinafter mean ± standard error), a freeboard of 9 ± 1 cm, and a snow depth of 13 ± 6 cm, along with the coldest (–1.7 ± 0.1 °C) and saltiest (5.8 ± 0.2) ice observed. In comparison, sea-ice concentrations were higher in north- and south-central Hudson Bay with more medium and thinner FYI being present (Figure 2c and d). Mean ice thickness was lowest in north-central Hudson Bay (75 ± 7 cm and freeboard of 5 ± 1 cm) where negative freeboard was observed at a few floes. Sea-ice concentration and thickness (128 ± 17 cm and a freeboard of 16 ± 2 cm) was highest in south-central Hudson Bay where the ice was much more deformed, and several ice floes were thicker than 2 m. Additionally, the ice in central Hudson Bay was in an advanced melt stage with high melt pond coverage, and the ice itself was warmer and less salty with mean ice temperatures of –0.9 ± 0.1 °C and –0.8 ± 0.1 °C as well as mean bulk salinities of 3.6 ± 0.2 and 1.9 ± 0.2 in north- and south-central Hudson Bay, respectively.

The differences in ice thickness and state of decay directly impacted the optical parameters of the ice cover in the Narrows and central Hudson Bay. The observed decrease in $R¯(PAR)$ and increase in $T¯(PAR)$ throughout the sampling period matched the observed ice-surface melt progression (Figure 2). Although the areal fraction of more transparent melt ponds increased, $T¯(PAR)$ remained in the same range in south-central (0.01–0.40) compared to north-central (0.07–0.27) Hudson Bay due to the thicker ice cover.

Differences in the water column structure between the regions are presented as potential temperature-salinity diagrams of the vertical CTD profiles (Figure 3) and along transects in the three regions (Figure 4). The Narrows were characterized by a surface water layer with S_m and T_m of 31.9 ± 0.4 °C and –1.5 ± 0.1 °C, respectively (Figure 3a). The deep water (>100 m) was saltier, but in the same temperature range, resulting in a weakly stratified water column with a Z_m at 22 ± 5 m. In western Hudson Bay, the mixed layer shoaled from 20 ± 4 m, measured in the center of the NW polynya, to 10 ± 2 m, measured inshore. This surface mixed layer was characterized by S_m and T_m of 31.9 ± 0.3 °C and 0.4 ± 0.3 °C, respectively, and was ice-free (<15% ice concentration) for an average of 25 days prior to sampling. The deep water in the center of the NW polynya was the coldest and saltiest observed in the entire Hudson Bay with S > 33.1 and T < −1.7 °C below 100 m (Figure 3b). In central Hudson Bay, the observed vertical salinity gradient followed the seawater freezing point (Figure 3c). Similar S_m and T_m of 31.2 ± 0.2 °C and –1.4 ± 0.1 °C, respectively, were measured throughout central Hudson Bay, while the mixed layer of 23 ± 4 m was deeper in the north compared to the mixed layer of 13 ± 2 m in the south.

Nutrient concentrations in the euphotic zone and in the deep waters differed among the three regions and are shown along a transect in each region (Figures 1b and 4). Transect 1 in the Narrows extends across the mouth of Foxe Basin and the strait between Southampton Island and Coats Island. Transect 2 in western Hudson Bay extends from the western shore of Hudson Bay across the area of open water (NW polynya) and into the western edge of the ice pack. Transect 3 in central Hudson Bay extends from the outer Nelson River estuary into the thicker ice pack of central Hudson Bay.

Along Transect 1, the NO₃ + NO₂ concentration ranged from 3.34 to 9.09 μmol L^–1, with the highest concentrations in the bottom waters of the Narrows (measured 10 m above the sea floor; Figure 4e and Table 1). Concentrations of PO₄ and Si(OH)₄ ranged from 0.85 to 1.07 μmol L^–1 and 9.03 to 16.1 μmol L^–1, respectively (Figure 4f and g). Overall, nutrient concentrations in the euphotic zone were higher in the Narrows with mean surface N:P and N:Si molar ratios of 5.05 and 0.44, respectively, compared to the euphotic zone across Hudson Bay with mean surface N:P molar ratios between 0.16 and 1.73 and N:Si molar ratios between 0.03 and 0.22 (Table 1).

In the euphotic zone of the NW polynya along Transect 2 in western Hudson Bay, concentrations of NO₃ + NO₂, PO₄, and Si(OH)₄ ranged from 0.01 to 3.22, 0.41 to 0.90, and 0.01 to 9.01 μmol L^–1, respectively (Figure 4l–n; Table 1), with inshore Si(OH)₄ concentrations near the detection limit. The nitracline depth (NO₃ + NO₂ < 1 μmol L^–1) largely tracked the depth of the mixed layer (Figure 4l), extending to 30 m at station 28 in the polynya and shoaling toward the ice edge. Nutrient concentrations in the deep waters below 100 m remained high along the transect with NO₃ + NO₂, PO₄, and Si(OH)₄ concentrations ranging from 4.47 to 12.9, 0.75 to 1.62, and 9.21 to 35.7 μmol L^–1, respectively. Concentrations of NO₃ + NO₂ in the deep water at stations 20, 21, and 24 in western and central Hudson Bay were even higher than the observed nitrogen inventory of the deep waters in the Narrows (Figure 4e and l).

Along Transect 3 in central Hudson Bay, the nitracline depth exhibited a similar pattern to Transect 2, being deepest in the open water and shoaled toward the ice edge (Figure 4s). In the euphotic zone, integrated NO₃ + NO₂ concentration increased from 0.52 μmol L^–1 at station 46 to 2.05 μmol L^–1 at station 38. Concentrations of NO₃ + NO₂, PO₄, and Si(OH)₄ in the ice-covered euphotic zone ranged from 0.44 to 5.73, 0.64 to 1.01, and 3.16 to 13.0 μmol L^–1, respectively (Figure 4s–u; Table 1). Concentrations of NO₃ + NO₂ and Si(OH)₄ in the bottom water were difficult to compare across Transect 3 as water depth varied greatly from 31 m at station 32 to 178 m at station 38. In general, concentrations of NO₃ + NO₂ (Si(OH)₄) of 13.1 (38.5) μmol L^–1 in the deepest waters were comparable to those observed in western Hudson Bay (Figure 4l, n, s, and u).

Concentrations of NO₃ + NO₂, PO₄, and Si(OH)₄ in the bottom ice of the mobile ice cover were higher in the Narrows compared to central Hudson Bay (Table 1). Ice stations within central Hudson Bay further showed a spatial gradient of 2-fold higher bottom-ice nutrient concentrations in the north compared to the south. Overall, bottom-ice nutrient concentrations were low with smaller N:P and N:Si molar ratios of 1.04 and 0.50, respectively, than N:P and N:Si molar ratios in the underlying surface water (Table 1).

The spatial distribution of TChl a and PP in the water column largely reflected the vertical gradients of nutrient concentration in the different regions with high production estimates being associated with low nutrient concentrations. Within the Narrows, TChl a was low in the euphotic zone with values < 1 mg m^–3, although Z_eu reached a depth of 41 ± 7 m (Figure 4d and Table 2). Integrated daily NPP values for phytoplankton in this region were the lowest observed during the study with a mean value of 98.4 ± 18.2 mg C m^–2 d^–1.

In western Hudson Bay, a strong SCM was observed between 9 and 50 m, usually between Z_m and Z_eu (Figure 4k; Table 2). The strongest SCM with TChl a between 2.6 and 4.7 mg m^–3 was observed in the center of the NW polynya, resulting in a higher K_d0(PAR) and a slightly shallower Z_eu of 38 ± 4 m compared to that of the Narrows (labeled “Integration Depth” in Table 2). However, inshore stations 19 and 22 were characterized by a low K_d0(PAR) of 0.12 m^–1, a deep Z_eu of 49 m, and low TChl a of <1 mg m^–3. Inshore station 46, which was located near the Nelson River estuary, differed from these characteristics with a K_d0(PAR) of 0.23 m^–1, a Z_eu of 24 m, and TChl a of >1 mg m^–3 attributed to a phytoplankton bloom in the Nelson River estuary (Jacquemot et al., 2021). NPP varied largely in the open water between 170 mg C m^–2 d^–1 at station 22 and 803 mg C m^–2 d^–1 at station 17 with a mean NPP of 460 ± 70 mg C m^–2 d^–1 (Table 2).

Phytoplankton TChl a measured beneath the ice cover in central Hudson Bay exceeded 0.5 mg m^–3 throughout the euphotic zone with highest concentrations between 2.7 and 4.0 mg m^–3 at station 25. These high TChl a values in the north-central Hudson Bay resulted in high K_d0(PAR) similar to those in the NW polynya (Table S1) and the highest estimated NPP of 1,400 mg C m^–2 d^–1. Mean phytoplankton NPP in central Hudson Bay was calculated at 414 ± 146 mg C m^–2 d^–1 with the lowest NPP in the south due to TChl a that only exceeded 1 mg m^–3 in the shallow surface mixed layer (Figure 4r). However, NPP increased further into the ice-covered area from 128 mg C m^–2 d^–1 at station 40 to 391 mg C m^–2 d^–1 at station 36.

NPP measured in the ice bottom in the Narrows and central Hudson Bay and in the evolving melt ponds on the ice surface within Hudson Bay were minimal due to low TChl a (Table 2). Mean NPP was highest in the ice bottom in the Narrows. Overall, the combined contribution of ice algal and melt pond communities to late spring PP in Hudson Bay accounted for less than 1% during this study. In contrast, the observed subice diatom M. arctica contributed 30% to late spring production (Table 2). M. arctica was observed visibly as long strands attached to the ice bottom at stations in the Narrows and north-central Hudson Bay, and microscopically as small chains in the ice-bottom samples from south-central Hudson Bay, except at station 38 where long strands were observed visibly. Samples of M. arctica for biological analysis were collected at stations 25 and 38. The measured TChl a of 13.7 ± 0.8 mg m^–2 and an assumed M. arctica mat thickness of 5 cm, estimated from videos taken of the bottom ice through an ice hole, resulted in a NPP of 378 ± 119 mg C m^–2 d^–1. This estimate does not account for the observed patchiness of M. arctica aggregates due to limited sampling. However, much of the M. arctica biomass sloughed from the ice bottom upon extraction of ice cores (Figure S1), making our NPP estimate conservative. Because of these contrasting issues, we have kept the NPP value as the best current estimate but strongly suggest that future work seek to refine the contribution of M. arctica to spring production in Hudson Bay.

CHEMTAX results, which calculate the relative contribution of each algal group to Chl a, and results from the inverted light microscopy suggest a flagellate-dominated phytoplankton community in the Narrows, with a particularly high relative abundance of unclassified flagellates (including prymnesiophytes, raphidophytes, and choanoflagellates) in the deeper water layers of the euphotic zone between 16 and 50 m (Figure 5). Diatoms made up less than 33% of the relative contribution to the main protist groups in the Narrows. Within Hudson Bay, the open and ice-covered water column was dominated by diatoms with a relative contribution of more than 61% on the surface and 64% in the deeper layers of the euphotic zone in the CHEMTAX analysis. The relative contribution of unclassified flagellates decreased to less than 17% in surface water between 0 and 15 m and 13% between 16 and 50 m. Cryptophytes, chrysophytes, and prasinophytes were present in relative contributions below 14%, 7%, and 19%, respectively, at all stations, while chlorophytes were only sparsely detected in the calculated pigment ratios. The microscopic analysis showed that centric and pennate diatoms were similarly abundant in the SCM with 28% and 26%, respectively, in western Hudson Bay and 19% and 24%, respectively, in central Hudson Bay. In western Hudson Bay, the most abundant centric diatoms were Chaetoceros gelidus (3.9 × 10⁶ cells L^–1) and Thalassiosira nordenskioeldii (2.3 × 10⁶ cells L^–1), while Fragilariopsis cylindrus (5.0 × 10⁶ cells L^–1) and Fragilariopsis oceanica (6.5 × 10⁶ cells L^–1) were the most abundant pennate diatoms. In central Hudson Bay, the most abundant centric diatoms were Chaetoceros spp. (0.8 × 10⁶ cells L^–1) and T. nordenskioeldii (0.7 × 10⁶ cells L^–1). The most abundant pennate diatoms were F. cylindrus (1.1 × 10⁶ cells L^–1) and Nitzschia frigida (1.2 × 10⁶ cells L^–1).

Bottom-ice algal communities within Hudson Bay were dominated by diatoms with a mean relative contribution of 92% (Figure 5) to Chl a of the major algal groups in the CHEMTAX analysis. The microscopic analysis of the bottom-ice community revealed a similar high mean relative abundance of diatoms at 82% of all cells enumerated. Pennate diatoms were especially abundant with a mean relative abundance of 66%. The most abundant pennate diatom was N. frigida (261.8 × 10⁶ cells L^–1), while Chaetoceros spp. (75.3 × 10⁶ cells L^–1), Thalassiosira spp. (29.0 × 10⁶ cells L^–1), and M. arctica (31.2 × 10⁶ cells L^–1) were abundant centric diatoms. Melt pond communities were also dominated by diatoms with a relative contribution of 53% in the CHEMTAX analysis but were overall more diverse with a larger relative contribution of cryptophytes (18%), unclassified flagellates (15%), and prasinophytes (14%) compared to their relative contributions to Chl a in the bottom-ice algal communities. Sampled melt ponds were not connected to the underlying water column, and salinities were between 0.2 and 4.1.

The P-E parameters varied between the microalgal communities in the different habitats (Table 3). $PmaxB$ of phytoplankton in the open water was significantly higher (F_1,75 = 4.53, p < .05) than that beneath the ice cover. E_k and α^B were not significantly different, and photoinhibition was only observed in a few under-ice surface samples. Depth influenced all three P-E parameters with significantly higher $PmaxB$ (F_1,75 = 5.55, p < .05), significantly lower α^B (F_1,75 = 5.29, p < .05), and significantly higher E_k (F_1,75 = 36.49, p < .001) in the surface water. E_k was significantly higher in the open (p < .05) and in the ice-covered surface water (p < .001) compared to those of the deeper water layers in the respective environments.

P-E parameters of M. arctica were in the same range of phytoplankton in the ice-covered surface water. However, $PmaxB$ and α^B of bottom-ice algae and melt pond communities were 3–10 times lower compared to under-ice phytoplankton. Only E_k of the bottom-ice algae was similar to that of the under-ice communities. Melt pond communities at the ice surface showed the highest and lowest E_k and α^B, respectively, and high β^B, which was not measured in other ice-associated communities.

Water column POC:Chl a ratios (wt:wt) were significantly different between different environments (F_1,94 = 5.27, p < .05) and depths (F_1,94 = 8.15, p < .01). Mean open water ratios were 175 ± 38 (median = 134; Figure 6) in the surface and 134 ± 39 (67) in the deeper layer. Under-ice POC:Chl a ratios were significantly lower at 87 ± 9.0 (72) in the surface and 67.4 ± 11.3 (52.2) in the deeper layer. The regression of POC versus Chl a (data not presented) showed no statistically significant relationships between POC and Chl a in western Hudson Bay. The y-intercept of POC versus Chl a relationships of stations in south-central HB, although significantly different from zero, were low or even negative. Mean POC:Chl a ratios of 341 ± 123 (153) for ice bottom and 401 ± 105 (312) for melt ponds were not significantly different.

The ratio of photoprotective to photosynthetic carotenoids (PPC:PSC, wt:wt) was not significantly different between phytoplankton in the open water and ice-covered water column (Figure 6). However, ratios decreased significantly with depth (F_1,97 = 43.5, p < .001) with measured ratios of 0.28 ± 0.04 (median = 0.27; Figure 6) in the open surface water and 0.17 ± 0.01 (0.17) in the deeper water. Mean ratios of under-ice phytoplankton were 0.24 ± 0.10 (0.23) in the surface and 0.16 ± 0.01 (0.16) in the deeper water. PPC:PSC ratios of ice-associated communities were significantly higher (t_1,11 = –7.14, p < .001) in melt ponds with a mean ratio of 1.63 ± 0.30 (2.11) compared to the ice bottom with a mean ratio of 0.27 ± 0.03 (0.27). Furthermore, bottom-ice algal communities had a higher mean PPC:PSC ratio than under-ice phytoplankton communities.

Time series of Chl a fluorescence at the lower SCM depth (28–32 m) was recorded by mooring AN01 (Figure 1) in southwestern Hudson Bay to gain more information about the timing of PP in the water column (Figure 7). In 2017, the ice cover (>8/10 concentration, CIS ice charts) was present until early July. Chl a fluorescence had already begun to increase in the fully ice-covered surface water layer in the beginning of June. During the following open water season of the same year, Chl a fluorescence decreased, which could have been related to the formation of a deeper SCM observed at 37 m in late June 2018. In 2018, Chl a fluorescence also increased while the ice cover was still present. However, maximum Chl a fluorescence was measured in the open water column following an earlier ice breakup in early June at the mooring location.

The observed large differences in spring PP, biomass (TChl a), and phytoplankton community composition between the Narrows and western and central Hudson Bay are in line with previous observations during summer and fall. The main factors influencing these various regions are differences in freshwater input, nutrient concentrations, light conditions, and distance from shore (Bursa, 1961; Anderson and Roff, 1980; Harvey et al., 1997; Ferland et al., 2011; Heikkilä et al., 2014).

In western Hudson Bay, surface phytoplankton communities benefited from a continuously open latent-heat polynya in early May, that thereby increased underwater light availability and promoted strong surface stratification through solar heating, as well as contributions from ice melt. Relatively high surface Chl a (>1.2 mg m^–3) was observed by satellite in late May 2018 within the first 3 weeks after the ice breakup (Barbedo et al., 2020). At the time of sampling in mid-June 2018, the region had been ice-free for 25 days, providing more than enough time for a surface bloom to nearly deplete NO₃ + NO₂ and Si(OH)₄ in the surface mixed layer and form a strong SCM (Figure 4). PO₄ was still available throughout the euphotic zone in the entire Hudson Bay following the Redfield ratio of 16 N:1P (Redfield, 1963) and thus was not limiting algal growth anywhere.

In early spring, PP in the NW polynya benefits from replenished surface nutrient concentrations brought up by vertical mixing during the winter months (Tremblay et al., 2019). The enhanced ice formation and brine production in the NW polynya (Bruneau et al., n.d.; Landy et al., 2017; Kirillov et al., 2020) can overcome stratification and deeply mix the water column to depths of 100 m by the end of winter in the region (Prinsenberg, 1986; Saucier et al., 2004). Indeed, deep water in the center of the NW polynya was the coldest and saltiest observed during our study. These waters were further characterized by high concentrations of inorganic nutrients (Table 1), which likely accumulate in the deep interior of the Bay due to the small water exchange with the adjacent marine water bodies and the long residence time of deep waters between 4 and 14 years within Hudson Bay (Pett and Roff, 1982; Tremblay et al., 2019). With the deep winter mixing potential, this pool of nutrients can likely help to increase surface production within the NW polynya.

Rivers draining into western Hudson Bay, with the largest contributors being Chesterfield Inlet in the northwest and Churchill, Hayes, and Nelson rivers in the southwest, have not been shown to supply substantial additional inorganic nutrients during late spring to summer (Déry et al., 2011; Tremblay et al., 2019). During our study, several coastal stations (17, 18, 19, 22, 46; Figure 1b) lay within 30 to 75 km from shore and were influenced by the large cyclonic coastal buoyancy current that carries freshwater along the coast (Prinsenberg, 1983; Granskog et al., 2007; Déry et al., 2011; St-Laurent et al., 2011). Salinities decreased in the surface mixed layer from 32.3 in the north (station 18) to 31.6 and 29.8 at southern stations (22 and 46, respectively). This boundary current reaches up to 100 km offshore and creates a fresh, thick (5–25 m) summer mixed layer overlaying a colder subsurface layer formed during winter mixing (Granskog et al., 2009). The investigated coastal stations during this study were characterized by a shallow and fresher mixed layer of 12.0 ± 2.9 m thickness and very low nutrient concentrations. In the center of the NW polynya, the mixed layer was 22.5 ± 5.3 m. Thus, riverine input likely decreases the potential for coastal PP in this region by adding a buoyant, nutrient-depleted surface layer, particularly after phytoplankton deplete surface nutrients originally replenished via winter mixing processes.

However, several studies reported an inshore-offshore gradient of higher biomass found inshore with values between 0.2 and 1.0 mg Chl a m^–3 versus lower biomass found offshore with values between 0.1 and 0.5 mg Chl a m^–3 in summer (Anderson and Roff, 1980; Roff and Legendre, 1986; Harvey et al., 1997; Granskog et al., 2007; Ferland et al., 2011). During summer and fall, strong tidal and wind-driven mixing can weaken surface stratification and, in combination with the entrainment of deeper salt water and accompanying nutrients into the freshwater plume via estuarine circulation, lead to increased production inshore (Kuzyk et al., 2009; Ferland et al., 2011). In late spring, the inshore-offshore gradient was reversed with a lower TChl a between 0.3 and 1.4 mg m^–3 in the euphotic zone inshore compared to higher TChl a between 1.4 and 4.9 mg m^–3 in the euphotic zone of the NW polynya.

Although TChl a in the SCM was high in the center of the NW polynya, late spring PP was driven by phytoplankton in the surface layer. Production in the SCM, which generally occurred below the mixed layer depth near the nitracline and was associated with the 0.2%–1% optical depth at 40 ± 4 m, only contributed 1%–9% to total production, assuming an SCM thickness of 5 m. A well-developed SCM, often found at similar optical depths and between 20 and 60 m, is characteristic of central Hudson Bay in the summer and fall (Roff and Legendre, 1986; Harvey et al., 1997; Granskog et al., 2007; Ferland et al., 2011; Lapoussière et al., 2013). However, estimated late spring PP of 460 mg C m^–2 d^–1 of the diatom-dominated phytoplankton community in western Hudson Bay was higher than the estimated production of 322 mg C m^–2 d^–1 in summer (Ferland et al., 2011) and of 100 mg C m^–2 d^–1 in fall (Lapoussière et al., 2013), which was dominated by smaller cells (0.7–5 μm). We conclude that the bloom was likely past its peak, although integrated phytoplankton biomass and PP in the NW polynya were still greater than those in the Narrows and central Hudson Bay (Table 2).

Phytoplankton production in the Narrows and central Hudson Bay was driven by the formation of open water through ice export in the Narrows and by the sea ice melt and increasing melt pond formation at the ice surface in central Hudson Bay, which contributed largely to the increase in under-ice light levels, a deepening of the euphotic zone and surface stratification. Phytoplankton communities within the Narrows appeared to be in a prebloom to early bloom stage with observed low biomass and NPP, which were likely the result of density instabilities in surface waters due to freezing air temperatures (Oziel et al., 2019) and stronger tidal mixing at the southern end of Foxe Basin (Drinkwater and Jones, 1987). This early stage had little impact on surface nutrient concentrations in the Narrows, which remained relatively high throughout the water column. Later in the season, after increasing air temperatures and sea ice melt produce a more stabilized surface mixed layer (Drinkwater and Jones, 1987), these relatively nutrient-replete waters create favorable conditions for a phytoplankton bloom (Ferland et al., 2011). Previously observed late summer NPP of 371 mg C m^–2 d^–1 (Ferland et al., 2011) in the Narrows was four times higher than our measured early June NPP of 98.4 mg C m^–2 d^–1. Furthermore, the late summer production presented by Ferland et al. (2011) was also driven by a diatom-dominated community, while the spring phytoplankton community observed in our study contained a large fraction of flagellates, particularly in the water column below 15 m, which is more typical of a prebloom stage (Norrbin et al., 2009).

In central Hudson Bay, a diatom-dominated under-ice phytoplankton bloom was observed. TChl a was high throughout the euphotic zone with no distinct SCM as nutrients were still available in the surface layer with NO₃ + NO₂ concentrations just below 2 μmol L^–1. Under-ice NPP in the euphotic zone was highly variable with a greater NPP of 612 mg C m^–2 d^–1 in north-central Hudson Bay compared to 215 mg C m^–2 d^–1 in south-central Hudson Bay. Several environmental conditions may have caused these regional differences. North-central ice stations (16, 21, 24, 25, and 29) were in proximity of the incoming polar surface and Atlantic water through Foxe Basin and Hudson Strait, respectively, which represents an external nutrient source for an ice-edge/under-ice bloom. Satellite observations suggested moderate surface Chl a (0.2–0.5 mg m^–3) immediately after the breakup followed by a decreasing trend as the season progressed (Barbedo et al., 2020). This pattern is consistent with our in situ observations.

The south-central ice stations 32 and 34 were only 44–65 km away from shore and were characterized by a shallow mixed layer with a low surface salinity (Table S1) indicating the influence of the previously mentioned coastal buoyancy current. This current also carries an elevated CDOM concentration, particularly in the south (Granskog et al., 2009), which could explain the observed high PAR attenuation of 0.19 and 0.27 m^–1 at stations 32 and 42, respectively. The high K_d0(PAR) in combination with measured low $T¯(PAR)$ of 0.01 through ice floes thicker than 2 m resulted in a shallow Z_eu and, subsequently, low NPP. Barber et al. (2021) further described a vast area of thick (>10 m), heavily deformed sediment-laden sea ice in this region, which, with its thickness, prolongs ice melt until August and could limit light availability and ultimately PP during spring and summer in this area. Additionally, the calculated low values for integrated surface nutrient concentrations over the euphotic zone in south-central Hudson Bay (Figure 4) indicate an overall lower potential for under-ice production compared to the ice-covered northern region. This phenomenon of low surface nutrient concentrations at the beginning of the sea ice melt could be a function of the localized cyclonic circulation of water with lower nutrient concentrations in this region (Ridenour et al., 2019).

Previous studies on landfast ice in southeastern Hudson Bay reported the formation of under-ice blooms after freshwater from snow and ice melt stabilized the water column in late May (Legendre et al., 1981; Runge et al., 1991; Michel et al., 1993). The observed blooms reached maximum Chl a between 1.5 and 2.7 mg m^–3 in the surface water, which is similar to the TChl a of 1.8 mg m^–3 that we observed in the ice-covered surface mixed layer in June 2018. The Chl a fluorescence sensor attached to mooring AN01 at 30 m detected an increase in Chl a fluorescence at the beginning of June 2017, highlighting an early onset of under-ice PP (Figure 7). A similar trend of an under-ice Chl a accumulation was observed at the mooring site in 2018. However, the ice broke up a month earlier in early June 2018, which fueled a phytoplankton bloom in the open water at the ice edge (Barbedo et al., 2020). Under-ice blooms occur in Hudson Bay as evidenced from our study and the historical record of blooms beneath landfast ice. However, considering the calculated mean integrated TChl a of 35.10 mg m^–2 over the ice-covered euphotic zone in central Hudson Bay (Table 2), phytoplankton biomass was comparable to the central Arctic Ocean, but much lower than under-ice blooms in the Arctic shelf regions (Ardyna et al., 2020).

The investigation of the state of photoacclimation of the phytoplankton communities showed that communities displayed greater light (shade) acclimation near the surface (deeper waters). In the open water, surface communities synthesized more PPC, displayed in the significantly greater PPC:PSC ratio, that dissipate excess light energy via nonphotochemical quenching (Hill et al., 2005; Alou-Font et al., 2016; Joy-Warren et al., 2019; Kauko et al., 2019) compared to communities in the deeper layer of the euphotic zone. However, the significantly greater POC:Chl a in the open water surface layer cannot necessarily be attributed to a lower amount of light-absorbing pigments due to the potential for the increased contribution of detritus to POC during late bloom stages. Nevertheless, these acclimation mechanisms help to explain the greater $PmaxB$ of surface communities than that in the SCM. Our observations are also consistent with Huot et al. (2013) who found decreasing $PmaxB$ with depth in the Beaufort Sea and in the Canadian Archipelago.

In the ice-covered surface layer, the phytoplankton community was acclimated to the reduced light; however, with increasing light levels this surface community displayed a greater E_k, lower α^B, and a higher PPC:PSC ratio compared to phytoplankton found in the deeper ice-covered water. $PmaxB$ as well as POC:Chl a ratios of the surface community were similar to those observed during the large under-ice phytoplankton bloom in the Chukchi Sea (Palmer et al., 2013; Arrigo et al., 2014). This observed capacity to acclimate to the variable light conditions in the different environments demonstrates considerable plasticity of the photosynthetic apparatus of phytoplankton over large spatial scales and is in line with observations of Arctic phytoplankton by Palmer et al. (2011) and Lewis et al. (2019).

In late spring, three ice-associated communities, namely, melt-pond algae, bottom-ice algae, and the subice algae with varying contributions to the late spring PP, were identified in central Hudson Bay. NPP of bottom-ice algal and melt pond communities were insignificant compared to NPP of the subice algae M. arctica, which was found in large, but patchy quantities growing attached to the bottom of ice floes mainly in the north-central region and the Narrows. This subice algal species benefits from relatively high light transmission through melting sea ice while having access to surface water nutrients through its long filaments. It has also been found to significantly increase local PP in the otherwise marginally productive central Arctic (Gutt, 1995; Gosselin et al., 1997; Fernández-Méndez et al., 2014, 2015). Subice algae could play a similar key role in carbon export in central Hudson Bay, given that our conservative estimates of NPP were on the same order of magnitude as the rates of the observed under-ice phytoplankton bloom. Filament samples showed a biomass of 13.7 ± 0.8 mg Chl a m^–2, which corresponds to the lower end of the Melosira aggregate biomass range of 14–44 mg Chl a m^–2 sampled in the central Arctic (Fernández-Méndez et al., 2014). However, due to the sporadic sampling, quantifying the biomass and production of M. arctica in central Hudson Bay was not possible. Images from the ice edge showed extensive coverage, however, highlighting the need for future investigation of the role of M. arctica during the spring bloom (Figure S1).

Bottom-ice algal communities had a much lower biomass compared to previous observations in landfast sea ice (Gosselin et al., 1986; Welch et al., 1991; Michel et al., 1993) and were likely already in a postbloom state with partial biomass loss through ice bottom melt, reflected also in the relatively high sea-ice temperatures. The measured low molar nutrient ratios of N:P and N:Si of 1.04 and 0.50, respectively, in the ice bottom as well as the high POC:Chl a ratios of ice-algal cells suggest a strong nitrogen depletion (Gosselin et al., 1990; Campbell et al., 2016; Dalman et al., 2019) and are typical of a postbloom scenario in the Canadian Arctic (Niemi and Michel, 2015). Much higher biomass was observed previously between March and May in landfast ice with Chl a between 27.0 and 170.0 mg m^–2 in northwestern Hudson Bay near Chesterfield Inlet (Bergmann et al., 1991; Welch et al., 1991) and between 23.6 and 39.7 mg m^–2 in southeastern Hudson Bay (Freeman, 1982; Gosselin et al., 1986; Michel et al., 1993), suggesting that sea ice can play an important role in the overall carbon budget of Hudson Bay.

Despite being in a nutrient-limited postbloom stage, bottom-ice algae were well acclimated to the high light levels during melt pond formation at the ice surface. Elevated concentrations of PPC with PPC:PSC ratios even higher than in under-ice phytoplankton communities in the surface water as well as a significantly higher E_k were found throughout the sampled mobile ice cover and correspond to acclimated Arctic ice algal communities during advanced melt stages (Michel et al., 1988; Mundy et al., 2011; Galindo et al., 2017). However, the observed mean PPC:PSC ratio of 0.27 during the postbloom stage was much lower than previously reported postbloom ratios in the Canadian Arctic (up to 1–3.5 in Alou-Font et al., 2013; up to 0.81 in Galindo et al., 2017) and were only found in the melt pond samples. As our bottom-ice algal communities were not photoinhibited despite the relatively high under-ice light levels, we conclude that ice algae have the opportunity to photoacclimate and reduce susceptibility to photoinhibition (Michel et al., 1988; Juhl and Krembs, 2010). Nevertheless, $PmaxB$ was much lower than previously observed in the landfast ice (Gosselin et al., 1985; Gosselin et al., 1986; Michel et al., 1988; Bergmann et al., 1991) and could be explained by an additional nitrogen limitation (Campbell et al., 2016).

Melt pond communities were subject to even higher light levels near the ice surface and, therefore, showed the highest E_k and PPC:PSC ratios of all microalgal communities in this study and high photoinhibition. $PmaxB$ and E_k were in the range of P-E parameters measured in melt pond algae in the Arctic Ocean (Lee et al., 2012; Fernández-Méndez et al., 2015). During melt pond formation, sea-ice algae can get trapped at the surface where they need to adapt rapidly to the changing conditions of high light levels, variable salinities, and potential nutrient limitation as observed elsewhere (Mundy et al., 2011; Fernández-Méndez et al., 2015; Sørensen et al., 2017). Mundy et al. (2011) further observed a high abundance of flagellates, which overlaps with the findings of this study. Overall, the contribution of melt pond communities to late spring PP in Hudson Bay was inconsequential due to low biomass and low NPP. This conclusion is in contrast to observations on MYI and FYI in the central Arctic, where measured melt pond algal biomass was up to eight times higher with daily production rates of 0.8–60 mg C m^–2 (Fernández-Méndez et al., 2015), resulting in a contribution of 1%–10% to total NPP in the central Arctic (Lee et al., 2012; Fernández-Méndez et al., 2015).

Figure 8 summarizes existing data on ice algal and phytoplankton production in the open and ice-covered water column from direct field measurements and satellite-derived Chl a (data can be found in Table S5). Total particulate annual production of microalgal communities was estimated at 72 g C m^–2 yr^–1 in Hudson Bay and represents the sum of seasonal production in early spring (March to May) and during the spring melt (June) and ice-free period (July to November; Table S6). Growth season of bottom-ice algae in the peripheral landfast sea ice starts in March, while an increase in under-ice phytoplankton Chl a was measured in May. PP during the sea-ice melt is driven by phytoplankton in the open and ice-covered water column with a significant contribution of M. arctica in central Hudson Bay. Our estimate shows that approximately 32% of annual biomass is produced during the 34-day melt period.

Seasonal production in the ice-free water represents 57% of annual production and is supported by a lengthening of the growth season to 146 open water days between 2008 and 2018 compared to an estimated growth season of 120 days in previous annual PP estimates (Ferland et al., 2011). This finding is in line with the observation of an increase in PP on a pan-Arctic scale (Arrigo and van Dijken, 2015). Satellite-derived daily production rates were not included in the estimation of seasonal production as these rates seem to largely underestimate production in the open water (Figure 8). Overall, our updated estimate of annual production is almost twice as high as annual estimates of 24–39 g C m^–2 yr^–1 based on postbloom summer and fall measurements (Roff and Legendre, 1986; Jones and Anderson, 1994; Ferland et al., 2011) and satellite-derived annual rates of approximately 20–25 g C m^–2 yr^–1 for the open water season (Bélanger et al., 2013), but in the range of modeled annual PP of 50–80 g C m^–2 yr^–1 (Sibert et al., 2011).

This study has revised the total estimated annual PP in Hudson Bay from 21.5–39 g C m^–2 yr^–1 to 72 g C m^–2 yr^–1 by including the first measurements of PP in late spring. This estimate also includes the first scientific observations of the subice diatom M. arctica in Hudson Bay. Removing the contribution of M arctica due to uncertainties in its estimate still results in an annual PP estimate of 59 g C m^–2 yr^–1 for Hudson Bay. The diatom-dominated spring bloom is driven by phytoplankton production in the surface layer beneath the melting ice cover in central Hudson Bay and in the SCM in the open water of western Hudson Bay. The measured high production rates in the ice-covered and open water thereby highlight the considerable plasticity of phytoplankton photosynthetic performance in the variable light environment of the Hudson Bay Complex. However, capturing the peak in production and biomass by the different microalgal communities is challenging due to the spatiotemporal variability in the environmental factors. In this study, we were not able to quantify the contribution of bottom-ice algae to PP in central Hudson Bay because by the time we reached the sampling area in mid-June, the ice algal community was already in a postbloom state. Instead, our observations have shown that the thin mobile ice cover in the north-central region provides a favorable habitat for M. arctica, which has the potential to significantly contribute to spring PP in Hudson Bay.

Climate-induced trends toward earlier sea-ice breakup and delayed freeze-up will likely have a negative impact on habitat availability for ice-associated communities such as M. arctica and may shift peak production earlier in Hudson Bay. An extended open water season will further increase the amount of light and heat received in the surface water in spring and will lead to changes in the timing of the phytoplankton bloom. While the spring bloom may develop earlier in the year, the longer open water season in fall combined with the projected increase in wind speeds in the Hudson Bay region (Steiner et al., 2013) could enhance mixing and result in greater access to the deep nutrient pool in the Bay. Freshwater discharge into Hudson Bay is projected to increase considerably, particularly in winter and spring, due to increased air temperature and precipitation (Stadnyk et al., 2019). This freshwater addition in winter counters the addition of brine from sea-ice formation in polynyas and leads (Eastwood et al., 2020), resulting in reduced mixing and thus a reduced replenishment of the surface nutrient inventory during winter. Ultimately, such a change could lead to a decrease in spring PP in the NW polynya, which we have shown is the largest regional contributor to annual production in Hudson Bay. This freshwater addition also highlights the possibility to use Hudson Bay as a small-scale model system for the entire Arctic Ocean to investigate the interplay of increasing freshwater buoyancy input and the increase in turbulent mixing processes caused by an intensification of storms, strong tides, and brine rejection during sea-ice formation and their impacts on future nutrient availability and PP potential in Arctic surface waters.

To gain more knowledge about the response of microalgal communities to the rapidly changing environmental conditions, the marine environment of Hudson Bay needs to be monitored more frequently with annual resolution. In the future, more autonomous observing systems such as moorings, autonomous underwater vehicles (e.g., gliders), or drifting buoy systems could be deployed in the three key regions presented here to collect year-round and multiannual data sets of biogeochemical cycles, especially in the winter–spring and summer–fall transitions when sea ice is present.

Data sets within this article have been archived in the Canadian Watershed Information Network (CanWIN): http://lwbin-datahub.ad.umanitoba.ca/organization/baysys; http://lwbin-datahub.ad.umanitoba.ca/dataset/lisa-matthes-thesis (DOI: http://dx.doi.org/10.34992/pdgt-ss90).

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

Table S1. Environmental variables of 23 stations sampled in Hudson Bay in June 2018.

Table S2. Initial and final (after CHEMTAX optimization) pigment to chlorophyll a ratios for sea-ice algae.

Table S3. Initial pigment to chlorophyll a ratio for each phytoplankton group.

Table S4. Final (after CHEMTAX optimization) pigment to chlorophyll a ratio for each phytoplankton group.

Table S5. Historical and measured seasonal primary production of microalgal communities in Hudson Bay presented in Figure 8.

Table S6. Seasonal and annual primary production in Hudson Bay.

Figure S1. Melosira arctica growing attached to the bottom of first-year sea ice in central Hudson Bay (photo credit: L. A. Dalman).

This study is a contribution to the Arctic Science Partnership (ASP). The authors thank the officers and crew of the Canadian Research Icebreaker, CCGS Amundsen for their excellence in support of this study. They further thank Rachel Hussher, Julie Major, Jonathan Gagnon, and Atreya Basu for their support during sample collection. They thank Marcel Babin for the permission to use his lab equipment for the incubation experiments. They also thank the two anonymous reviewers for their valuable comments on this article.

This study is a contribution to the Natural Sciences and Engineering Council of Canada (NSERC) Collaborative Research and Development project: BaySys (CRDPJ 470028-14). Funding for this study, including field studies, was provided by NSERC, Manitoba Hydro, ArcticNet, Ouranos, Hydro Quebec, the Canada Excellence Research Chair (CERC), and the Canada Research Chairs (CRC) programs. Further equipment funding was provided by ABOL-CFI to C. J. Mundy and J. K. Ehn. Individual support from NSERC has been provided to D. G. Barber, D. G. Babb, C. J. Mundy, and J. K. Ehn. I. Peeken was funded by the Polar Regions and Coasts in a Changing Earth System (PACES II) Program of the Helmholtz Association. Additional support was provided to L. C. Matthes and M. Harasyn by the University of Manitoba Graduate Student Fellowship and to D. G. Babb by the Canadian Meteorological and Oceanographic Society.

The authors declare that they have no conflict of interest.

Contributed to the original data acquisition: LCM, LAD, DGB, MH, SK, JL.

Contributed to the analysis and interpretation of the data: LCM, LAD, DGB, CJM, JKE, IP, JL, SK.

Drafted the article: LCM, CJM, JKE, IP.

Revised the article: All authors.

Approved the submitted version for publication: All authors.

How to cite this article: Matthes, LC, Ehn, JK, Dalman, LA, Babb, DG, Peeken, I, Harasyn, M, Kirillov, S, Lee, J, Bélanger, S, Tremblay, J-É, Barber, DG, Mundy, CJ. 2021. Environmental drivers of spring primary production in Hudson Bay. Elementa: Science of Anthropocene 9(1). DOI: https://doi.org/10.1525/elementa.2020.00160

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

Knowledge Domain: Ocean Science

Part of an Elementa Special Feature: BaySys