Seasonality of zooplankton communities in the Hudson Bay system

Arctic and sub-Arctic seas are characterised by strong seasonal cycles that have major implications for our understanding of these ecosystems. The Hudson Bay system is a large sub-Arctic region at the southern limit of the Canadian Arctic that includes the shallow Hudson Bay at its heart, the Foxe Basin in the north (also very shallow) and the Hudson Strait connecting it to the North Atlantic (the Labrador Sea) in the east.

Despite the high density of marine predators such as polar bears, ringed seals, beluga whales, bowhead whales and seabirds (e.g., Strong, 1989; Lunn et al., 1997; Dueck and Ferguson, 2009; Chambellant, 2010; Gaston et al., 2012; Castro de la Guardia et al., 2017; Smith et al., 2017), the Hudson Bay ecosystem has been known, at least until recently, to have low primary and secondary production giving the impression of a top-heavy food chain (Harvey et al., 2001; Ferland et al., 2011; Sibert et al., 2011). However, a recent study showed that primary production in the Hudson Bay system is higher than previously estimated when the spring bloom is included in the calculation (Matthes et al., 2021). Similarly, studies that estimated low zooplankton biomass in the Hudson Bay system were conducted when the bay was largely ice-free (Harvey et al., 2001; Estrada et al., 2012) and might have missed a significant proportion of the zooplankton biomass earlier in the season.

Coastal zooplankton communities in late summer and autumn are characterised by high abundances of small boreal copepod species in the east and meroplankton in the west, while low numbers of large Arctic-adapted copepods of the genus Calanus are found in Hudson Strait and in the centre of Hudson Bay (Estrada et al., 2012). In autumn along the Nunavik coast copepod assemblages are dominated by the large and lipid-rich Calanus spp. in Hudson Strait, by a mix of relatively small Arctic and boreal species in the middle section, and by euryhaline species such as Acartia longiremis and Centropages hamatus in the southern part (Harvey et al., 2001).

In this study, we characterised the mesozooplankton assemblages of the Hudson Bay system using data collected between June and October in five different years spread between 2005 and 2018. Because data on zooplankton composition and abundance in the Hudson Bay system are scarce and zooplankton have never been sampled throughout a complete seasonal cycle in any given year, we assembled existing data to form this unique data set. Although the data set remains somewhat sparse, it is the most extensive zooplankton data set from the Hudson Bay system to date, in both spatial and temporal dimensions. The seasonal succession of the environmental conditions in the system is triggered directly or indirectly by the sea ice break-up, which itself is characterised by a marked interannual variability (median ice break-up date from 1980 to 2014 had a standard deviation of 14 days with a directional trend of earlier ice break-up in recent years; Andrews et al., 2018). Previous studies have shown that there are regional variations in zooplankton assemblages within the Hudson Bay system, but these might also be influenced at least partially by the different timing of ice retreat in different regions. Here we hypothesise that the ice break-up date is the main environmental driver for the seasonal distribution of mesozooplankton assemblages in the Hudson Bay system. We performed multivariate analyses to examine how the seasonal succession of mesozooplankton communities may relate to this and other significant environmental drivers.

The ice-free season in the Hudson Bay system is becoming longer and the ice break-up and freeze-up dates are becoming more variable from year to year (Hochheim and Barber, 2010; Andrews et al., 2018). Therefore, characterising the seasonal variation in zooplankton assemblages is needed to establish a baseline against which future changes can be compared.

In the highly seasonal Hudson Bay system, environmental variables affecting the zooplankton assemblages change acutely from spring to autumn. Notwithstanding its low latitude, the Hudson Bay system goes through an annual cryogenic cycle and is completely ice-covered from late November to early July (Prinsenberg, 1986). The ice starts to melt along the west coast of the bay in the area of the Kivalliq polynya in May (Figure 1; Bruneau et al., 2021) and the remaining ice completely melts by mid-July in the south-east (Saucier et al., 2004). River discharge into the Hudson Bay system contributes, on average, the equivalent of a 0.8 m layer of freshwater annually, most of which flows into the bay in spring and summer (St-Laurent et al., 2011). Starting in spring, surface water in the Hudson Bay system, but especially in Hudson Bay and James Bay, is highly influenced by the high volume of warm and low salinity riverine input, making the bays more stratified as the summer progresses (Drinkwater and Jones, 1987; Ferland et al., 2011; Ridenour et al., 2019a). The Hudson Strait in the open-water season is more mixed than the Hudson Bay due to the stronger southeastward and cross-channel currents and intense tidal mixing (Drinkwater, 1986). Tidal mixing in the Hudson Bay is restricted to coastal regions, and vertical mixing of the surface and deep waters is rare in the inner Hudson Bay in mid-summer (Roff and Legendre, 1986). Deep water enters Hudson Bay from the Atlantic along the northern coast of Hudson Strait and from the Arctic, through Foxe Basin and via the Roes Welcome Sound (Drinkwater, 1988; Stewart and Lockhart, 2005; Defossez et al., 2010; Ridenour et al., 2021). Deep water circulation is slow and moves in a cyclonic direction (Prinsenberg, 1986). Surface circulation changes between spring and autumn. In spring, a weak double gyre system pushes surface water through the centre of the bay and some of it is then circulated back into the bay through a cyclonic gyre in the northeast. The circulation becomes simpler and largely cyclonic in autumn (Figure 1a; Ridenour et al., 2019b).

Five expeditions of the research icebreaker CCGS Amundsen were conducted in the Hudson Bay system between 2005 and 2018 (Figure 1b). The expeditions were carried out in different months (Table 1) and followed different cruise paths with some spatial overlap between years.

At 65 stations, mesozooplankton were sampled immediately after deployment of a CTD-rosette system (described below). The sampler, consisting of four nets each having a square aperture of 1 m², was lowered into the water column and then towed vertically upwards at a speed of 30 m min⁻¹ from a depth of 10 m above the bottom to the surface. The mesozooplankton samples used in this study came from one of the four nets, a 4.5 m long, square-conical net with a mesh size of 200 µm and a rigid, detachable cod end. Onboard, zooplankton samples were preserved in a borax-buffered seawater solution of 4% formaldehyde.

At each station, a CTD-rosette system, equipped with an SBE19plus^® conductivity, temperature, depth sensor (CTD), a Seapoint^® chlorophyll fluorometer, and 24 Niskin bottles, was deployed from the surface to 10 m above the seafloor. Chlorophyll a (Chl-a) concentrations in the water samples from the Niskin bottles were used to calibrate the in-situ Chl-a measured via fluorimetry (detailed methodology in Matthes et al., 2021).

The mixed layer depth (MLD) was estimated using a relative variance method. This method compares the standard deviation of salinity along the whole depth profile normalised by the maximum variation. MLD is determined as the point along the depth profile where the relative variance is at a minimum (Huang et al., 2018). The MLD estimation was done on R version 4.1.2 using the “oce” package (Kelley and Richards, 2021).

The annual date of ice break-up was determined for each station by extracting data from regional ice charts produced by the Canadian Ice Service. The first Monday of the week during which the ice concentration dropped below 50% was used as the ice break-up date. Then, for each zooplankton sampling station, the number of days between the ice break-up date and the zooplankton sampling event was calculated as the variable “days of open water” (DOW). This variable can be negative in the cases where the icebreaker managed to access areas with ice concentration above 50%.

In the laboratory, mesozooplankton were separated into a small fraction (200–1000 µm) and a large fraction (>1000 µm), resuspended in water, and each fraction was split using a Motoda-type plankton splitter until around 300 zooplankton remained in the sub-sample. Each individual in the sub-sample was then identified under a stereomicroscope to the species or to the lowest possible taxonomic level. Calanus species were differentiated based on prosome lengths specific to the developmental stage (Forest et al., 2011). The copepod developmental and nauplii stages were also determined for biomass calculation. Advanced stages of Copepoda spp. nauplii (N3–N6) were identified to level of order (Cyclopoida) or genus (Calanus, Microcalanus, Metridia, Oithona, Paraeuchaeta, Pseudocalanus). Young nauplii (N1–N3) were identified as Copepoda spp. nauplii. The taxon Copepoda spp. also includes other copepod individuals of all stages that were degraded and impossible to identify. The abundance was estimated in terms of individuals (ind) m⁻² by multiplying the number of individuals found in the subsample by the appropriate fraction that was taken from the complete sample.

The mass of each copepodite (in grams C) was estimated from its development stage using published species-specific empirical relationships between stage, prosome length and carbon content (Table S1; Mumm, 1991; Hopcroft et al., 2005; Forest et al., 2011). In 2018, two size cohorts of Pseudocalanus were observed, suggesting the presence of at least two distinct species. Therefore, the prosome lengths of Pseudocalanus spp. were measured to allow for more accurate biomass estimates. Lengths of Chaetognatha were measured in the laboratory, and their biomass was estimated from their total length and published empirical relationships between total length and carbon content (Table S2). Other macrozooplankton taxa were rare, partly because macrozooplankton species (such as Euphausiacea and large Amphipoda) are able to avoid vertically trawled nets; these were enumerated during taxonomical analyses but excluded from biomass calculations (Table S3).

The species abundances per station were converted into proportional abundances per station. The proportions were arcsine square-root transformed to reduce the impact of dominant species and give more weight to the less abundant species. This data transformation helps to better detect trends in distribution and abundance. To analyse the similarities in species composition between stations, the stations were grouped into clusters according to their species composition using a Bray-Curtis dissimilarity index and then plotted in a dendrogram. To test whether the environmental variables were significantly different between clusters, a Kruskal-Wallis rank sum test was done on each environmental variable. When a significant effect was detected (p-value < 0.05), a pairwise Wilcox test was performed to identify significantly different clusters.

The taxa most associated with each cluster were determined by calculating the indicator values of each taxon. The calculation for indicator values takes into account the relative frequency of occurrence and relative average abundance of a taxon in the cluster based on analysis proposed by Dufrene and Legendre (1997). The equation used in this study, through the R package “labdsv,” was adapted by Roberts (2019). For each cluster, the taxa with high indicator values are more frequently found and are more abundant in the stations from that cluster compared to the rest of the stations.

To analyse the relationship between taxa and environmental variables, we performed a redundancy analysis (RDA) on the taxa proportions within the whole community data set. The RDA uses multiple linear regressions between taxa and environmental variables that, as a result, rank the environmental variables in terms of importance in explaining the variation of the taxa (Kenkel, 2006). First, environmental variables were tested for co-linearity and when detected (bivariate relationship with a correlation coefficient larger than 0.75), the variable that was of lesser interest to the scope of the study was removed. Taxa that appeared in fewer than three stations were removed from the analysis to simplify the data set as our primary interest was in taxa with strong correlations to environmental variables. The taxa proportions were then Hellinger-transformed, and an RDA was performed with all the environmental variables shown in Table 2, plus the DOY of each station. As a post-hoc test to get the most parsimonious model, forward selection of each environmental variable was done using the adjusted R² as a measure of goodness of fit. The R package “vegan” (Oksanen et al., 2020) was used to perform the RDA. The version of R (R Core Team, 2021) used was 4.1.2 for all of the analysis.

The physical properties of the water column differed among years and among months (Table 2). As the shallow and deep stations were sampled in every month (although only two stations were sampled in August), the median depth sampled is similar in the different months. Temperature and salinity varied between months, respectively increasing and decreasing as the summer progressed, with June having significantly lower temperature and salinity than most other months (Kruskal-Wallis test = 34.6.5, df = 4, p-value < 0.001; Kruskal-Wallis test = 27.0, df = 4, p-value < 0.001, respectively). Again, August stations were an exception due to their high temperatures and low salinity. The MLD was shallower as the season progressed, but only significantly greater in June (Kruskal-Wallis test = 10.2, df = 4, p-value < 0.05). Chl-a varied but not significantly, and there was no clear trend from June to October. The number of DOW before sampling was significantly greater in the later months (Kruskal-Wallis test = 35.1, df = 4, p-value < 0.001), but there was a lot of variation in DOW between stations sampled within the same month. This variation was most noticeable in July, which was sampled in 2010, 2017 and 2018. Some July stations in 2010 had been ice-free for much longer than July stations in 2017 and 2018. This difference resulted in greater variability in salinity, temperature and Chl-a in July than in the other months sampled, which justifies our approach of considering carefully the DOW as a potential explanatory variable. These results show that different environmental conditions were encountered in different months, especially in June. They do not show a progression from month to month, which can be attributed to different months being sampled in different years at different locations.

A total of 92 zooplankton taxa were captured and identified in the Hudson Bay system over the five years of sampling (Table S3). The most diverse group were copepods with 36 taxa identified. Cnidaria and Amphipoda were also relatively diverse with 13 and 10 taxa, respectively.

The number of taxa per station varied between a minimum of 11 in northeast Hudson Bay, close to Hudson Strait in 2018, and a maximum of 35, along the east coast of the Hudson Bay in 2017. The most abundant taxa in terms of numbers were Oithona similis, Pseudocalanus spp., Microcalanus spp., Triconia borealis and Cirripedia. Only O. similis and Pseudocalanus spp. occurred in 100% of the stations sampled. The copepods Calanus glacialis, T. borealis, Microcalanus spp. and Metridia longa and the chaetognath Parasagitta elegans appeared in over 90% of the stations. Seventy-six of the 92 taxa identified (84%) appeared in fewer than 50% of the stations, and 54 taxa (59%) appeared in fewer than 10% of the stations, meaning that a majority of the taxa appeared in only a few stations.

The five most important taxa in terms of biomass were, in order of total biomass: Paraeuchaeta glacialis (which only appeared in 11 of the 65 stations but is a large species—females are at least twice as large as Calanus glacialis), C. hyperboreus, Metridia longa, C. glacialis and Pseudocalanus spp. Although few in number, large copepods constituted most of the mesozooplankton biomass. The overall average biomass per station was 2.9 g C m⁻², ranging from 3.5 × 10⁻³ to 109.2 g C m⁻² (or an average of 0.02 g C m⁻³, ranging from 1.2 × 10⁻⁴ to 0.5 g C m⁻³).

Cluster analysis of the 65 stations based on the proportional abundance of zooplankton species in each station revealed seven clusters at 40.2% dissimilarity. The clusters include three large clusters named after their main characteristic: deep cluster (DC), coastal-spring cluster (CSC), and coastal-autumn cluster (CAC) containing 58 of the stations; and four minor clusters (MC1, MC2, MC3 and MC4) containing seven outlier stations (Figure 2a).

The stations within the DC seemingly did not have much in common in terms of geographical location, sampling month and DOW. However, DC stations were significantly deeper than those from the other clusters, with depth being the only environmental variable that varied significantly from both the CSC and CAC clusters (median depth of 177 m; Kruskal-Wallis test = 18.7, df = 2, p-value < 0.001; Figure 2a and Table 3). The stations in the CSC and CAC were all coastal stations, but the stations in the CSC were sampled early in the season, about three weeks after ice break-up, while stations in the CAC were sampled late in the season, with a median DOW of 108 days (Kruskal-Wallis test = 19.3, df = 2, p-value < 0.001; Figure 2a and Table 3). Moreover, the CSC had the highest Chl-a concentration, which was the only variable other than depth that differed significantly between the stations in the DC and the CSC (Kruskal-Wallis test = 13.4, df = 2, p-value < 0.01). The DC had a large range of Chl-a concentrations, while the CAC had the lowest Chl-a of the three major clusters.

In summary, the stations in the DC were deep, while the stations in the CSC and CAC were shallow and coastal. The stations in the CSC varied from the rest of the stations by having higher Chl-a concentration. The stations in the CAC were sampled late in the season and had a significantly lower surface salinity and a significantly higher surface temperature than the rest of the stations (Kruskal-Wallis test = 13.1 and 15.9, df = 2 and 2, p-value < 0.01 and < 0.001, respectively).

When testing for co-linearity of the environmental variables, salinity and temperature were found to be highly co-linear (correlation coefficient: 0.81). We decided to keep salinity over temperature because salinity was of higher interest for this study of the Hudson Bay system. While surface temperature reflects the marked seasonal change in heat in the system, variations in salinity also integrate information on river discharge and sea-ice dynamics. The most parsimonious RDA model was the one including all the variables. No variables were excluded in the post-hoc analysis. The RDA model was statistically significant (Figure 2b; F = 3.9, df = 6, p-value = 0.001). Of the six RDA components in the model, RDA1 and RDA2 explained a statistically significant amount of variability in taxa between stations (RDA1: F = 9.6, df = 1, p-value = 0.001 and RDA2: F = 7.6, df = 1, p-value = 0.003). The RDA showed that the day of the year was the variable that explained the most variation in mesozooplankton composition between stations. The day of the year was highly correlated to the first RDA axis with a biplot score of 0.89. DOW was the second most important variable to explain variation between taxa, with a biplot score on the first RDA of 0.61. With a biplot score of 0.93 on the second RDA, depth was the third variable to explain a notable amount of variation. Salinity, MLD and Chl-a explained very little variation in mesozooplankton assemblages among stations. The three main dendrogram clusters were clearly separated on the RDA plot. CSC and CAC stations separated along the first RDA, while stations in the DC were grouped on the positive end of the second RDA.

The indicator taxa of the DC were Microcalanus spp. (MISP) and Triconia borealis (TRBO). On the RDA (Figure 2c) these indicator taxa, especially Microcalanus spp., were associated with deeper stations that were sampled in spring and early summer when salinity and Chl-a concentrations were high and the MLD was still deep. The indicator taxon of the CSC was Calanus spp. nauplii (CASP), which, like the stations in the CSC, was associated with the negative end of the first RDA, indicating that it was more commonly found in stations that were sampled early in the season soon after ice break-up. The indicator taxon of the CAC, Acartia sp. (ACSP), was not well correlated with the first two RDA axes. However, Acartia longiremis (ACLO), Pseudocalanus spp. (PSSP) and Bivalvia (BI) were associated with stations sampled later in the year at the positive end of the first RDA. The large copepods Metridia longa (MELO), Calanus glacialis (CAGL) and Calanus hyperboreus (CAHY) were closely associated with the second RDA, indicating that they were more common in deeper stations. Cirripedia (CI) was the only taxon that was highly associated with shallow stations sampled early in the year.

The DC contained 72 different taxa with a mean of 19 taxa per station. The most abundant taxa were Oithona similis, Pseudocalanus spp., Microcalanus spp. and Triconia borealis (Figure 3a). Among the DC, CSC and CAC, the DC had the highest average abundance and biomass at 220,695 ind m⁻² (ranging from 64,458 to 544,120 ind m⁻²) and 1.5 g C m⁻² (ranging from 0.1 to 7.0 g C m⁻²), respectively. Despite their relatively low numbers, the copepods Calanus hyperboreus, Metridia longa and Calanus glacialis contributed to most of the biomass (Figure 3b). These same species were most abundant in stations sampled in June (in 2018). In fact, of the ten stations with a biomass higher than 1.5 g C m⁻², only two were sampled in September (in 2005). Stations located in Hudson Strait had a higher biomass than stations located in Hudson Bay.

Fifty-two taxa were found in the CSC but, on average, the stations were more taxa-rich than the ones in the DC with a mean of 24 taxa per station. Pseudocalanus spp., Oithona similis and Cirripedia were the most abundant taxa, while Calanus glacialis, C. hyperboreus and Pseudocalanus spp. contributed the most to the biomass (Figure 4). Another large copepod, Calanus finmarchicus, was more common in the CSC than in the stations in the other clusters. The CSC had the highest abundance of copepod nauplii which included Cyclopoida and Metridia longa nauplii. The average abundance and biomass per station was 118,608 ind m⁻² (ranging from 12,747 to 281,834 ind m⁻²) and 0.54 g C m⁻² (ranging from 0.08 to 1.0 g C m⁻²).

The CAC had a total of 57 taxa and an average of 18 taxa per station; however, Pseudocalanus spp. and Oithona similis were much more abundant than any other taxa (Figure 5a). Two stations with a high abundance of Paraeuchaeta glacialis were in the CAC, and the presence of this species inflated the mean biomass per station to 8.9 g C m⁻² (ranging from 3.6 × 10⁻³ to 109.2 g C m⁻²; Figure 5b). However, after removing P. glacialis, the mean biomass per station was only 0.47 g C m⁻², slightly lower than the average biomass in the CSC. Unlike in stations in other clusters, Calanus glacialis, C. hyperboreus and Metrida longa did not constitute a large amount of the biomass. In most stations, the majority of the biomass was contributed by Pseudocalanus spp.

The community of Arctic-adapted and boreal mesozooplankton in the Hudson Bay system changes markedly over the seasons. In spring and early summer Arctic-adapted taxa such as Calanus glacialis and C. hyperboreus are widespread and contribute significantly to the zooplankton biomass. In autumn, when the sea ice has melted several weeks prior and the water column is stratified and oligotrophic, sub-Arctic and boreal taxa, such as Pseudocalanus spp. and Paraeuchaeta spp., become dominant. Although there were similarities in the dominant taxa at most stations (Oithona similis and Microcalanus spp. were the most ubiquitous taxa), there was a marked difference in the distribution of the other taxa among stations. In fact, most of the taxa occurred in fewer than 10% of the stations sampled. The difference in zooplankton assemblages between spring (June of 2018) and late summer and autumn (September–October of 2005) was significantly explained by the day of the year and the timing of the sea ice break-up (measured by the DOW). Surprisingly, salinity did not appear as a major driver of the mesozooplankton assemblage composition, even though the system is hydrologically dominated by freshwater processes. Other local forcings such as bathymetry appeared of secondary importance compared to the day of the year and the sea ice break-up.

While our data set covers different seasons sampled in different years and at different locations, several elements support the hypothesis that the differences in zooplankton composition among clusters reflect a seasonal succession driven by the sea ice break-up and are not simply interannual changes. First, a significant amount of variation in the taxa–environment relationships was explained by the day of the year and by DOW, both variables being essentially independent of the year of sampling. Second, the species assemblages recorded during the ice-melt period and in autumn in different years in the Hudson Bay system are similar to the seasonal succession described in other studies designed to capture the seasonal succession of sub-Arctic seas (Questel et al., 2013; Weydmann et al., 2013; Lischka and Hagen, 2016; Barth-Jensen et al., 2022). Third, while the assemblages observed in stations sampled in July 2010 were more similar to stations sampled in September and October of 2005 than to those sampled in July of 2017 and 2018, the mean ice break-up date of 2010 stations was almost three weeks earlier (May 20) than that of July 2017 and 2018 stations (June 7 and 6, respectively). A seasonal succession tied to ice break-up seems to explain the zooplankton composition observed in 2010.

The day of the year was the variable that described most of the variability in the distribution and abundance of mesozooplankton, with the number of days of open water before sampling a close second. One could argue, however, that DOW is the most important factor in the Hudson Bay system because, first, both day of the year and DOW are highly correlated although not co-linear. Second, early in the season, the sea ice break-up triggers spring primary productivity and the subsequent changes in temperature, salinity and stratification through the summer months. As the snow cover melts, the light field becomes sufficient for pelagic primary producers to start growing under the thinning sea ice and the peak of the phytoplankton bloom is reached only a few weeks after the sea ice break-up (Barbedo et al., 2020; Matthes et al., 2021).

Early sea ice break-up conditions, with high surface salinity and low temperature, were present at stations sampled in June (2018). These stations constituted most of the deep cluster and all of the coastal-spring cluster. The taxa that contributed the most to the total biomass in these stations, Calanus glacialis and C. hyperboreus, usually go through diapause at depth in winter and come back to the surface to start a new production cycle, in synchrony with the sea ice break-up and the phytoplankton bloom that follows (Daase et al., 2013). The high proportions of Copepoda spp. and Calanus spp. nauplii in June and July that were captured despite the 200 µm mesh size, indicates that copepods are reproducing in the Hudson Bay system during the period of ice break-up. In fact, even though they were likely under-sampled, Calanus spp. nauplii were the indicator taxa for coastal stations sampled in June and July. However, this finding does not exclude the possibility that some individuals of C. glacialis and C. hyperboreus remained active throughout the winter. Remaining active might be a beneficial overwintering strategy (Hobbs et al., 2020; Barth-Jensen et al., 2022), as long as the species as a whole succeeds in taking advantage of the spring phytoplankton bloom at the seasonal peak in reproductive activity.

Copepods are not the only mesozooplankton taxa of importance in the Hudson Bay system. As in other ice-covered seas, Cirripedia nauplii are very abundant early in the season in coastal regions (Kosobokova et al., 1998; Weydmann et al., 2013; Stübner et al., 2016). The strong presence of meroplankton has already been recorded in the western part of Hudson Bay by Estrada et al. (2012), but in this study, we found that Cirripedia were highly abundant in coastal stations all over Hudson Bay and Hudson Strait in the coastal-spring cluster. These stations have probably been sampled around the time of spawning for the benthic organisms producing such larvae in 2018 and 2017.

A few weeks after 50% ice break-up, the abundance of large copepods decreased, resulting in a lower biomass except for the stations where Paraeuchaeta glacialis was present. Of the taxa contributing to most of the biomass in coastal stations in spring, Pseudocalanus spp. were the only ones that remained relatively abundant in late summer and autumn. These relatively small copepods are neritic and therefore often found in coastal areas with high freshwater influence throughout the Arctic (Wassmann et al., 2015), meaning that Pseudocalanus spp. can persist in the stratified, chlorophyll-depleted waters, even months after sea ice break-up has occurred. This ability allows Pseudocalanus spp. to thrive in most sub-Arctic, seasonally ice-covered systems like Hudson Bay, but also the Chukchi sea (Questel et al., 2013) and fjords in the European Arctic (Walkusz et al., 2009; Barth-Jensen et al., 2022; Søreide et al., 2022). In several Arctic and sub-Arctic seas, late summer to autumn has been recorded as the reproductive period for Pseudocalanus spp. (Lischka and Hagen, 2005; Kosobokova and Hirche, 2016). Given the relatively low abundance of Pseudocalanus spp. in spring, followed by an increase throughout summer and the presence of Pseudocalanus spp. nauplii, its reproductive period in the Hudson Bay system also likely occurs during the short Arctic summer.

Taxa associated with stratified waters of widely varying salinity and temperature, such as Acartia spp., Decapoda and the amphipods Themisto spp., have been shown to increase in abundance later in the autumn (Rochet and Grainger, 1988; Baker et al., 1993; Harvey et al., 2001; Darnis et al., 2008). Acartia sp. were the indicator taxa of coastal stations sampled in autumn and Acartia longiremis was one of the taxa most correlated to day of the year and DOW. In our study, the abundance of these euryhaline species in September and October might also have been favoured by the anomalously high freshwater influx that occurred during the autumn of 2005 (Lukovich et al., 2021). The dominant meroplankton taxon changed in late summer and autumn from Cirripedia to Bivalvia. Bivalvia reproduce late in the season, and possibly throughout winter. They are known to be more common in autumn in seasonally ice-covered Arctic fjords (Weydmann et al., 2013; Coguiec et al., 2021).

Carnivore and omnivore taxa such as Cnidaria, Chaetognatha, Centropages hamatus and Paraeuchaeta glacialis (Saage et al., 2009; Darnis et al., 2019; Grigor et al., 2020) were abundant in autumn. The presence of these taxa should impose some predation pressure on the Arctic-adapted copepods that dominated the zooplankton abundance and biomass earlier in the season, which could partially account for the low abundances of Calanus glacialis and Calanus hyperboreus by the end of the productive season, also observed in previous studies in the Hudson Bay system (Rochet and Grainger, 1988; Harvey et al., 2001; Estrada et al., 2012).

The physical condition of the surface water could also be impacting the distribution of large copepods. Small copepods such as Oithona similis and Pseudocalanus spp. were the predominant taxa in coastal waters in autumn, where the surface water had a significantly higher temperature and lower salinity than the surface coastal waters in spring and in central Hudson Bay. Oithona similis is known to be very abundant in higher temperatures (Balazy et al., 2021). Furthermore, small copepods are generally less sensitive to temperature differences, while large copepods tend to escape higher surface temperatures via diel vertical migration (Eiane and Ohman, 2004; Darnis and Fortier, 2014). In the shallow coastal areas, the cold refuge provided by deep water is absent, and in Svalbard Calanus copepods were observed closer to the coast in colder years (Balazy et al., 2018). In the Barents Sea, a direct correlation between biomass of C. glacialis and C. hyperboreus and temperature was observed (Aarflot et al., 2018). However, over several decades, the mesozooplankton stock in the Barents Sea stayed stable, despite fluctuations in temperature, because the effects of temperature were counteracted by predation (Dalpadado et al., 2020). This finding demonstrates the importance of evaluating the effects of multiple variables on trends in zooplankton distribution and biomass.

In contrast with other Arctic and sub-Arctic seas, the Hudson Bay system has a higher abundance and biomass of mesozooplankton in the spring than in mid-summer and autumn (Ashjian et al., 2003; Questel et al., 2013; Lischka and Hagen, 2016). The difference in biomass between spring and autumn in the Hudson Bay system could also be due to the disappearance in autumn of the gyre system in the northeast of the bay, which might result in planktonic organisms being flushed out of the bay (Ridenour et al., 2019b). The copepods that have their reproductive period in winter or spring, such as Calanus glacialis and C. hyperboreus, would have more time to be advected and eventually flushed out of Hudson Bay.

Because we compare seasons in different years, the lower zooplankton biomass observed in autumn is possibly a result of interannual differences. However, mesozooplankton assemblages in September and October have also been sampled in 1993 and in the early 2000s in northern and eastern Hudson Bay (Harvey et al., 2001; Estrada et al., 2012). Our data from September and October (2005) match with these previous studies that reported low abundances of lipid-rich copepods in late summer and autumn. In contrast, the relatively higher biomass recorded in June and July of 2017 and 2018 cannot be checked against previous observations from the same region as they are the only ones available. The biomass in some of these stations was nonetheless comparable to the Beaufort Sea shelf (Darnis et al., 2022). While 2017 and 2018 might have been two successive anomalous years, the difference in biomass recorded between years was more likely due to sampling different seasons. Previous studies in which the seasonality of Arctic and sub-Arctic waters was observed through multiple years revealed more pronounced differences among seasons than among years (Arendt et al., 2013; Ershova et al., 2015). Therefore, we hypothesise that seasonal changes in predation, surface temperature and advection could result in the lower abundance of Calanus glacialis, C. hyperboreus and Metrida longa observed in late summer and autumn.

As observed in previous studies (Harvey et al., 2001; Estrada et al., 2012) the deep parts of Hudson Bay and Hudson Strait had consistently higher biomass due to a higher abundance of Calanus glacialis, Metridia longa and C. hyperboreus. Because the relationships between depth and C. glacialis, M. longa and C. hyperboreus were not impacted by sea ice break-up date nor any of the other variables in the RDA, these relationships could involve biotic interactions such as predation. The greater depths in the inner Hudson Bay might allow for effective diel vertical migration of Calanus copepods to avoid predators roaming the shallow depths during the daytime (Aarflot et al., 2019; Aarflot et al., 2020), especially when ice creates a shading effect improving predator avoidance (Fortier et al., 2001). The presence of ice, in combination with other factors such as predator-prey relationships, has an important impact on diel vertical migration of zooplankton (Wallace et al., 2010). Indeed, the visual range of predators is estimated to quadruple after the ice melts (Langbehn and Varpe, 2017). The ice in the centre of the bay differs from the coastal ice as it remains intact for a longer period of time and is thicker and less fragmented (Lukovich et al., 2021). The hypothesis that the distribution of large copepods is partly controlled via top-down processes is also supported by the fact that the planktivorous fish larvae are more abundant in the coastal areas of the bay than in its centre (Schembri et al., 2019; Florko et al., 2021). In particular, Arctic cod larvae (Boreogadus saida) that are known to favour a diet of C. glacialis copepodites (Bouchard and Fortier, 2020) hatch predominantly in coastal areas of the Hudson Bay system and remain close to their hatch location in their first couple of months (Schembri et al., 2021).

The seasonal succession we observed in the Hudson Bay system can be summarized as follows: 1) high biomass of the large Arctic Calanus congeners in the spring when there is still ice algal production as well as phytoplankton production; 2) high abundance of Cirripedia in spring; 3) high abundance of Copepoda spp. nauplii around the time of the peak of the phytoplankton bloom; 4) higher relative abundance of the smaller taxa like Pseudocalanus spp. in late summer and autumn, coupled with a decrease in relative abundance of Calanus spp. and a lower total mesozooplankton biomass; and 5) higher abundances of carnivore and omnivore taxa and Bivalvia meroplanktonic larvae towards the end of the productive season.

Despite the low abundance of zooplankton from mid to late summer, the Hudson Bay system supports an ecosystem with many large predators, including polar bears, seals, and beluga and bowhead whales among others, which demonstrates the importance of the highly productive period during and right after ice break-up in spring. A large variability in environmental conditions is constitutive of Arctic and sub-Arctic marine ecosystems. Moreover, since 1980, the ice break-up date in the Hudson Bay system has been moving forward by an average of 0.58 days per year (Andrews et al., 2018). Unfortunately, the dependence of the base of the Hudson Bay ecosystem on the short sea ice break-up season could increase its vulnerability to events that may push ice conditions outside of acceptable boundaries, potentially resulting in a mismatch between primary and secondary production (Dezutter et al., 2019). However, this unfortunate scenario may materialize only if detrimental environmental conditions were to occur on successive years, as key taxa like Calanus glacialis and C. hyperboreus have adapted complex life cycles spreading over several years, in part to overcome the inherently large variability in sea ice conditions. In fact, in other areas of the Arctic, the loss of sea ice has not resulted in loss of abundance of these copepods yet (Daase et al., 2013; Ershova et al., 2021). As a result, conducting consistent sampling of this sub-Arctic system along the seasons and over the years is critically important to determine whether its observed dynamics are essentially local, or a window open to a possible future for Arctic marine ecosystems.

Data sets within this article have been archived in the Canadian Watershed Information Network (CanWIN). DOI: https://doi.org/10.34992/yfb0-r626.

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

Table S1. Carbon index to convert Copepoda abundance (individuals m⁻²) into biomass (µg C m⁻²)

Table S2. Equations used to convert length of Chaetognatha into biomass and references

Table S3. List of all taxa recorded in the Hudson Bay system

Thank you to Louis Fortier, who supervised SS and TP when this analysis was in its early stages, before his untimely passing. Thanks to the crew and scientists on the CCGS Amundsen, especially Cyril Aubry, for their help with fieldwork. Cyril Aubry also helped with lab work and taught zooplankton species identification to SS and TP. We would also like to thank David Babb for extracting data from regional ice charts. Finally, thanks go to Gérald Darnis for guidance to SS and TP on anything zooplankton.

This research is a contribution to the BaySys project funded by the Natural Science and Engineering Council of Canada (NSERC) and Manitoba Hydro.

The authors have declared that no competing interests exist.

Contributed to conception and design: SS, CB, FM.

Contributed to acquisition of data: SS, TP.

Contributed to analysis and interpretation of data: SS, TP, CB, FM.

Drafted and/or revised the article: SS, CB, TP, FM.

Approved the submitted version for publication: SS, CB, TP, FM.

How to cite this article: Schembri, S, Bouchard, C, Pontbriand, T, Maps, F. 2023. Seasonality of zooplankton communities in the Hudson Bay system. Elementa: Science of the Anthropocene 11(1). DOI: https://doi.org/10.1525/elementa.2022.00113

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

Associate Editor: Maxime Geoffroy, Fisheries and Marine Institute, Memorial University of Newfoundland, St. John’s, Canada

Knowledge Domain: Ocean Science

Part of an Elementa Special Feature: Four Decades of Arctic Climate Change: A Tribute to Louis Fortier