The archipelago Novaya Zemlya is located between the Barents and Kara Seas and consists of two major islands (Southern and Northern) and several hundred smaller islands. The coastline of the archipelago is rugged with numerous fjords and inlets. All the investigated fjords are located on the eastern coast of Novaya Zemlya (Figure 1). Hereafter in the text, we refer to the investigated fjords as “bays,” as this is a common geographical name. 12 km-long Blagopoluchiya Bay has a single basin with depths up to 180–200 m separated from the open sea by a 40 m sill. Two rivers originating in the Nally glacier provide a significant runoff here. Sedova Bay is a 20 km-long narrow fjord with a simple bottom relief (Figure 1). Unlike other fjords described in this study, Sedova Bay has no sill at the inlet. There is a 150–200 m deep central basin that meanders smoothly into the Novaya Zemlya Trough. It also has the weakest terrigenous runoff. Tsivolki Bay has a length of 38 km and is the longest of the fjords studied (Figure 1). It has three major basins. The inner semi-isolated basin with a maximum depth of 150 m is separated from the second basin by a 60 m sill. The second and third basins, both 80–100 m deep, are partially separated by sills with numerous islands and connected by the straits of varying depths. Thus, the complex system of islands, shoals, moraines, and depressions, formed due to the interaction of tectonic and glacial processes, prevents the isolation of deep water layers. In contrast, sills only reduce the intensity of water exchange but do not prevent it entirely. The Serp-i-Molot glacier, located at the apex of the fjord, largely determines the characteristics of runoff and sediment accumulation. 24 km-long Oga Bay is very close to Tsivolki Bay and has similar geomorphology (Figure 1). There is a semi-isolated inner basin with a maximum depth of 130 m, where the Goluboy glacier terminates. The second (central) basin is separated from the slope of the Novaya Zemlya Trough by islands and shoals. The depth here mostly does not exceed 60–80 m, but there are several extended depressions deeper than 100 m. 13 km-long Stepovogo Bay is the shallowest fjord studied (Figure 1). It includes two basins separated by a 25 m deep sill: the inner basin with a maximum depth of 60 m and the central one—35–45 m deep. A similar sill isolates the central basin from the outer slope. There is no glacier here, though Stepovogo river provides the fjord with the stream runoff.

All bays are high-arctic and experience pronounced seasonal fluctuations in temperature and salinity. In summer, the surface water temperature in the main parts of the bays and on the outer slope is +4–5.5°C, but in the areas adjacent to the glaciers, it drops to +1.8°C. Surface salinity varies depending on the volume of river/glacier runoff: from 21 to 23 in the inner periglacial areas to 31.2–32.5 in the outer parts of the bays. A transition layer with sharp changes in salinity and temperature is observed at a depth of 10–40 m, depending on the season and runoff rate. Below this zone, the temperature and salinity distribution are similar throughout the bays. Moreover, since the sills between the basins are located deeper than the seasonal thermocline, they do not prevent free water exchange between the open parts and inner basins (Stepanova and Nedospasov, 2017).

The volume of glacial runoff, the intensity of snowfall, prevailing winds and wave action influence the characteristics of suspended matter fluxes. Water turbidity in the bays varies greatly both seasonally (during the year) and between years. Its values are determined by glacial runoff (as in Oga and Tsivolki Bays) as well as river runoff, which in turn may be glacial-derived (as in Blagopoluchiya and Sedova Bays) or mixed (as in Stepovogo Bay). In the inner periglacial parts of Oga and Tsivolki Bays, turbidity reaches its maximum (Chava et al., 2017; Udalov et al., 2019).

Bottom sediments in the inner parts of the bays and on the sills to a depth of 45–50 m frequently consist of debris (gravel, pebbles, shales) silted to varying degrees. In the inner and central parts of the bays from 50 m to maximum depths, sediments are represented by silts (Udalov et al., 2016, 2018, 2019, 2020; Chava et al., 2017).

Sampling was carried out during five expeditions of the RV “Akademik Mstislav Keldysh” and the RV “Professor Shtokman” in 2013–2016. Fifty stations were collected in five bays (Blagopoluchiya, Sedova, Tsivolki, Oga, Stepovogo) and on the adjacent slope. In all bays, stations were taken along the bay axis from the inner part to the outer slope of the Novaya Zemlya Trough (Figure 1). The “Okean” grab or the Van Veen grab (0.1 m²) were used for sampling (3–5 replicates per station) (Eleftheriou and McIntyre, 2005). Samples were washed through a 0.5 mm mesh sieve and then fixed with a buffered solution of 6% formalin. In the laboratory, animals were sorted, identified to the lowest possible taxonomic level, counted, weighed (wet weight), and preserved in 70% ethanol. Molluscs were weighed with shells, and polychaetes were weighed without tubes. A detailed description of the benthic fauna of each bay is available in earlier publications (Udalov et al., 2016, 2018, 2019, 2020; Chava et al., 2017). Complete station data are given in Supplementary Table 1.

Vertical profiles of the water column were obtained using SBE 911 Plus CTD. Eighty-seven CTD measurements of temperature, conductivity (salinity), dissolved oxygen and turbidity were carried out. We analyzed the data measured over the entire area of the bays (not only at benthic stations), which allowed an adequate estimation of turbidity. The volume ratio of the gravel fraction in the sample was used as a descriptive sediment parameter. According to Wentworth grain size classification, the gravel fraction included pebble and cobble (Eleftheriou and McIntyre, 2005).

The following hydrological data were used in the statistical analysis: near-bottom salinity, temperature and dissolved oxygen. Since we had no direct measurements with sediment traps, we used mean turbidity in the water column as an approximation of terrigenous (i.e., glacial) flux and sedimentation rate to the bottom:

Further, we analyzed turbidity distribution as a function of distance from the apex of the bay. Detailed data on hydrological characteristics are given in Supplementary Table 1.

In all subsequent analyses, the species respiration rate R was used as a measure of abundance, estimated as:

where N_i—the abundance of a species, B_i—biomass, and k_i—a taxon-specific coefficient. Then, the relative respiration (metabolic) rate for a species i was calculated as a proportion of the total: r_i = R_i/(∑R_i) (Kucheruk and Savilova, 1985; Azovsky et al., 2000). The relative respiration rate characterizes the relative species contribution to community production and metabolism and shows only the percentage of each taxa. By combining both abundance and biomass of each species, this measure provides a balanced contribution of small but abundant species and large ones with low abundance but high biomass (Azovsky et al., 2000; Vedenin et al., 2015). Large species represented by single individuals with extremely high biomass (several sponges, ascidians, echinoderms) were excluded from the analysis as their distribution cannot be adequately estimated using small grabs. Prior to multivariate analyses, data were fourth root transformed to reduce the role of dominant species (Clarke et al., 2008). Finally, the similarity between stations was estimated using Bray–Curtis similarity index based on species respiration rates data.

Non-metric multidimensional scaling (nMDS) and hierarchical clustering (UPGMA) were performed to identify key patterns in community structure. First, a similarity profile test (SIMPROF type 1) was applied to detect clusters with non-random structure (p < 0.05) (Clarke et al., 2008). The resulting clusters were analyzed using similarity percentage routine (SIMPER) to determine which species were responsible for the intra-group similarities and between-group dissimilarities (Clarke, 1993). In addition, number of taxa, total abundance, biomass and univariate diversity measures [Shannon index (H’loge), Margalef species richness (d), Hurlbert rarefaction (ES100), Pielou’s evenness (J)] were calculated for each cluster (McCune et al., 2002). A shade plot was constructed to visualize the differences in species structure between the clusters (respiration rate data). According to the identified clusters, the key species that provided a high proportion of within-group similarity in station groups (by SIMPER) were selected for analysis. Taxa were grouped in clusters using UPGMA algorithm based on the index of association, coherent species sets were identified by the Type 3 SIMPROF tests. Permutational analysis of variance (PERMANOVA) on the Bray–Curtis similarity matrix (respiration rate data) was used to confirm the differences in species structure between the groups of clusters, identified in the shade plot. Next, we examined the differences in the quantitative characteristics of these groups using the non-parametric Kruskal-Wallis test.

Relationships between the environmental variables and macrofauna distribution patterns were analyzed using DISTLM (Distance-based linear model) with forward selection of predictor variables (McArdle and Anderson, 2001; Anderson et al., 2008). Four independent parameters were used in the analysis: depth (DEPTH), temperature (T), oxygen concentration (O₂), mean turbidity (TURBIDITY). Salinity was excluded because of a strong correlation with temperature. The following variables were used as indicators: SEDIMENT (two levels: presence/absence of stone fraction) and BAY (categorical variable with five levels which corresponded to five fjords investigated). In addition, we performed DISTLM for each fjord separately to identify any patterns related to environmental factors within fjords that the joint analysis could potentially mask.

Statistical analyses were performed in PRIMER V7 with the PERMANOVA+ add-on package (Clarke and Gorley, 2015) and PAST 3 (Hammer et al., 2003) software.

The benthic zone deeper than 40–50 m was exposed to the same temperature and salinity throughout the entire sampled area in all the bays (Table 1). The distribution of dissolved oxygen was similar to the distribution of the main hydrophysical parameters (temperature and salinity). The lowest values were registered in the deepest parts of Oga and Tsivolki Bays (>100 m), especially in the periglacial basins, although even in these areas the oxygen content near the bottom did not drop below 6.06 ml/l (72.8%). At depths of up to 50 m, the average saturation level of dissolved oxygen in the near-bottom waters was 91%.

Analysis of turbidity distribution as a function of distance from the apex of the bay allowed us to distinguish three types of fjords (Figure 2). The first type includes the bays with a glacier terminus (Oga and Tsivolki). There is only one source of the suspended matter here (the glacier), turbidity values are extremely high (up to 24.5 FTU—the upper limit of measurement rate) and they decrease with distance from the source as power function (R² = 0.854 and R² = 0.823 for Oga and Tsivolki, respectively). The second type comprises Sedova and Stepovogo Bays with depleted runoff and lower turbidity (TURB = 0.46 FTU). The negative correlation between turbidity levels and distance from the inner part of the bay is much weaker than in the first type (R² = 0.386 and R² = 0.485, respectively). Blagopoluchiya Bay is of the third type and is characterized by high turbidity values throughout the bay (R² = 0.069). This is most likely due to the presence of several (though weaker than in Oga and Tsivolki Bays) sources of runoff from the rivers of glacial origin, which “feed” different parts of the bay. In the outer parts of all the bays and on the slope of Novaya Zemlya Trough, turbidity values become more or less equal (Figure 2).

The distribution of bottom sediments varied among the different bays. In general, the substrate at most stations was mostly silt; hard sediments were mainly found in shallow waters, especially in Stepovogo Bay (Table 1).

In total, we identified 246 different taxa. Among them 77 species were rare, i.e., represented by only one or two individuals. The share of rare species was proportional to the total number of species encountered at the station or in the bay. The four most common species (bivalves Mendicula ferruginosa, Yoldiella solidula and polychaetes Scoletoma fragilis, Micronephthys minuta) contributed 36% to the total benthos abundance. Two species (brittle star Ophiopleura borealis, sipuncula Golfingia margaritacea) accounted for 31% of total biomass.

The hierarchical clustering revealed 11 distinct non-random clusters (SIMPROF test, p < 0.05) (Figure 3). Two stations located on the stony-gravel bottom were not included in any cluster: 125–45 (Blagopoluchiya Bay) and 5259 (Stepovogo Bay). In most cases, one cluster included stations from different bays, except for clusters 3 and 5 (Stepovogo Bay), 7 (Blagopoluchiya Bay), and 8 (Sedova Bay) (Figure 3). Cluster 1 and cluster 2 included the stations with severely depleted fauna from the inner parts of the bays (terminal stations), located in a wide range of depths (from 20 to 150 m) on the muddy bottom with increased sedimentation rate. Bivalve Portlandia arctica dominated at three terminal stations assigned to cluster 1 (Table 2 and Figure 4). These stations were located either on the slope of periglacial basins (Oga and Tsivolki Bays) or in the area of intense runoff of the glacial river (Blagopoluchiya Bay). Isopod crustacean Saduria sabini was the dominant species in cluster 2, which included two stations located in the deepest parts of inner basins at the depths of 53 m (Stepovogo Bay) and 143 m (Tsivolki Bay). All five stations were characterized by low species richness, total abundance and biomass and had the lowest values of diversity indices (Table 2). Shallow-water stations (up to 40–50 m depth), located mainly in Stepovogo Bay in areas with high debris content, were assigned to three different clusters (clusters 3, 4, and 5), although 32 species were common to all three clusters, representing 82% of the total abundance. Cluster 3 included five stations from the inner part of Stepovogo Bay. Sipuncula Golfingia margaritacea dominated here, followed by bivalve Bathyarca glacialis and polychaete Scoletoma fragilis (Table 2 and Figure 4). Cluster 4 (four stations) combined stations from the inner and outer parts of Stepovogo Bay and the sill areas of Stepovogo and Blagopoluchiya Bays. Although its species composition was similar to cluster 3 (Figure 4), the dominant species were two bivalves, Astarte elliptica and Astarte borealis, which mainly were absent at the stations belonging to cluster 3 (Table 2 and Figure 4). Cluster 5 included only two stations located on the stony-gravel bottom, where polyplacophoran mollusc Tonicella marmorea and polychaete Pista maculata dominated (Table 2 and Figure 4).

Cluster 6 included nine stations from the inner basins of Blagopoluchiya, Tsivolki and Oga Bays located on muddy bottom deeper than 45–50 m. Although quantitative characteristics (total abundance and biomass) were similar to those of other groups (except for cluster 1 and cluster 2), diversity was significantly lower (Table 2 and Figure 4). The dominant taxa were bivalves Portlandia arctica, Ennucula tenuis, and brittle star Ophiopleura borealis (Table 2). Each of clusters 7 and 8 included stations from only one bay. Cluster 7 consisted of three stations from the sill slope in Blagopoluchiya Bay (60–130 m). Among the dominant species were sipuncula Golfingia margaritacea and bivalve Bathyarca glacialis (Table 2 and Figure 4). Cluster 8 included four stations from the deep central part of Sedova Bay (130–200 m). The most abundant species were brittle star Ophiopleura borealis and polychaetes Scoletoma fragilis and Nothria hyperborea (Table 2 and Figure 4). Seven stations from the outer parts of Sedova, Oga, Tsivolki and Blagopoluchiya Bays formed cluster 9, where bivalves Astarte crenata, Ennucula tenuis, Yoldiella solidula, and Yoldiella lenticula dominated. These stations were located on the muddy bottom at a depth of more than 50 m. Clusters 10 and 11 were assigned to the outer slope of the bays. Cluster 10 included five stations dominated by Astarte crenata bivalves, located on the slope of the Novaya Zemlya Trough in front of Tsivolki, Stepovogo and Sedova Bays (50–125 m deep). These stations were characterized by the maximum species diversity (Table 2). Four stations located in the lower part of the outer slope of Tsivolki and Stepovogo Bays (75–220 m) formed cluster 11. The most dominant taxa were Nephasoma spp. sipunculids, Scoletoma fragilis polychaetes and Astarte crenata bivalves. Diversity indices, the mean number of species per station and ES100 were lower compared to other clusters and similar to those described for cluster 6 (stations from the deep inner basins). Total abundance (563 ind/m²) was also 2–4 times lower (Table 2).

Using SIMPER-analysis, we evaluated the differences in the species composition of the identified communities to answer whether each community is assembled from different species or whether they are similar in species composition but differ only in dominants. It was found that the same species contributed most to the intra-group similarity and that the difference between the groups was due to the dominance of a small number of species (Figure 4). For example, Scoletoma fragilis polychaetes contribute significantly to the intra-group similarity in clusters 3, 6, 7, 8, 9, 10, and 11 and Ennucula tenuis bivalves in clusters 6, 7, 8, and 9 (A detailed SIMPER-analysis is given in Supplementary Table 2).

Thus, there were several clearly distinguished groups of species in all of the bays studied, typical for particular habitats (Figure 4). First, there was a group of species characteristic only for the periglacial areas of the bays and for the inner basins (A), such as Saduria sabini isopods, Portlandia arctica bivalves and Cossura longocirrata polychaetes. Another important group (B) included species found almost everywhere, in all the bays studied: polychaetes Aricidea hartmanae, Scoletoma fragilis, Micronephthys minuta, Cirratilidae gen. sp., bivalve Thyasira dunbari and tanaid Acanthophoreus gracilis. The following group of species (C) (brittle star Ophiopleura borealis, polychaetes Artacama proboscidea, Scoloplos armiger, Aglaophamus malmgreni, and bivalves Mendicula ferruginosa, Yoldiella lenticula, Yoldiella solidula) was also found everywhere, except for the stony substrates in Stepovogo and Blagopoluchiya Bays. The set of species occurring mainly on the deep outer slope of the bays (D) was distinguished separately. The next group (E) included species present in mixed sediments both in Stepovogo Bay and on the outer slope of the bays. In addition, the group confined to stony substrates in Stepovogo and Blagopoluchiya Bays (F) was identified. This set of species, living on the hard substrate, was also found on the outer slope of the bays, where gravel and boulders are noted, but in significantly smaller amounts.

Results of shade plot showed that five fjord zones differed in the structure of benthic communities: I zone represented the depleted inner periglacial areas (clusters 1–2), II—the subtidal belt of hard substrates at shallow depths (clusters 3–5), III—the deep inner semi-isolated basin (cluster 6), IV—the outer non-isolated basins and upper slope (clusters 7–10), and V- the lower slope (cluster 11) (Figure 4). These differences appeared statistically significant in the PERMANOVA analysis (p(MC) = 0.0001). The distinguished zones differed both in community structure and in basic integral characteristics (Figure 5 and Table 3). The lowest abundance, biomass, and diversity values were in the depleted inner periglacial areas (I) and in the lower slope (V). Relatively low diversity was also in the deep inner basin (III). The highest diversity was observed in the zone with the presence of hard substrates (II), the outer non-isolated basins and upper slope (IV). Thus, our data reliably show the association of certain benthic communities to different parts of the fjord (such as sill, periglacial basin, outer slope, etc.).

DISTLM model demonstrated that only half of the total biotic variance (42.1%) was explained by the selected set of environmental variables (Table 4). The forward selection procedure indicated that a higher percentage of explanation was observed for factor BAY (25.2%), followed by SEDIMENT (6.8%) and TURBIDITY (5%). The effects of these factors were significant (p < 0.001). The influence of DEPTH, T, O₂ was not significant.

If we consider each of the fjords separately, we found two different models explaining macrofauna distribution patterns. In the fjords with intensive glacial runoff (Oga, Tsivolki and Blagopoluchiya) TURBIDITY explained a significant percentage of variations in community structure (55.0–18.9%). In contrast, in the fjords with low runoff (Stepovogo and Sedova), its influence was weaker, and TURBIDITY was excluded from the model by the forward selection procedure. Thus, the other factors (DEPTH and SEDIMENT) played an essential role in structuring the benthic communities in those fjords. The results of DISTLM are summarized in Supplementary Table 3.

When testing the general DISTLM model for three fjords with intensive runoff (Blagopoluchiya, Tsivolki and Oga) (with the exclusion of bay-specific clusters), the contribution and the order of factors changed with the most considerable contribution of TURBIDITY (p = 0.0001). As a result, the effect of BAY became weaker (Table 5).

Fjords are semi-enclosed marine inlets that remain under strong terrestrial influences and are considered to be strongly dependent on regional species pools of neighboring open shelf seas (Pearson, 1980). Traditionally, fjord communities are perceived as range extensions of shelf communities (Buhl-Mortensen and Høisæter, 1993; Josefson and Hansen, 2004) with a subset of the regional species pool and reduced diversity. The diversity clines in fjords are often explained by either the barrier hypothesis, related to possible colonization barriers posed by fjord geomorphological features (as sills) or simply the distance to the main species pool located on the shelf and/or the disturbance hypothesis, i.e., environmental deterioration of the fjord habitats (Buhl-Mortensen and Høisæter, 1993; Buhl-Mortensen, 1996; Włodarska-Kowalczuk et al., 2012). The barrier hypothesis assumes that shelf seas serve as local species pools for fjord biocenosis.

In the fjords of Novaya Zemlya, all communities consist of species widely distributed in the Kara Sea in a similar depth range (Jorgensen et al., 1999; Denisenko et al., 2003; Vedenin et al., 2015; Azovsky and Kokarev, 2019). In our study, the percentage of rare species (singletons or doubletons found only in one fjord) was proportional to the total number of species encountered at the station or in the fjord and hence was explained by statistical reasons rather than by the specificity of the fauna of a particular fjord. A large number of studies have shown the same pattern. A similar set of species has been shown for fjords, bays and adjacent offshore areas in the Canada Basin (Aitken and Fournier, 1993); northern Norway (Holte and Gulliksen, 1998; Oug and Høiśter, 2000); Greenland (Thorson, 1934; Sejr et al., 2000). In a comparative analysis of the different bays and fjords of the Barents Sea, it has been shown that the main pool of species in the bays is represented by common, widespread boreal and arcto-boreal species. At the same time, although these species are quite common, they are rarely dominant in the open sea shelf communities and usually represent a group of subdominant species (Britayev et al., 2010).

Contrary to this, in several cases, the specificity of fjord fauna was demonstrated. The comparative analysis of west Spitsbergen fjords and offshore Barents Sea species lists showed that 30% of all species recorded occurred only in fjords (Włodarska-Kowalczuk et al., 2012). Basin-specific communities were also described for the deepest fjord of northern Norway—Tysfjord multi-basin system (Jorda Molina et al., 2019). The general explanation is that the silted sub-arctic or arctic Barents Sea fjords are a refuge for Arctic benthic species, which appear to have local populations in the inner parts of the fjords (Nordgaard, 1905; Soot-Ryen, 1951). The relative isolation of these basins facilitated this process, maintaining low temperatures due to reduced water exchange with an open sea. The same factor is responsible for the development of an arctic community dominated by Portlandia arctica in the deep part of the White Sea (Naumov and Fedyakov, 2000).

In the fjords of Novaya Zemlya, the depth of the sills and vertical stratification of the water masses prevents the isolation of the inner parts of the fjords. The seasonal thermocline is located above the sill, and water exchange between the open sea and the inner fjord basins is intense (Stepanova and Nedospasov, 2017). As a result, the inner parts of the fjords are quite similar to the open sea in terms of their hydrological characteristics (salinity, temperature, oxygen), which prevents the appearance of isolated faunal complexes.

Besides, fjord communities of the Novaya Zemlya include the same species in various combinations. Therefore, the difference between communities is usually determined not by the presence of initial community-specific groups of species but by the predominance of a particular species (or several species) from the general pool of species (Figure 4). This leads to the formation of local specific communities, similar in species composition but differing in the number and order of dominant species.

However, we can see differences between the investigated fjords in terms of benthic composition, as indicated by the significance of the BAY factor in the overall DistLM model (Table 4). First of all, the influence of the BAY factor can be explained by the presence of bay-specific environmental zones inhabited by the bay-specific communities. In particular, cluster 7 was confined to the sill in Blagopoluchiya Bay, and cluster 8 to the deep central part in the non-isolated Sedova Bay, opened to the Kara Sea slope. At the same time, most stations in clusters 3–5 were confined to the shallow Stepovogo Bay (Figure 3) due to insufficient sampling on gravel substrates at the depths of 30–50 m in other bays (caused by the risks of ship operations). Yet, similar communities were identified in Blagopoluchiya Bay (stations 125–45 and 125–38, Figures 3, 4). In contrast, when considering the same geomorphological zones of different fjords (e.g., the deep inner periglacial basins), we found the same communities with the similar species pool and integral characteristics (Figures 3, 4 and Table 2).

The second reason is the differences between the fjords in the level and intensity of runoff (Figure 2). This factor defines the set of benthic communities and their distribution in the fjords with significant glacial runoff (Oga, Tsivolki and Blagopoluchiya) (Figure 4). We see the presence of communities belonging to clusters 1, 6, and 9 successively replacing each other along the fjord axis (Figures 3, 4). In this case, we can distinguish clear benthic zonation from the depleted inner periglacial areas with the highest sedimentation rate toward the outer parts of the bay and the lower slope (Figure 6). Within this scheme, we see five main zones in which local specific variants of communities can be implemented.

Two zones (Zone II and IV) are the most diverse and abundant in quantitative characteristics of macrofauna (Table 2). The first of them is the sublittoral belt at 20–50 m depths (Zone II). This zone is characterized by environmental heterogeneity related to seasonal variation of temperature and salinity, intense hydrodynamics, the presence of macroalgae and various sources of organic matter, high content of debris material. Benthic communities in this zone have a patchy distribution and vary in dominants. But all of them are characterized by the presence of a single pool of species preferring hard substrates (Figure 4), which is typical for the eastern coast of Novaya Zemlya (Uschakow, 1931). The outer and central parts of the bays and the upper part of the slope (Zone IV) are the most diverse, considering the soft-sediment fauna. Benthic communities in this zone are similar to the benthic communities of the southwestern part of the Kara Sea (Jorgensen et al., 1999; Denisenko et al., 2003). The successive impoverishment of benthic communities forced by an increase of terrestrial and glacial runoff was observed along the axis of the fjord toward the inner part. As a result, fjord-specific communities formed in Zones III and I (Figure 6). Toward the open sea (Zone V), we also detected community change and the decrease of quantitative and diversity characteristics related to the vertical zonality of the slope of the Novaya Zemlya Trough (Galkin et al., 2010). In the fjords without glacial terminus, we can find community types similar to Zones II and IV communities. However, their distribution is more heterogeneous and depends on the depth and sediment structure (Table 5).

The presence of strong gradients related to sedimentation rates has been noted in previous studies of benthic communities in glacial fjords (Holte and Gulliksen, 1998; Włodarska-Kowalczuk and Pearson, 2004; Włodarska-Kowalczuk et al., 2005; Krishnapriya et al., 2019). Depletion of the innermost fjord basins of both biomass and number of species has also been commonly reported (Görlich et al., 1987; Włodarska-Kowalczuk et al., 2005, 2012), and this is clearly evident in our data as well. The presence of intense glacial runoff leads to the development of communities 1, 2, and 6 (in Oga, Tsivolki, Blagopoluchiya Bays). At the same time, the main feature—glacial runoff—may be provided either by direct terminal glacier outflow into the bay (Tsivolki Bay) or by runoff from a glacier located at a distance in the inland part of the island (Blagopoluchiya Bay).

The key community developing under such conditions is a depleted community (both quantitatively and in species diversity) dominated by Portlandia arctica—Ophiopleura borealis—Ennucula tenuis. As sedimentation rates increase, the number of species drops to 8–9, and the dominance of P. arctica reaches 75%. P. arctica community is widely represented in the innermost parts of many Arctic bays and fjords (Spark, 1933; Thorson, 1933, 1934; Dale et al., 1989; Syvitski et al., 1989; Rozycky, 1992; Aitken and Fournier, 1993; Aitken and Gilbert, 1995; Włodarska et al., 1996; Holte and Gulliksen, 1998; Sejr et al., 2000; Renaud et al., 2007). This eurybiontic species can withstand considerable warming (up to +5–6°C and even +12°C) and sea salinity (Filatova, 1951). In the White Sea, this species may survive in warm shallow waters (10–30 m) in semi-isolated bays (Naumov and Fedyakov, 2000). In contrast to the White Sea or coastal areas of the southern Kara Sea (Ob-Yenisei shallows), where P. arctica inhabits deeper adjacent seabed areas with constant subzero temperatures, in the fjords of the eastern coast of Novaya Zemlya P. arctica forms completely isolated populations. In Blagopoluchiya Bay, this species inhabits the inner basin, reaching depths of 120 m, where it dominates together with E. tenuis and O. borealis. In Tsivolki Bay, P. arctica is found on silts in the central basin at a depth of about 60 m, where it is second in biomass (8%) after E. tenuis (74%). At the same time, according to our data, P. arctica is absent on the outer slope of Novaya Zemlya. It is likely that isolated parts of the species’ areal may be relict post-glacial enclaves (Filatova, 1951), preserved in sufficiently large basins. In particular, in Blagopoluchiya Bay, the presence of a large isolated population in the uppermost basin, where water temperatures are constantly negative, allows this species to inhabit the shallower parts of the bay successfully. However, P. arctica is absent in Stepovogo Bay, which is probably due to a lack of suitable biotopes, although the inner basin has both permanent subzero temperatures and silted sediment.

Several fundamental issues can be highlighted when analyzing the distribution of the P. arctica community across the Arctic region. The first one is that this community is absent in almost all of the studied fjords of the west coast of Spitsbergen, including Kongsfjorden (Włodarska et al., 1996; Kendall et al., 2003; Włodarska-Kowalczuk and Pearson, 2004; Wêsławski et al., 2012). Among more than ten investigated fjords, this species is found only in the inner glacial basins of Van Mijenfjord (Holte and Gulliksen, 1998; Renaud et al., 2007) and Hornzund (Rozycky, 1992). The authors attribute this to the warm-water character of the Western Svalbard fjords, where P. arctica communities survive exclusively in the coldest parts of the fjords.

It is widely agreed that near-glacier waters are oligotrophic due to the inorganic sedimentation and low content of available organic matter (Görlich et al., 1987). Most of the suspended matter in the bays with the presence of outlet glaciers is represented by a fine mineral fraction. Non-selective surface-deposit-feeding species and filter feeders are therefore at a disadvantage, as they cannot handle the intense flow of mineral particles and the advantage goes to organisms with the ability to select food particles. As an example, changes in the species composition and quantitative characteristics of bivalves with distance from the glacier are commonly cited (Syvitski et al., 1989; Aitken and Gilbert, 1995). According to this pattern, the mobile protobranch bivalves, which have the ability to select food particles with labial palps (Ward and Shumway, 2004), can survive exceptionally high rates of sedimentation and are the first to inhabit the near-glacial areas of the fjord (Syvitski et al., 1989). Apart from P. arctica, it is also E. tenuis, which is widespread in the fjords of Greenland and Baffin Land (Thorson, 1933; Syvitski et al., 1989; Aitken and Fournier, 1993), northern Norway (Oug and Høiśter, 2000), Franz Josef Land (Dahle et al., 2009a) and is noted among the eight most frequent species in the fjords of Spitsbergen (Wêsławski et al., 2012).

Another impoverished community (cluster 2), which we also observed in our study, was characterized with very low species richness and dominance of a large mobile predator fauna—Saduria sabini, which is quite tolerant to changes in salinity and temperature (Percy, 1985). However, this community does not appear to be specific to glacial areas. It develops in unstable conditions, which may be related to frequent sediment slides and gravity flows, resulting in the formation of unconsolidated, easily eroded sediments, where there are no suitable habitats even for protobranchs tolerant enough to increased sedimentation. This is the case of the communities in the first basin of Stepovogo Bay in particular.

As for the polychaetes present in glacial bays and fjords, deposit feeders, carnivores/omnivores and predators play a significant role. Mass polychaete species of Novaya Zemlya bays are burrowing deposit-feeding family Cirratulidae, Cossura longocirrata, and carnivore Scoletoma fragilis, Aglaophamus malmgreni (Jumars et al., 2015). They are also common in the fjords of the Svalbard and Northern Norway (Gromisz, 1992; Holte et al., 1996; Włodarska et al., 1996; Włodarska-Kowalczuk et al., 1998; Włodarska-Kowalczuk and Pearson, 2004). Furthermore, the bays of Novaya Zemlya are inhabited by species indifferent to disturbance and widespread in the coastal and estuarine areas of the Kara Sea and the Laptev Sea, where an increased runoff of organic matter has been noted. These are deposit-feeders Aricidea hartmanae, Scoloplos armiger, Capitella spp., and predator Micronephthys minuta (Averintsev, 1990; Zhirkov and Paraketsova, 1996; Vedenin et al., 2015; Gagaev, 2021).

Thus, small-bodied surface-deposit feeding species with low biomass that dominate in near-glacial habitats are better adapted to resist disturbance, which is consistent with the classical successional models for organic pollution (Pearson and Rosenberg, 1978) and physical disturbance of sediments as a result of dredging (Rhoads et al., 1978).

The state of the Kara Sea benthic ecosystem was considered relatively stable until recently (Kozlovskiy et al., 2011; Vedenin et al., 2015; Gerasimova et al., 2021). However, in 2012, the snow crab Chionoecetes opilio was found in the south-western part of the Kara Sea (Zimina, 2014) and in 2013 in the south- and northeastern parts of the sea (Sokolov et al., 2016). Snow crabs are benthic omnivores, and their diets typically comprise a wide variety of prey types, including macrobenthos, fish and fishery discards, and other crustaceans (Squires and Dawe, 2003). Thus, decapods are largely involved in predator-prey interactions and are known to compete with other species over food and grounds (Boudreau and Worm, 2012; Gebruk et al., 2021). A number of studies showed that predation could shape benthic communities and production processes (Grosholz et al., 2000; de Rivera et al., 2011), and significant shifts in the ecosystem’s trophic structure are described due to increased abundance of benthic decapods (Baum and Worm, 2009).

The current work describes the state of the benthic ecosystem before the mass development of the snow crab, so it can be regarded as a starting point for further research of invasive predator impact on the structure of benthic communities. Before the invasion of C. opilio, large invertebrate predators were almost absent in the bottom communities of the Kara Sea, and benthic fish capable of consuming C. opilio are less common than in the neighboring Barents Sea (Dolgov and Benzik, 2016). The introduction of C. opilio into the Kara Sea provides a rare opportunity to trace the development of an invasive alien species and the response to it of an ecosystem previously unaffected.

The question of whether food webs are bottom-up (resource) or top-down (predation) controlled is one of the most fundamental research questions in ecology (Lynam et al., 2017). Marine ecosystems, originally thought to be mainly steered by bottom-up control, periodically exhibit top-down control, particularly when predators are introduced into the ecosystem (Grosholz et al., 2000) or, conversely, when predators are eaten by fish (Llope et al., 2011). An important issue is how the changes in quantitative characteristics of individual species relate to the functioning of benthic communities, their role in organic matter transformation and their long-term sustainability under unstable environmental conditions (Górska and Włodarska-Kowalczuk, 2017). Fjord ecosystems are rather fragile and have lower functional redundancy of macrofauna in comparison with the open sea, especially in the upper part with intensive run-off. Global climate change results in retreating glaciers discharging meltwater with substantial amounts of particulate material, leading to an increase in mineral particles sedimentation (Wêsławski et al., 2011). The elimination of functional groups sensitive to glacial sedimentation leads to the simplification of the functional diversity, accompanied by a decline in species redundancy, i.e., the number of species performing similar functions in a system. In turn, reduction of functional complexity can indicate the higher sensitivity of fjord ecosystems to different disturbances (Włodarska-Kowalczuk et al., 2012). This is particularly important regarding the process of the ongoing snow crab Chionoecetes opilio invasion, which expanded over the entire western Kara Sea shelf in less than 5 years after the initial records (Zalota et al., 2018) and now have nursery areas in the bays of the eastern coast of the Novaya Zemlya Archipelago (Zalota et al., 2019).