Spatial distribution of epifaunal communities in the Hudson Bay system: Patterns and drivers

The spatial distribution of benthic community structure is predominantly related to variables such as water depth, salinity, substrate type, and food supply (Piepenburg, 2005; Grebmeier et al., 2006a; Cusson et al., 2007; Witman et al., 2008; Roy et al., 2014). In Arctic waters, sea ice cover is an additional variable that influences primary production and thus the efficiency of pelagic-benthic coupling (Piepenburg, 2005; Renaud et al., 2007; Boetius et al., 2013; Roy et al., 2015; Olivier et al., 2020). Most Arctic marine ecosystems are currently responding to climate-induced changes to environmental and ecological variables, such as changing precipitation, river discharge, sea ice cover, and marine biota (Déry et al., 2016; Bring et al., 2017; Osborne et al., 2018; Derksen et al., 2019). In some Arctic regions, such as the Northern Bering Sea, the decline of the sea ice cover has resulted in decline of the clam populations (Grebmeier et al., 2006b; Grebmeier, 2012). Despite a growing number of species inventories being collected by various research projects and programs around the Arctic (Piepenburg et al., 2011; Link et al., 2013; Roy et al., 2014, 2015), baseline knowledge of some Arctic regions is still limited, preventing accurate predictions of how species richness and distribution will respond to climate change (Piepenburg et al., 2011).

Because of its unique attributes, the Hudson Bay System (HBS) nested within the Canadian Arctic (Figure 1) has been identified as one of the most sensitive regions to climate change (Gagnon and Gough, 2005; Tivy et al., 2011; Derksen et al., 2019). The HBS sea ice season has already grown shorter (Andrews et al., 2018) and is projected to continue to shorten, leading to an extended open water (OW) season (Derksen et al., 2019). River inflow is projected to increase up to 50% due to an earlier snowmelt in the surrounding drainage basins and an overall increase in precipitation (Gagnon and Gough, 2005; Bring et al., 2017). These different water inputs lead to variations of nutrient concentrations and salinity within the system that in turn affect biological processes (Déry et al., 2011, 2016). However, little is known about the environmental drivers of benthic communities in the HBS compared to other Arctic regions, such as the Canadian Arctic Archipelago. To date, most of our knowledge on benthic communities comes from relatively old or low spatial resolution diversity data and is based mainly on grab sampling (Atkinson and Wacasey, 1989a, 1989b; Cusson et al., 2007; Kenchington et al., 2011; Piepenburg et al., 2011; Jørgensen et al., 2016; Roy and Gagnon, 2016; Pierrejean et al., 2019; Wei et al., 2019). Most of these studies focused on infaunal organisms, and none related benthic community structure to environmental drivers.

This study aimed to assess the influence of various environmental factors on the structure of epibenthic communities in the study area. The specific objectives of this study were to (i) characterize the epifaunal diversity and density patterns according to environmental variables, (ii) delineate epibenthic communities, and (iii) determine which abiotic factors may drive spatial distribution of the epibenthic communities in the study area. We hypothesized that benthic assemblages will differ along a coast-to-offshore gradient and that salinity, food supply, and sea ice cover will be strong drivers of the epibenthic community structure.

The HBS is composed of four regions: James Bay, Hudson Strait, Hudson Bay (HB), and the Foxe Basin (Figure 1) occupying an area of 1.3 million km². Water masses from the Arctic Ocean enter the HBS from the Canadian Arctic via Fury and Hecla Strait and from Baffin Bay via Hudson Strait (Drinkwater, 1986; Prisenberg, 1986). Within HB, water is rotated cyclonically around the Bay and eventually exported through Hudson Strait (Saucier et al., 2004). However, the major water input for the HBS is river discharge, with around 760 km³ discharged per year from 42 rivers (Déry et al., 2011, 2016). Over half of this river discharge enters the southern and eastern portions of HB, with the largest contributions from La Grande (84.22 km³ y^–1) and the Nelson River (102.70 km³y^–1; Déry et al., 2016). Furthermore, the HBS is covered by a dynamic seasonal ice cover for most of the year (Hochheim and Barber, 2014). Freeze-up progresses from northwest to southeast across the HBS and during recent years has begun in November and formed a complete ice cover by the end of December (Andrews et al., 2018). Within HB, sea ice breakup generally starts in the northwestern and eastern parts of the Bay between May and June and progresses toward the southern region where the last ice typically remains until late July (Andrews et al., 2018; Kirillov et al., 2020). Within the dynamic seasonal ice cover of the HBS, offshore winds generate numerous coastal latent heat polynyas, biologically active areas of OW and thin ice in the dead of winter (Barber and Massom, 2007). The largest polynya is in northwestern HB, but several smaller polynyas are located close to the Nelson River Estuary, the Belcher Islands, Coats and Mansel Islands, and along the coast of Quebec (Barber and Massom, 2007).

Benthic organisms were sampled at 46 stations in HB, Hudson Strait, and Ungava Bay (i.e., study area), between May and July in 2010, 2017, and 2018 (Figure 1). These samples were taken onboard the Canadian scientific icebreaker CCGS Amundsen as part of Arcticnet, the Hudson Bay System study (BaySys) and the Bridging Global Change, Inuit Health and the Transforming Arctic Ocean project (BriGHT). Stations were scattered throughout the HBS in geographically and biologically defined regions with depths ranging from 10 to 322 m (Figure 1; Barber and Massom, 2007; Wilkinson et al., 2009; Kenchington et al., 2011). Epifauna samples were collected with an Agassiz trawl (aperture of 1.5 m and net mesh size of 5 mm), with an average trawling time of 3 min and speed of 1.5 knots, respectively. Four coastal stations were sampled with an epibenthic trawl (aperture of 1 m and net mesh size of 3 mm) with an average trawling time and speed of 3 min and 1.3 knots, respectively. Samples were sieved through 2-mm mesh to retain only macrofauna and megafauna, and identifications were made onboard to the lowest possible taxonomic level. Unidentified taxa were fixed with 4% formaldehyde solution and later identified to the lowest possible taxonomic level under a dissecting microscope. Vertebrates (e.g., Actinopterygii) and planktonic invertebrates (e.g., Chaetognatha and Euphausiacea) collected by the trawl were removed from the analyses. Some taxa were only identified to the phylum level because no complete identification keys exist for HB waters (e.g., Nemertea, Nudibranchia, and Porifera), and hence, taxonomic richness could be underestimated in this study. Taxonomic names were checked and updated using the World Register of Marine Species Editorial Board (2020).

At each sampling station, a conductivity–temperature–depth probe recorded bottom temperature (°C), bottom dissolved oxygen (µM), and bottom salinity (Table 1). Particulate organic carbon (POC; mg m^–3) content measured at the surface of the water column over multiple years (1998–2008) and mean annual surface primary production (PP; mg C m^–2 y^–1) measured over multiple years (2006–2010) were extracted from interpolated environmental data layers generated at the global scale (Global Marine Environment Datasets), as well as in the Eastern Canadian Arctic and Sub-Arctic regions (https://data.mendeley.com/datasets/zmwyjs222s/2), using the package “raster” in R (Table 1; Basher et al., 2018; Beazley et al., 2019; Hijmans et al., 2020). The substratum type at each station was classified into three separate qualitative classes based on either visual observation from trawls or substratum data presented by Henderson (1989) and Pelletier (1986). The three classes of the substrate are “coarse,” referring to stations composed mostly of gravel, sandy gravel, and pebbles (grain size > 2 mm); “mixed,” referring to stations containing particles ranging in size from silt to boulders; and “mud,” referring to stations characterized by fine-grained sediment (grain size < 0.06 mm).

The timing of sea ice breakup and freeze-up and, therefore, the duration of the OW period at each station were extracted from regional ice charts produced weekly by the Canadian Ice Service. Ice charts were produced through expert manual interpretation of remotely sensed imagery, which since 1996 have been based primarily on imagery provided by RADARSAT-1 and RADARSAT-2 (Tivy et al., 2011). In the current study, we defined the sea ice breakup (freeze-up) by the week that the total ice concentration at the study site fell below (surpassed) one-tenth.

Epibenthic characteristics determined for each of the 46 stations were wet biomass (g m^–2), density (ind m^–2), and three biodiversity metrics: taxonomic richness (S, number of taxa), Shannon–Wiener’s diversity index (H’), and Pielou’s evenness index (J ’). The indices H ’ and J ’ were calculated based on biomass data, including colonial taxa. Bryozoa, Porifera, and Cnidaria were excluded from density analysis because whole organisms were not collected by the sampling method. The nonparametric Chao2 index, which represents the expected number of taxa in the study area, was calculated using the “vegan” package (Oksanen et al., 2017). The Chao2 estimator is defined as $Chao2=Sobs+Q122Q2,$ where S_obs is the total number of observed taxa, Q₁ is the number of taxa found at only one station, and Q₂ is the number of taxa found at exactly two stations (Chao, 1987).

Principal component analysis (PCA) was applied to the environmental variables from 32 stations to reduce the dimensionality of the data set using the “FactoMineR” package (Pages, 2004). Fourteen stations were removed from this analysis due to missing data in primary production and bottom dissolved oxygen (Table 1). Environmental variables used in the PCA were depth, bottom salinity, bottom temperature, bottom dissolved oxygen, mean annual surface primary production, surface POC, and duration of the OW season. The first resulting components (PC1, PC2, and PC3), representing a set of environmental variables, were used in linear regressions to model the relationships between community characteristics (density, biomass, and taxonomic richness) and environmental variables. Homogeneity of variance and normality of residuals had been verified using regression diagnostic plots and the Shapiro–Wilk test on residuals.

The list of taxa identified at each station were downgraded to the family level (158 taxa) due to an incomplete set of organisms restraining the identification at the species level for the analysis of the epibenthic community composition. Bray–Curtis dissimilarity was calculated on the fourth-root transformed data for the biomass matrix in order to include colonial taxa. The fourth-root transformation was chosen to balance the effects of high and low biomass taxa (Clarke and Warwick, 2001). To define distinct communities in the study area, the dissimilarity matrix was subjected to a hierarchical cluster analysis using Ward’s (1963) minimum variance method. Community clusters were determined by selecting a distance where stations were fused in well-defined clusters. The geographical distribution of these communities in the study area was then mapped using the “ggplot2” package (Wickham, 2016). Similarity percentage analysis (SIMPER) determined which taxa contributed to the dissimilarity between community clusters based on the biomass matrix (cutoff at 70% similarity; Clarke, 1993).

Epibenthic community characteristics determined at each station were density (ind m^–2), wet biomass (g m^–2), and the alpha (α) and beta (β) diversities. Alpha diversity (α) is determined as the mean number of taxa at a given station. Beta diversity (β), also called turnover diversity, provides an indication of species replacement between habitats or along an environmental gradient and indirectly indicates the habitat diversity. The latter is calculated as the ratio between gamma (γ, total number of taxa in a given community) and alpha (α) diversities (Whittaker, 1960).

The relationship between epibenthic community composition based on wet biomass and the environmental variables was evaluated for 32 stations using a multivariate method of constrained ordination, canonical correspondence analysis (CCA; ter Braak and Verdonschot, 1995). To down-weight rare species, we used a stepwise method that removes the species with occurrence 1, then the species with occurrence 2 and so on. At each step, we compared eigenvalues and the total inertia of the analysis. When a marked decrease of ca 4% was observed in these values, we preserved the previous analysis (Legendre and Legendre, 2012). Removing rare species reduced the total number of taxa from 158 to 101, but this reduction had no impact on the outcome of the analyses. Environmental variables entered into the model were similar to the multiple regressions. To avoid redundancy in the model, we tested linear dependencies among constraints. An analysis of variance for CCA was used to assess the significance of variables and axes. We performed these analyses using the “vegan” package (Oksanen et al., 2017) in R (R Core Team, 2020).

We identified 380 taxa across the 46 stations sampled and 265 to the species level. The Chao2 index reached a value of 539 taxa, exceeding the number of observed taxa. Epifaunal biomass ranged from 0.02 to 45.2 g m^–2, density from 0.11 to 29.5 ind m^–2, taxa richness from 5 to 71 taxa per station, diversity index H ’ from 0.26 to 2.64, and evenness index J ’ from 0.09 to 0.88 (Figure 2).

PC1, 2, and 3 accounted for about 47%, 18%, and 15%, respectively, of the variance in the selected environmental variables, for an approximate combined 80% of the variance (Table 2; Figure S1). PC1 strongly correlated with depth, bottom salinity, bottom dissolved oxygen, surface-water POC, and mean annual primary production (coordinates > 0.6 or < –0.6; Table 2). PC1 reflected an environmental gradient, with low PC1 scores related to shallow depths and low bottom salinity, high POC content, high bottom dissolved oxygen, and low mean annual primary production, whereas high PC1 scores related to deeper waters, high bottom salinity, low POC content, low bottom dissolved oxygen, and high mean primary production. Conversely, PC2 had the strongest correlation with the duration of OW (coordinates > 0.6), while PC3 had the strongest correlation with bottom temperature (Table 2).

Linear models for epibenthic characteristics showed that variables included in the PC1 were significant (P value < .05; Table S1). Biomass, density, and taxonomic richness increased with PC1 scores (Table S1 and Figure 3), which means with deeper water, high bottom salinity, low POC content, low bottom dissolved oxygen, and high mean primary production. The stations located near the Nelson Estuary and in the southern part of HB presented the lowest PC1 scores and the lowest values of biomass, density, and taxonomic richness (Figures 2 and 3). The same observation was noted for Ungava Bay (Figures 2 and 3). Stations located in western HB, northern HB, around the Belcher Islands, and in Hudson Strait presented the highest values of PC1 scores and benthic characteristics (Figures 2 and 3). In contrast, the middle of the Bay was characterized by low values of biomass, density, and taxonomic richness (Figures 2 and 3). Neither PC2 nor PC3 was significant for epibenthic characteristics (Table S1).

Cluster analysis of all 46 stations based on Bray–Curtis dissimilarity of biomass data highlights three communities independent of the sampling year and the sampling gear (Figure 4A). These communities are distributed along a coastal-to-offshore gradient (Figure 4B). Community 1 is located along coastal areas, whereas Community 3 is mainly located in the middle part of the Bay and Community 2 is generally located between the other two (Figure 4B). Across all of the communities, Echinodermata was dominant with different proportions of Echinoidea, Ophiuroidea, and other Echinodermata (Figure 5). A higher proportion of Echinoidea, predominantly sea urchins (e.g., Strongylocentrotus sp.), was observed in Community 1 (23%), whereas Ophiuroidea presented the opposite pattern to Echinoidea with a higher proportion in Community 3 (63%). Community 2 presented an intermediate value of Ophiuroidea and similar proportions of Echinoidea and other Echinodermata (34%, 13%, and 12%, respectively). Arthropoda decreased from Community 1 to Community 3 (19% and 6%, respectively). Mollusca presented a higher proportion in Communities 1 and 3 (22% and 15%, respectively) than in Community 2 (2%). The latter community presented a higher proportion of Porifera (23%). Other taxa, including tunicates (e.g., Ascidiacea), soft corals (e.g., Nephtheidae), and sea anemones (e.g., Actiniaria), were present in notable quantities in Communities 1 (12%), 2 (4%), and 3 (4%). The proportion of Annelida was low across all communities, which reflects reflecting the nature of sampling by trawl. Thirteen taxa were responsible for the differences among the communities, indicating that a particular set of taxa and their respective biomass are discriminant of community dissimilarity (SIMPER analysis, Table S2).

The highest α diversity was observed in Community 2 (28.2 ± 3.3 taxa, n = 11), whereas the lowest value was found in Community 1 (20.5 ± 1.9 taxa, n = 9). The highest value of turnover (β) diversity was found in Community 1 (4.87) compared to the others (3.84 and 4.26 for Communities 2 and 3, respectively). Community 1 was also characterized by a low density and an intermediate biomass (3.07 ± 1.52 ind m^–2 and 6.81 ± 2.58 g m^–2; n = 9). Community 2 had the highest biomass and intermediate density (7.42 ± 1.66 g m^–2 and 5.29 ± 6.61 ind m^–2; n = 11), whereas Community 3 showed the highest density and the lowest biomass (6.65 ± 1.48 ind m^–2 and 4.22 ± 0.89 g m^–2; n = 12).

For all data, the first two axes generated by CCA explained approximately 43% of the taxa biomass–environment relationship, with the first axis explaining the highest variation of 26% (Table 3). The variables coarse substrate, bottom dissolved oxygen, bottom temperature, POC content, and the duration of OW correlated positively with the first axis. Mud substrate, depth, and mean primary production were inversely correlated with these factors and with the first axis (Table 3; Figure 6A). The mixed substrate was located along the second axis (Table 3; Figure 6A). Among the environmental variables considered, only substrate type, salinity, and mean primary production were significantly correlated to the communities (P value < 0.01; Table 3).

The arrangement of the samples on the CCA biplot showed three primary aggregates driven by the substrate type (Figure 6). Stations hosting Community 1 were associated with shallower waters with high values of bottom dissolved oxygen, longer duration of OW, and coarse substrate. These stations were more associated to deposit feeders (e.g., Strongylocentrotidae), filter-feeder bivalves and barnacles (e.g., Pectinidae, Balanidae), and opportunist-predator decapods (e.g., Thoridae and Oregoniidae; Figure 6B). Stations hosting Community 2 were spread along the second axis indicating a high variability of the environmental variables. These stations were found at different depths and temperatures with mixed substrate (Figure 6A). Mixed substrate was linked to filter-feeder sponges, basket stars (e.g., Porifera and Gorgonocephalidae), and soft corals (e.g., Nephtheidae; Figure 6B). Conversely, Community 3 was strongly correlated with high values of mean primary production, low dissolved oxygen, longer duration of ice cover, and deeper waters with mud substrate (Figure 6A). This community was associated with deposit- and filter-suspension feeder bivalves (e.g., Yoldiidae and Astartidae, respectively) and opportunist-predator brittle stars (e.g., Ophiuridae; Figure 6B).

Based on bottom trawl bycatch surveys in the HBS, Jørgensen et al. (2016) estimated a benthic megafauna species richness ranging from 3 to 59 species per station. In our study, the macro-megabenthic taxa richness ranged from 5 to 71 taxa per station. Despite a difference in the gear used to sample river mouths, our results showed no influence of the different types of sampling gear. Moreover, our work increases the observed richness of the macro-megabenthic taxa of the study area (i.e., HB, Hudson Strait, and Ungava Bay) to 380, which is 71% of the taxa expected in the study area, indicating that about one third of the expected species pool is still unrecorded. Unlike previous studies, which have concluded that the HB was a less productive ecosystem than other Arctic regions, our results demonstrate that present values of epibenthic characteristics (i.e., biomass, density, and taxonomic richness) are comparable to other Arctic regions (Grebmeier et al., 2006a; Piepenburg et al., 2011; Roy et al., 2015).

Food supply, salinity gradients, and freshwater discharge are generally considered to be significant environmental drivers for both pelagic and benthic organisms (Remane and Schlieper, 1971; Mayer and Piepenburg, 1996; Piepenburg, 2005; Cusson et al., 2007; Palmer et al., 2011). Univariate characteristics related to biodiversity (density, taxonomic richness, biomass, and Shannon Index) highlighted differences among the 46 stations sampled across the study area. As shown by the PCA, epibenthic biomass, density, and taxonomic richness were driven by the first environmental principal component (PC1), which reflected the following environmental variables: salinity, POC content, depth, dissolved oxygen, and mean annual primary production. Salinity and POC content contributed the most to PC1 and can be considered the main drivers of these benthic community characteristics. Thus, an increase in biomass, density, and taxonomic richness was linked to an increase in bottom salinity and mean annual primary production as well as a decrease in the POC content.

Previously, Pearson and Rosenberg (1978) found that enrichment of organic material (e.g., POC content) generally reduced the number of species and increased both biomass and abundance. Our study showed that high POC content, which we attribute to strong terrestrial contributions, decreased the number of taxa, as well as biomass and density. Stations located in the Nelson River Estuary and along the south coast of HB showed the highest POC content and the lowest values of benthic community characteristics. This area is known to be strongly influenced by river discharge (Prisenberg, 1986; Granskog et al., 2007), which causes a reduction in salinity (to values that range from 28.92 to 32.38) and promotes the deposition of fine riverine particles (Duboc et al., 2017), which alters substrate composition and thus can impact benthic structure significantly (Harrison et al., 2007).

Moreover, the relationships between salinity and primary productivity and benthic biodiversity have already been observed in the Estuary and Gulf of the St. Lawrence River and in the Arctic, especially in the Beaufort Sea (Witman et al., 2008). A drop in salinity is assumed to be responsible for the reduction in benthic richness observed in the Nelson River Estuary and along the south coast of HB. Similarly, stations located in eastern HB showed mainly low values of benthic characteristics. Despite the occurrence of polynyas in the area, increased freshwater runoff and a longer ice season (i.e., late breakup; Kirillov et al., 2020) limit food supply in the area (Granskog et al., 2007; Lapoussière et al., 2009; Ferland et al., 2011; Sibert et al., 2011). Highly biologically productive areas, such as polynyas, generally promote benthic systems in terms of taxa richness, biomass, and secondary production (Ambrose and Renaud, 1995; Link et al., 2011). Despite relatively lower POC content compared to the other geographical areas, stations located near polynyas showed higher diversity, biomass, and density (Figures 2 and 3), hence reflecting strong pelagic-benthic coupling in these areas (Lapoussière et al., 2009; Hochheim et al., 2010).

Central HB is deeper and less productive than both western and eastern HB due to, among other factors, local hydrodynamics (e.g., strong haline stratification in summer) and a longer ice season (Ferland et al., 2011; Kenchington et al., 2011; Sibert et al., 2011). In central HB, stations were associated with lower values of biomass, density, and taxonomic richness. Previous studies have shown that Hudson Strait was more productive than HB during summer and fall (Lapoussière et al., 2009; Ferland et al., 2011). However, our results suggest that HB could be as productive as Hudson Strait in spring because of higher stratification in Hudson Strait confining phytoplankton to the euphotic zone and thus limiting carbon export toward benthic organisms (L Matthes, personal communication). Despite a limited number of stations sampled in Hudson Strait, we observed diversity and biomass similar to HB. Remarkably, stations located in Ungava Bay showed low values of epibenthic characteristics that were similar to those observed near the Nelson River Estuary and southern part of HB. This pattern could be because Ungava Bay is subject to a higher stratification due to large freshwater inputs and intense tidal mixing (up to 17 m) that ultimately reduce food supply by flushing out local production (Markham, 1986; Drinkwater and Jones, 1987).

Our results highlight that epibenthic characteristics are spatially segregated by salinity and food supply (i.e., surface-water POC content and mean annual primary production). Thus, we have shown that epibenthic biomass, density, and diversity increase further off the coast as salinity increases and POC content decreases, until the waters become too deep and these epibenthic variables decline again in central HB.

Salinity, food supply, depth, temperature, and substrate type are known as critical environmental factors that explain the distribution and composition of benthic macrofauna (Mayer and Piepenburg, 1996; Cusson et al., 2007; Bluhm et al., 2009; Palmer et al., 2011; Roy et al., 2014). Data on substrate type are relatively scarce, rather dated and at low spatial resolution in the HBS (e.g., Pelletier, 1986; Henderson, 1989; Misiuk and Aitken, 2020). Despite the difficulty in obtaining recent information on substrate type in the study area, the CCA revealed noticeable changes in diversity and composition of benthic macro-megafauna along an environmental gradient mostly driven by the substrate type, thus also reflecting different feeding types. Epibenthic community structure was divided into three communities based on the substrate (coarse to soft) that characterizes their habitat. Soft substrate mainly characterized the deeper and more productive stations of the study area and the middle of the Bay with low POC content, whereas coarse substrate mainly characterized the shallower and oxygenated stations along the coasts with high POC content.

Community 1 was associated with coarse substrate and dominated by deposit feeders such as echinoderms (e.g., Echinoidea and Strongylocentrotidae) and filter feeders such as bivalves (e.g., Pectinoidae), arthropods (e.g., Balanidae), and ascidians (Ascidiacea). Generally, coarse substrate occurs in areas with strong currents that resuspend sediments and therefore provide food for filter and suspension feeders (Bluhm et al., 2009). This type of benthic community (Community 1), typically attached to rocks and cobbles, was found in coastal areas where river discharge maintained high POC content and high dissolved oxygen, along with strong currents and a longer OW season. Moreover, this community had the lowest density and α diversity values but the highest value for β diversity. Turnover diversity (β) provides an indication of species replacement and, indirectly, of habitat diversity (Whittaker, 1960; Cusson et al., 2007). A high β value indicates a large difference in community composition among the stations in this community. We attribute this difference to the environmental heterogeneity found in this coastal community. These coastal stations encompass a wide range of depths (from 10 to 127 m), salinity (from 24.4 to 33), dissolved oxygen (from 298 to 369 µM), and POC content (from 171 to 774 mg m^–3), reflecting the specific characteristics of each river in the HBS. The strong influence of rivers causes salinity to vary between stations, leading to low values of α diversity (24.2 ± 2.4 species, n = 9) and density (3.07 ± 1.52 ind m^–2, n = 9).

Community 3 was associated with soft substrate (i.e., clay and mud) and dominated by deposit and suspension feeders such as bivalves (e.g., Yoldiidae, Astartidae, and Arcidae), deposit feeders such as holothurians (e.g., Myriotrochus rinkii and Eupyrgus scaber), and brittle stars (e.g., Ophiuridae). Ophiuroids are common in Arctic shelf and slope habitats compared to coastal habitats (Piepenburg, 2000, 2005; Piepenburg et al., 2011). Not surprisingly, brittle stars dominated deeper communities in this study. Moreover, stations characterized by Community 3 were strongly influenced by mean primary production and presented high epifaunal density. However, the prolonged ice cover observed in this area indicates a dominant sympagic system and delayed pelagic primary production (Sibert et al., 2011), suggesting a pulsed food supply. Holothurians are known to exploit fresh phytodetrital pulses on soft sediment (Bluhm et al., 2009; Boetius et al., 2013; Kirillov et al., 2020). Their presence and the presence of deposit and suspension feeders as well as the low biomass observed (4.22 ± 0.89 g m^–2, n = 12) are consistent with a discontinuous food supply in this area. Furthermore, this community was found on average in deeper waters than the coarse substrate community (i.e., Community 1) and presented a high α diversity. As the number of taxa generally increases at depths between 0 and 1,000 m (Levin et al., 2001), the difference in mean depth between Communities 1 and 3 (58 ± 8 m, n = 9, and 160 ± 16 m, n = 12, respectively) could explain the difference in α diversity.

Community 2 was associated with mixed substrate (i.e., particles ranging from coarse to soft sediment) and combined taxa from the other two communities without exhibiting a specific dominant group or class. For example, filters and suspension feeders, such as sponges and soft corals (i.e., Porifera and Nephtheidae), have been associated with poorly sorted sediment (coarse to soft substratum; Hogg et al., 2010). As a consequence, substrate heterogeneity that includes pebbles, cobbles, and/or boulders likely explains why sponges and soft corals were present more in Community 2 than in the other communities. Deposit feeders, however, such as Echinoidea (i.e., Strongylocentrotidae) and Ophiuroidea (e.g., Ophiuridae and Gorgonocephalidae), presented intermediate values of biomass. Moreover, stations supporting these communities were mostly located in deep waters (163 ± 23 m, n = 11) and near recurrent polynyas, such as the large polynya in northwestern HB (Barber and Massom, 2007; Landy et al., 2017). These areas are generally characterized by high food supply and strong pelagic-benthic coupling (Kenchington et al., 2011). Unfortunately, due to the substrate type and sampling method, quantifying sediment pigment concentrations in this study was not possible. Nevertheless, this community showed high α diversity, density, and biomass, highlighting niche diversification (coarse to soft substratum) and an important food supply to this benthic community.

This study presents results from the most recent survey of epibenthic organisms in part of the HBS. We identified 380 epibenthic taxa, representing 71% of the total taxa that are estimated to be present within the study area (i.e., HB, Hudson Strait, and Ungava Bay). The overall analysis of epibenthic characteristics showed that bottom salinity and surface-water POC content were the main environmental drivers. We showed that coastal waters, directly influenced by rivers, harbored the lowest epibenthic density, biomass, and taxonomic richness. These low values were located in shallower depths with low salinity and high POC content. In accordance with previous studies, we also showed that central HB was less productive than the other regions of HB productive than Hudson Strait and we showed that some areas of HB can be as productive as Hudson Strait (i.e., high density, biomass, and taxonomic richness). Further benthic sampling between these two regions would be necessary to confirm this result conclusively.

Based on biomass data, three epibenthic community clusters have been identified and broadly related to the substrate type reflecting food supply proxies (i.e., mean annual primary production and surface-water POC content). More specifically, coarse substrate along the coastlines hosted a higher biomass of filter and suspension feeders, whereas soft substrate in deeper water was mostly associated with deposit and suspension feeders. We showed that the low density and taxonomic richness occurring in the coarse substrate community could be attributed to high POC content and freshwater discharge from rivers. In contrast, the soft substrate community likely received a pulsed food supply because of the longer duration of ice cover and dominant sympagic system. The mixed substrate type did not show dominant taxa or classes and was characterized by large and diverse epibenthic organisms. This last community showed the highest biomass and diversity, which was attributed to a high food supply and strong pelagic-benthic coupling near polynyas.

Benthic organisms respond to natural and anthropogenic changes occurring in their environment, leading to changes in their distribution. Given our results, high POC content of overlying surface water and low salinity bottom water lead to a decrease in biodiversity and change in epibenthic community composition in the study area. Projections toward a longer OW season and increased river discharge as a result of climate change may have a major impact on these epibenthic communities due, in particular, to increases in OW primary production, freshwater inflow, and inputs of terrestrial organic matter from permafrost thaw and forest growth upstream of the HBS.

All data are accessible at the Polar Data Catalogue (https://www.polardata.ca/pdcinput/) and will be made public prior to publication.

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

Figure S1. Principal component analysis (PCA) plots among 32 stations. (A) PCA in terms of environmental variables. (B) PCA in terms of individual stations. OW indicates duration of open water, and PP indicates primary production based on interpolation from satellite data.

Table S1. Results of linear regression tests of density, biomass, and taxonomic richness related to principal component analysis (PCA) axes.

Table S2. Results of similarity of percentage (SIMPER) analysis of the taxa contributing to 70% of the dissimilarity in epibenthic composition based on biomass (g m^–2) of the three identified communities.

This project is part of the Natural Sciences and Engineering Research Council (NSERC)-Manitoba Hydro funded Collaborative Research and Development program known as Hudson Bay System study. Data collection for this research would not be possible without the support and hospitality of the CCGS Amundsen crew during the 2018 field season. D. G. Babb was supported through an NSERC PGS-D, the Canadian Meteorological and Oceanographic Society, and D. Barber’s Canada Research Chair funding. NSERC discovery grant funding (DI) supported research specific to this work, as did the Northern Scientific Training Program. This work is a contribution to the ArcticNet Networks of Centres of Excellence and the Arctic Science Partnership (asp-net.org). We would like to thank the Sentinel North Research Project BriGHT, Amundsen Science, and Québec-Océan for their contribution in sampling time, polar logistic, and scientific equipment. We thank all the participants in the Hudson Bay System study campaign for their contribution to the field work and data collection.

Funding for this research was graciously provided by Manitoba Hydro, the Natural Sciences and Engineering Research Council of Canada, Amundsen Science, and the Canada Research Chairs program.

The authors declare that they have no competing interests.

MP led the design of the study, performed the experiments and taxonomic identification, analyzed the data, prepared figures and tables, and authored and reviewed drafts of this article. DGB performed the sea ice extraction analysis and authored and reviewed drafts of this article. PA and CN led the design of the study, authored and reviewed drafts of this article, and obtained the funding. FM authored and reviewed drafts of this article.

How to cite this article: Pierrejean, M, Babb, DG, Maps, F, Nozais, C, Archambault, P, 2020. Spatial distribution of epifaunal communities in the Hudson Bay System: Patterns and drivers. Elem Sci Anth. 8: 1. DOI: https://doi.org/10.1525/elementa.00044

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

Associate Editor: Laurenz Thomsen, Jacobs University Bremen, Germany

Knowledge Domain: Ocean Science

Part of an Elementa Special Feature: BaySys