Patterns of Benthic Communities in Arctic Fjords (Novaya Zemlya Archipelago, Kara Sea): Resilience vs. Fragility

Despite a large number of studies, a detailed overall picture of benthic communities zonation in the Arctic fjords is currently lacking. Our study aimed to find out whether there is a universal model for the distribution of benthic communities based on the structural features of the fjords. We examined benthic macrofaunal communities in fjords with various environmental settings on the eastern coast of Novaya Zemlya Archipelago, Kara Sea. The material was collected during five cruises undertaken from 2013 to 2016. A total of 50 stations located in the five fjords were taken. In all five fjords, macrofauna had a similar composition assembled from a regional species pool, with a predominance of species tolerant to glacial sedimentation and fluctuations in temperature and salinity. Benthic communities changed consistently along the axis of the bay from the outer slope to the inner parts. Biodiversity and quantitative characteristics of the macrofauna decreased along the environmental gradient related to terrigenous and glacial runoff, consistent with patterns reported in other studies of Arctic glacial fjords. The most impoverished communities were dominated by bivalve Portlandia arctica and isopod Saduria sabini. At the same time, fjord walls and sills, characterized by low sedimentation rates, strong currents and the presence of ice-rafted debris, were inhabited by patchy distributed benthic communities dominated by species confined to hard substrates. In general, the distribution of communities corresponded to five zones: depleted inner periglacial areas, the upper subtidal belt with stony substrates, deep inner semi-isolated basin, outer non-isolated basins and upper slope transitioning to lower slope. Our study can provide a reference point for monitoring changes in fjord ecosystems in response to climate change and the potential impact of human activities.


INTRODUCTION
Global climate change and increased human activity can significantly and unpredictably alter highlatitude ecosystems (Bluhm et al., 2011;Wassman et al., 2011;Polyakov et al., 2020). Therefore, scientific interest in studying the coastal zone of the Arctic seas continues to grow steadily. Polar ecosystems are unique in several aspects, such as low temperatures, pronounced seasonal changes, and variable levels of organic carbon input (Grebmeier and Barry, 1991), and the species that inhabit them are often adapted to extreme environmental conditions (Peck et al., 2004). The ecosystems of Arctic fjords and bays are currently among the most vulnerable due to their position at the land-ocean interface (Bianchi et al., 2020) that is influenced by terrestrial runoff through both desalination and transport of suspended organic and inorganic matter. They are highly dependent on several factors such as geomorphological features, degree of isolation from the open sea, climatic parameters, the volume of terrestrial runoff, coastline type, etc. The environment described above, in turn, may vary significantly from one fjord to another, resulting in an impressive diversity of benthic communities (Pearson, 1980;Svendsen et al., 2002). Climate change, including the widespread decrease of outlet glaciers (Vieli et al., 2002;Thomas et al., 2009) and general variability in seasonality and runoff volume, leads to transformation in the hydrological structure of the water column. This results in increased stratification and isolation of the inner parts of semi-isolated fjords, and, on the contrary, to increased interaction between weakly isolated bays and the sea, with consequent effects on community faunal patterns (Beuchel et al., 2006;Renaud et al., 2007;Wêsławski et al., 2011). It remains unclear to what extent the ecosystems of fjords and bays are resilient to such impacts? Are there similar spatial patterns across different bays that would allow us to predict the simultaneity and the magnitude of changes, or is every bay unique and will respond in its own way in a changing environment? The study of ecosystems in a variety of fjords (e.g., with different runoff patterns) within the same regional fauna can identify key factors of community formation and, subsequently, predict ecosystem transformation due to climate change and associated shifts in driving factors.
The information on the coastal ecosystems of most Arctic areas (especially the Siberian Arctic) is still scarce. The most intensively explored region in the Arctic is the coastal zone and fjords of Svalbard (Görlich et al., 1987;Renaud et al., 2007), with Kongsfjorden being the most studied fjord in this area (Włodarska-Kowalczuk et al., 1998Kendall et al., 2003;Włodarska-Kowalczuk and Pearson, 2004;Wiencke and Hop, 2016;Mazurkiewicz et al., 2019;Jima et al., 2021). There are also several detailed surveys focused on the Greenland fjords (Thorson, 1933(Thorson, , 1934Sejr et al., 2000), Northern Norway (Larsen, 1997;Holte and Gulliksen, 1998;Oug and Høiśter, 2000;Holte et al., 2005;Jorda Molina et al., 2019;Kokarev et al., 2021), Canadian Arctic (Dale et al., 1989;Syvitski et al., 1989;Aitken and Fournier, 1993;Aitken and Gilbert, 1995) and Franz Joseph Land (Włodarska et al., 1996). At the same time, the coastal zone, fjords and bays of Novaya Zemlya archipelago have been beyond the scope of arctic ecological research since the early 1920s, when the first attempts to study the coastal benthic fauna were made (Zenkevich, 1925;Mesyatsev, 1931;Zenkevitch, 1963). As a result of these surveys, the regional fauna and benthic assemblages of Chernaya Bay (Gurjanowa and Uschakow, 1928) and Matochkin Shar Strait (Uschakow, 1931) were described in detail. Soon afterward, the archipelago was closed to research, and scientific work in the region was suspended. For more than half a century, Novaya Zemlya was a closed territory used for nuclear tests and disposal of radioactive waste and non-operational atomic vessels (Aibulatov, 2000). The studies were relaunched only in the 1990s as part of an environmental monitoring project (Pogrebov et al., 1997;Galtsova et al., 2004), but there was almost no information until recently about the benthic communities of Novaya Zemlya bays. Although in several fjords nuclear reactor compartments, a submarine and containers of low level radioactive waste were found, there is no indication that nuclear fuel from the dumped reactors or submarine has been releasing detectable quantities of radioactivity into the marine environment (Dahle et al., 2009b;Gwynn et al., 2016;Heldal et al., 2018).
Another important fact is that, due to the limited access and severe environmental conditions impeding dispersal, survival, growth and reproduction of many species, the Arctic has been considered a region of low risk for biological invasions (Vermeij and Roopnarine, 2008;Ruiz and Hewitt, 2009). However, climate change and anthropogenic pressure contributed to the occurrence of invasive species (Ware et al., 2016;Ricciardi et al., 2017;Spiridonov and Zalota, 2017). For example, at the beginning of the twenty-first century, the snow crab, Chionoecetes opilio (Decapoda: Oregonidae), invaded vast areas of the Barents and Kara Seas with unprecedented speed for a shelf species (Zimina, 2014;Sokolov et al., 2016;Spiridonov and Zalota, 2017;Zalota et al., 2018;Bakanev and Pavlov, 2020). At the moment, crab reaches the highest density in the fjords of Novaya Zemlya eastern coast (Zalota et al., 2019). The expansion of this species may lead to significant changes in the entire structure of bottom communities and, therefore, the fjord ecosystem as a whole.
The present study is based on data collected in five bays on the eastern coast of the Novaya Zemlya archipelago (Udalov et al., 2016(Udalov et al., , 2018(Udalov et al., , 2019(Udalov et al., , 2020Chava et al., 2017) before the invasion of the snow crab. Comparison among the fjords with various environmental settings (geomorphology, depth, presence/absence of the sills, sea/bay water exchange, glacial or terrigenous runoff, sediments, etc.) having the same set of regional fauna species allows us to raise several issues concerning the peculiarities of bottom ecosystem formation. This study aims to determine whether there is a universal scheme of fjord benthic communities spatial pattern based on structural features within the bays. If so, we can associate specific communities to different parts of the fjords (such as sill, periglacial basin, outer slope) and construct a general scheme of benthic communities distribution in the fjord-type bays. An alternative hypothesis states that each fjord and its communities are unique and cannot be combined. Our main objectives were to (1) compare benthic communities in several fjords with varying habitat conditions; (2) identify the main driving factors which influence bottom communities structure; (3) check the uniformity of benthic communities under certain conditions. With this study, we hope to contribute to our knowledge about the distribution patterns of macrobenthos in Arctic fjords and provide baseline data for tracking future ecosystem changes.

Study Area
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 Frontiers in Ecology and Evolution | www.frontiersin.org 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 glacialderived (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(Udalov et al., , 2018(Udalov et al., , 2019(Udalov et al., , 2020Chava et al., 2017).

Sampling and Sample Processing
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 2 ) 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(Udalov et al., , 2018(Udalov et al., , 2019(Udalov et al., , 2020Chava 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).

Data Analysis
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 . 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) . 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 2 ), 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.

Environmental Data
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 2 = 0.854 and R 2 = 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 2 = 0.386 and R 2 = 0.485, respectively). Blagopoluchiya Bay is of the third type and is characterized by high turbidity values throughout the bay (R 2 = 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).

Benthic Communities
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). Shallowwater 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 2 ) was also 2-4 times lower ( Table 2).

Community Structure
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), IIIthe deep inner semi-isolated basin (cluster 6), IV-the outer non-isolated basins and upper slope (clusters 7-10), and Vthe 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) Dominant species included only species with the respiration rate greater than 10%.
FIGURE 4 | Shade plot of untransformed species respiration rates data. The key species that provided a high proportion of within-group similarity (by SIMPER) in station groups (according to the 11 clusters identified) were selected for analysis. Species grouped using UPGMA clustering based on index of association. The colors and symbols indicating clusters are the same as in Figure 3. Five fjord zones differing in the structure of benthic communities are marked with the following symbols: I, II, III, IV, and V. 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 nonisolated 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.).

Environmental Drivers
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 2 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).

General Pattern and Driving Factors
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øisaeter, 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øisaeter, 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).

Fjord-Specific Communities and Their Resilience to Environmental Stressors
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, 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, 1933Thorson, , 1934Dale 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 surfacedeposit-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 .
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.
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).

Future Perspectives
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).

CONCLUSION
Commonly observed patterns of environmental heterogeneity include continuous variation in the form of spatially structured environmental gradients and patchy variation in which the environment is a series of nested habitat patches (Scheiner et al., 2011). In the case of the fjords of the Novaya Zemlya Archipelago, both patterns were presented. In all fjords, there was a consistent change in benthic communities along the axis of the bay from the inner parts to the outer slope related to coastal and glacial runoff. However, it was more pronounced in bays with glacier outlets (Oga and Tsivolki Bay), less pronounced in bays with runoff from a glacier located in the inland part of the island (Blagopolopoluchiya Bay) and the weakest in deep bays with moderate river runoff (Sedova Bay). At the same time, fjord walls and sills, characterized by low sedimentation rates, strong currents and the presence of ice-rafted debris, were inhabited by patchy distributed benthic communities assembled from taxa associated with hard substrates, including suspension-feeding species. Hence, the patchy local landscape was overlapped with environmental gradient, resulting in various combinations of macrobenthic communities composed from a regional species pool.
However, at different spatial scales, local communities are not random subsets of regional species pools but assemble from species having similar tolerance to environmental stressors and shared functional traits (Somerfield et al., 2009). Our comparative study showed that glacial runoff and associated changes in environmental variables are the most important driving factors that define the patterns of distribution and changes in the species composition and quantitative characteristics of fjord benthic communities. Species tolerant to increased sedimentation and fluctuations in temperature and salinity have an advantage. Results presented in our study provide a framework evaluating the level of natural variability in the fjord benthic ecosystem, which could further be used as a reference point when monitoring possible changes.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
AU and VM contributed to conception and design of the study. AU, AC, AV, and SS collected and processed the data. AU and MC analyzed the data and performed the statistical analysis and wrote the first draft of the manuscript. VM supervised the writing of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

ACKNOWLEDGMENTS
The study was carried out within the framework of the program "Ecosystems of the Siberian Arctic Seas." We would like to thank M.V. Flint-the leader and organizer of the expeditions and the crews of R/V "Professor Shtokmann" and "Akademik Mstislav Keldish." We are particularly grateful to A. Basin for his help in taxonomic identification of amphipods, and M. Simakov and V. Semin for their assistance in identification of polychaetes.