Decadal Stability of Macrobenthic Zonation along the Estuarine Gradient in the Ob Bay, the Largest Siberian Estuary

Abstract

The Siberian Arctic Shelf is an area of increasing anthropogenic activity and is predicted to be significantly impacted by climate change. Accordingly, there is an increased demand for baseline knowledge on the scales of natural spatiotemporal variation of different ecosystem components. The study aimed to investigate the spatial variability of macrobenthic communities and associated abiotic forcing in the Ob Bay, a major Arctic estuary. Four main zones of macrobenthic communities were identified: a brackish-water zone in the upper part of the estuary, which was divided into two subzones according to the dominant species; a transition zone at the mouth of the Ob Bay; and a marine zone. This zonation remained stable during the study period (2013–2019) and corresponded well with previous studies in the area. The large-scale variation in macrobenthic communities was related mainly to two independent drivers: salinity and sediment type. The within-zone variation increased with the number of coexisting species, but no temporal trends could be assessed. The study highlights the need to account for the small-scale heterogeneity of benthic communities to understand ecosystem functioning and long-term dynamics, particularly in areas where environmental conditions vary markedly.

1. Introduction

The structure and distribution of benthic fauna are determined by the spatial heterogeneity of the environment, which, in turn, is supported by a variety of both biotic and abiotic factors [1]. There are numerous studies on the relationships between spatial scale and community variability [2,3,4,5,6,7]. At large scales (tens or hundreds of kilometers), abiotic factors such as main currents, sediment composition, and salinity gradients drive benthic structure and distribution [1,8]. At small scales (meters and tens of meters), inter- and intraspecific interactions (competition, predation), life cycles of organisms, larval behavior, local substrate topography, and microscale hydrodynamics, and various interactions between these biotic and abiotic factors have a major influence [3]. The choice of an appropriate scale is decisive for research and the interpretation of results from single-scale studies is often insufficient [4,9,10]. The issue of spatial scale is widely covered in the recent literature; however, studies addressing such an important aspect as the relationship between temporal and spatial variability are considerably fewer [11,12,13,14,15,16,17].

Estuaries are considered to be systems with enormous variability in ecological processes at different spatiotemporal scales; they are transition zones between freshwater and marine systems with high levels of dynamism and productivity [18,19]. The mixing of sea and river water results in a strong gradient of oxygen, temperature, and salinity in the water column, while sedimentation of riverine particulate matter is a major factor determining sediment composition [20]. Salinity is considered to be a key driver of the distribution of estuarine organisms and, hence, of changes in the species composition of benthic fauna [21,22]. The salinity gradient acts as a physiological barrier for stenohaline marine and freshwater species and places environmental stress on euryhaline marine species [22].

Benthic fauna is the slowest-responding component that reflects the biological integrity of the ecosystem and plays an essential role in estuarine processes at the sediment-water interface, such as biogeochemical processes of organic matter and nutrient cycling [23,24]. Moreover, estuarine benthos provides trophic support for resident and migratory fishes and birds [25]. Considering the significance of macrofauna for benthic ecosystem functioning, accumulating knowledge on the scale and processes that shape the distribution of macrobenthic communities in estuarine systems and coastal areas influenced by freshwater runoff is of paramount importance. Evaluating the natural magnitude of spatial and temporal variability in benthic communities associated with seasonal dynamics and long-term climate trends helps to distinguish the long-term response of ecosystems to human activities, given that strong natural gradients in highly complex transitional systems make estimating anthropogenic stressors difficult [26,27,28,29,30]. This may be particularly true for the Arctic region, where changes in all climatic parameters and associated ecosystem responses have been documented [31], yet relatively little research has addressed the Arctic benthic ecosystems under the direct influence of freshwater runoff [32,33,34,35,36,37,38,39].

The Ob Bay is the largest estuarine area in the Arctic region. Until recently, it was considered a relatively “pristine area” with a low level of anthropogenic impact [40]. Based on hydrological parameters related to the dynamic characteristics of the water column, two main zones can be distinguished in this estuary: “freshwater” and “marine” [41]. The location of the frontal zone between them exhibits distinct variability due to large seasonal fluctuations in river discharge [42,43]. Conservative benthic communities show a more discrete and stable pattern of estuarine zonation and, consequently, reflect the long-term dynamics of abiotic parameters. The most intensive studies in the southern regions of the Kara Sea, including the Ob and Yenisei estuaries, were conducted in the 1990s [32,33,44,45]. According to the latest research in the Ob Bay [36], three benthic communities—marine, brackish water, and freshwater—successively replace each other along the salinity gradient. However, these investigations did not consider the long-term variability of benthic fauna. In the present study, we combine data from three surveys to provide a comprehensive description of macrobenthic communities in Ob Bay. In addition, we aimed to test the hypothesis that the spatial boundaries of benthic communities are stable on the decadal scale by comparing the original data of the present study (2013–2019) with the previous investigation in the study area.

2. Materials and Methods

2.1. Study Area

The Kara Sea is characterized by the strong influence of continental runoff [46]. Most of this volume comes from the runoff of the Yenisey (630 km³ yr⁻¹) and Ob (400–500 km³ yr⁻¹) rivers [46,47]. The Ob Bay extends from approximately 66° N to 73° N. Most of the bay is 10–15 m deep, reaching a maximum depth of about 28–30 m, and its width varies from 35 to 80 km [47]. It is a micro-tidal, highly stratified estuary [48,49]. The hydrological regime of the Ob shows large seasonal and annual fluctuations—low discharge is observed during the cold season (October–May), while flooding is common in the summer and early autumn [46]. In terms of hydrophysical and hydrochemical parameters, two zones are clearly distinguished: “riverine”, where the key driving factor is fresh river flow, and “marine” [41]. The position of the boundaries of these areas can vary significantly, according to [50], up to 200 km, both seasonally and from year to year. In addition, in the extended eastern part of Ob Bay, there is a pronounced area of nutrient accumulation where chemical and biological processes differ from those in the area of the main flow [51,52]. The river flow slows down in this part of the bay, and seawater inflow is possible [49]. Sediments are commonly represented by homogenous silts and clays replaced by silty sands in marine parts [53] and consist mainly of terrestrial material transported as suspended matter. Sedimentation rates within the boundaries of the Ob-Yenisey shallow water can be as high as 1 cm per year [54]. The concentrations of suspended particulate matter are much higher in the upper part of the Ob Estuary in comparison with open water outside of the bay [55].

2.2. Sampling and Sample Processing

Sampling was conducted in three expeditions: the 125th cruises of RV “Professor Shtokman” in September 2013, and the 66th and 76th cruises of RV “Akademik Mstislav Keldysh” in July 2016 and July 2019, respectively. In total, 35 stations were collected in the Ob Bay and on the adjacent shelf (Figure 1) (Supplementary Table S1). ‘Ocean’ grab or Van Veen grab (0.1 m²) were used for sampling (3–5 replicates per station) [56]. 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. The resulting species list was carefully checked for validity and unified according to the World Registry of Marine Species database (WoRMS, https://www.marinespecies.org). Molluscs were weighed with shells, and polychaetes were weighed without tubes. Complete station data are given in Supplementary Table S1.

Concurrently with macrobenthos sampling, hydrological parameters (salinity, temperature, and turbidity profiles of the water column) were measured at all stations with SBE911 CTD (Supplementary Table S1). Dissolved oxygen (mL/L, %) and nutrients (phosphate, dissolved inorganic silicon, nitrogen in ammonium, nitrite, and nitrate forms) were determined immediately after sampling in the onboard laboratory according to the accepted methods [57]. These parameters can be used both to trace freshwater input and as indicators of organic matter mineralization [58]. The granulometric composition of sediments was analyzed by an Analysette 22 MicroTec Plus particle-size laser counter. The volume ratio of the mud fraction (silt and clay, <50 μm) in the sample was further used in analyses.

2.3. Data Analysis

All environmental variables were normalized and subjected to principal components analysis (PCA) for ordination. Permutational analysis of variance (PERMANOVA) was used to confirm spatiotemporal variations in the environmental variables [59].

For multivariate analyses of biotic data, we used the species respiration rate R as a measure of abundance, estimated as follows:

$$\mathrm{R}={\mathrm{k}}_i{\mathrm{N}}_i^{0.25}{\mathrm{B}}_i^{0.75}$$

where N_i is the abundance of a species, B_i is the biomass, and k_i is a taxon-specific coefficient. Then, the relative respiration (metabolic) rate for a species i was calculated as a proportion of the total: ${\mathrm{r}}_i={\mathrm{R}}_i/(\sum {\mathrm{R}}_i)$ [17,60,61]. Relative respiration rate describes the relative contribution of species to the production and metabolism of a community and is expressed as the percentage of each taxon. Combining the abundance and biomass of each species, this indicator gives a balanced contribution of small but abundant species and large species with low abundance but high biomass [17,36,61]. Prior to multivariate analyses, data were square root transformed to reduce the role of dominant species [62]. Finally, the similarity between stations was estimated using the Bray–Curtis similarity index [63] based on species-transformed relative respiration rates data. Non-metric multi-dimensional scaling (nMDS) was applied to identify key patterns in community structure [64]. The correlation between environmental variables and the MDS axes was assessed (Pearson correlation > 0.6). The ANOSIM procedure was performed to test for statistically significant differences among groups of samples in terms of species structure [64]. The resulting groups of stations were analyzed using the similarity percentage routine (SIMPER) to determine which species were responsible for the intra-group similarities [64]. In addition, the number of taxa, total abundance, biomass, and univariate diversity measures (Shannon index (H’loge), the expected number of species per 200 individuals (ES200), Pielou’s evenness (J)) were calculated for each group of stations [64]. A shade plot was constructed to visualize the differences in species structure between the defined groups (respiration rate data) [62]. According to the identified groups of stations on the nMDS ordination, only the key species that provided a high proportion of within-group similarity (by SIMPER) were selected for the analysis. Taxa were grouped in clusters using the UPGMA algorithm based on the index of association, and coherent species sets were identified by the Type 3 SIMPROF tests.

The large-scale relationship between the species data set and measured environmental variables as well as the spatial distribution of macrobenthic communities was investigated using variation partitioning in R package Vegan ver. 2.6–4 [65,66]. Two groups of variables were used: “Spatial” (latitude and longitude) and “Environmental” (measured bottom water parameters, mud content, and depth). Several stations from 2013, for which a complete environmental dataset was not available, were excluded from the analysis. Prior to variation partitioning, three groups of factor variables were separately subjected to forward selection based on adjusted R² with the significance level as stopping criteria. We used dbRDA [67] based on Bray–Curtis dissimilarities derived from standardized square root transformed relative respiration rates both for forward selection and variation partitioning. Subsequently, separate variation partitioning was performed for environmental variables (measured bottom water parameters, mud content, and depth) to assess the proportion of unique variation associated with each group. The small-scale spatial heterogeneity was assessed by calculating multivariate dispersion from Bray–Curtis dissimilarities (PERMDISP, average distance to centroid) within stations based on individual grab samples.

Statistical analyses were performed in PRIMER V7 with the PERMANOVA+ add-on package [62] and in R (version 4.2.2R, Development Core Team, 2022).

3. Results

3.1. Environmental Setting

The northern part of Ob Bay is occupied by a highly dynamic estuarine water mass characterized by a two-layer water structure, with a warm freshwater layer on top and a cooler, saline layer near the bottom. Only in the southern part of the investigated area (St. 16, 5324), the distribution of the main hydrological characteristics is completely different, and we found a quasi-homogeneous structure throughout the entire water column. Salinity values were 0.17–0.28, and temperature changed gradually from 12–13 °C in the surface layer to 8–9 °C near the bottom. To the north, stratification of water masses was observed at most stations in all years studied, but the thickness of the near-bottom layer varied throughout the estuary and from year to year. In the brackish zone to the south of 72°30 N, a thick desalinated surface layer was observed with salinity values of 0.1–1.7. The near-bottom, more saline waters occupied only a thin 1–5 m layer. Salinity varied here within a large range from 13 to 29, and temperature from 5.3 to −1.1 °C. This area of the bay is extremely dynamic, and it is possible to detect significant differences in the characteristics of the near-bottom water in different years (even the complete absence of a saline layer) (Supplementary Table S1). These dynamics are most likely determined by the inflow of saline seawater through a 10–12 m-deep sill located at 72°30 N [44]. To the north, towards the sea, the depth increased, and the stratification and thickness of the near-bottom layer became stable. While surface salinity varied from 0.3 to 11 and became 24.2 only on the shelf (St. 5309, 6240), the near-bottom layer was rather stable: salinity was marine (28.5–32.5), the temperature was near zero or negative (−0.6–−1.5 °C), and the near-bottom layer had a thickness of 10–20 m. Therefore, no inter-annual changes in the characteristics of near-bottom waters occurred in this area, and the hydrological parameters were relatively stable.

Water oxygen saturation varied greatly in the near-bottom layer in time and space (Supplementary Table S1), but oxygen deficiency was not recorded anywhere. High values of oxygen saturation (80–90%) were observed in the near-bottom layer. The range of nitrite in the form of nitrogen concentrations varied from 0 to 1.2 μM. Nitrate concentrations varied markedly, with minimum concentrations of nitrate nitrogen occurring in the shelf area (0.2 to 2 μM), while concentrations near the freshwater ranged from 0.8 to 8 μM. Phosphorus concentrations varied from 0 to 3.6 μM (Supplementary Table S1). The lowest values for phosphate and silicate were in the marine part of the bay, which may be related to strong stratification and hence insufficient water exchange between the layers. Silica concentrations in the fluvial part ranged from 20 to 230 μM since silica concentrations in estuarine waters are higher than in marine waters.

Sediments in most of the Ob Bay generally consisted of various silts and clays with an average particle diameter of 2–20 μm (Supplementary Table S1). A significant fraction of sand was observed in the adjacent shelf area (St. 1, 5313, 5310, 6241).

PCA based on the environmental data showed that the first three principal components accounted for 68% of the total variability (Figure 2). The first axis (41.4% of variations) divided the stations according to the natural gradient from estuarine to marine environments; the variability was mainly explained by the salinity gradient, changes in temperature, depth, turbidity, and dissolved silicon. Variability along the second axis (15.4% of total variations) was mainly explained by changes in dissolved oxygen, phosphate, and NH₄. The third axis (11.2%) was related to the mud content in sediments. Inter-annual differences in environmental variables as well as spatial variations along the estuarine gradient from south to north were revealed. Significant inter-annual differences were found for PO₄, NH₄, and O₂, and no changes were observed for salinity, temperature, Si, NO₂, or NO₃ (Permanova routine, Supplementary Table S2).

3.2. Macrobenthic Communities

A total of 186 macrobenthic taxa were identified during the three surveys. The most diverse groups were polychaetes (61 species), followed by crustaceans (54 species) and bivalves (26 species). Multidimensional scaling of the entire species respiration rates dataset revealed four distinct groups of stations at the 30% level of Bray-Curtis similarity (ANOSIM test, R = 0.956; p = 0.01%) (Figure 3). The identified macrobenthic communities reflected changes in salinity, temperature, and mud content, and were named according to the spatial distribution of zones within the estuary: Brackish—“B1” and “B2”, Transitional—“T”, and Marine—“M”. Differences between zones M, T, and B1 and B2 were mainly driven by changes in salinity, while sediments played an important role in distinguishing zones M and T as well as B1 and B2. Brackish-water zones B1 and B2 were located in the upper part of the estuary. Near-bottom salinity ranged from 0.12 to 25 in zone B1 and from 3.9 to 29 in zone B2 (Supplementary Table S1). In zone B1, the bottom temperature was always above 0 °C and reached 9.64, while in zone B2, it varied from 1.78 to −1.07 °C. In the transition zone (T), salinity varied from 28 to 33, and temperature ranged from −0.55 to −1.5 °C. Further offshore, a marine zone (M) extends across the shelf; in this zone, salinity and bottom temperature are stable at around 32–33 °C and −1.5 °C, respectively. At stations in zone B2, the sediments were almost 100% composed of mud, while the sediments in the marine zone were mostly sands with a mud content of 17–39%. In the transitional zone (T), the proportion of mud was significant (80–99%); in zone B1, the sediments were heterogeneous and consisted of sand or mud at different stations (Supplementary Table S1).

On a large scale, the macrobenthic community structure was spatially constrained, with the majority of variation explained by the spatial factors (21%) being related to the environmental drivers (mainly salinity, mud content, and depth), reflecting the estuarine gradient. (

Table 1. Results of forward selection for two groups of variables to be included in variation partitioning. For each variable cumulative adjusted R² and p-value are shown. Significant (p-value < 0.05) variables are in bold.

Spatial			Environmental
	Adj. R2	p-Value		Adj. R2	p-Value
Latitude	0.188	0.001	Salinity	0.165	0.001
Longitude	0.251	0.004	Mud	0.270	0.001
			Depth	0.326	0.002
			T	0.351	0.029
			Turbidity	0.384	0.006
			O2 (mL/L)	0.400	0.083

; Figure 4A). However, part of the variation (18%) explained by the environmental variables was not attributed to the spatial variables, indicating more complex environmental forcing on the macrobenthic communities. All three groups of environmental factors (mud content, bottom water, and depth) explained an independent fraction of the variation in the dataset. Most notably, mud content in the sediments appeared to affect the macrobenthic community structure largely independently from depth and bottom water. The small-scale spatial heterogeneity (quantified as within-station multivariate dispersion) was lower in the brackish water part of the estuary (Figure 4B). Figure 5A shows the division of the Ob Estuary into four spatial zones, successively replacing each other along the estuarine gradient from south to north and differing both in community structure and basic quantitative characteristics (Figure 5B and Figure 6;

Table 2. Mean integral characteristics of the identified zones. Dominant species included only species with a respiration rate greater than 5%.

Zone	Year	Number of Stations	N, ind/m2	B, g/m2	S	J’	H’ (loge)	ES(200)	Dominant Species (by SIMPER)	Contrib (by SIMPER)
B1 Avg sim 46.77	2013	3	3879 ± 901	19.9 ± 8.4	7	0.50	0.89	4.96	Oligochaeta Marenzelleria arctia Ampharete vega Pontoporeia femorata	30.54 22.62 18.78 11.85
	2016	3	2563 ± 1321	11.2 ± 7.1	10	0.47	1.07	6.78
	2019	2	4148 ± 1280	8.9 ± 5.1	7	0.46	0.84	5.40
	avg	8	3453 ± 1232	13.9 ± 8	8	0.48	0.94	5.75
B2 Avg sim 50.54	2016	3	1206 ± 503	18.9 ± 9.5	13	0.68	1.75	11.19	Portlandia aestuariorum Marenzelleria arctia	55.30 16.02
	2019	2	503 ± 212	68.6 ± 74.2	12	0.63	1.53	10.35
	avg	5	924 ± 411	38.8 ± 46.5	12	0.66	1.66	10.86
T Avg sim 42.40	2013	8	1789 ± 1050	90.6 ± 79.2	43	0.62	2.28	25.72	Portlandia arctica Cistenides hyperborea Saduria sabini	47.27 10.85 9.44
	2016	2	897 ± 391	37.9 ± 18.3	21	0.72	2.17	17.06
	2019	6	914 ± 367	85.6 ± 66.1	26	0.68	2.13	20.78
	avg	16	1350 ± 881	82.1 ± 68.6	34	0.65	2.21	22.79
M Avg sim 43.04	2013	1	3456	478.0	73	0.73	3.13	37.03	Macoma calcarea Maldane sarsi Astarte borealis	20.63 8.12 6.40
	2016	3	2588 ± 1059	165.8 ± 142.6	48	0.71	2.76	31.43
	2019	2	2098 ± 35	111.3 ± 72.8	62	0.74	3.04	37.46
	avg	6	2569 ± 834	199.7 ± 168.8	57	0.73	2.92	34.37

).

The lowest biomass and diversity values were in the brackish zone B1, while abundance in this zone was the highest due to the extremely high density of small oligochaetes and the polychaete Marenzelleria arctia (up to 4000 ind/m² and 2460 ind/m², respectively). The overall abundance was the lowest in the brackish zone B2, which also had relatively low species diversity and biomass (Table 2, Figure 5B). The main difference between the zone B2 and the neighboring zone B1 is the dominance of the bivalve Portlandia aestuariorum in the former, which was only found at the stations in this area, accounting for up to 55% of the intra-station similarity contribution by SIMPER (Table 2, Figure 6). Its abundance varied from 63 to 370 ind/m² and biomass ranged from 10 to 106 g/m². The station 5326 occupied an intermediate position between zones B1 and B2 (Figure 6). Although P. aestuariorum was also detected at this station, its biomass was only 2.8 g/m² for a total of 57 ind/m². The similarity between this station and the neighboring B1 group was due to the relatively high abundance and biomass of species such as the polychaete M. arctia (470 ind/m²; 3.7 g/m²) and the priapulid worm Halicryptus spinulosus (27 ind/m²; 4.6 g/m²). In addition, unlike all other B2 stations, where the large isopod Saduria sabini was observed, at the station 5326, Saduria entomon was found, which is common in the B1 zone.

The transitional zone T is located at the mouth of the Ob Bay (Figure 5A). Species diversity in this area was higher than at the brackish-water stations but significantly lower than in the marine zone (Table 2; Figure 5B). The same can be stated for total benthic abundance and biomass values. The bivalve Portlandia arctica was the dominant taxon, accounting for almost half of the contribution to intra-group similarity by SIMPER (47%). Its abundance varied from 170 to 1090 ind/m² and biomass ranged from 2 to 111 g/m². The polychaete Cistenides hyperborea and the isopod Saduria sabini were also among the dominant species, contributing 11% and 9%, respectively, to the intra-group similarity by SIMPER (Table 2, Figure 6).

The greatest diversity was observed in the marine area M, with total benthic abundance and biomass also being the highest here (Table 2; Figure 5B). The species composition of this community was very heterogeneous, but the bivalve Macoma calcarea and the polychaete Maldane sarsi were dominant at all stations (Table 2; Figure 6). M. calcarea abundance varied from 23 to 90 ind/m² and biomass from 2 to 72 g/m², while M. sarsi abundance varied from 23 to 90 ind/m² and biomass from 2 to 72 g/m². In addition, at the station 1, a high abundance of the bivalve Astarte borealis was recorded in 2013 (290 ind/m²) with a biomass of 341 g/m², while in 2016 (station 5313), it was almost halved (113 ind/m² and 159 g/m²). Unfortunately, we have no data available for this station for 2019, but at the northernmost station 6240, only three specimens of this species were found.

4. Discussion

4.1. Environmental Filtering in Shaping Macrobenthic Communities

In areas with strong environmental gradients, environmental filtering rather than biotic interactions, population stochasticity, and dispersal determine community structure and species coexistence [68,69,70]. The presence of a strong gradient in the distribution of communities, both benthic and pelagic, depending on environmental conditions is a key characteristic of any estuarine ecosystem and has been widely described in the literature [12,19,22,71,72,73]. Our results corroborate the existence of an estuarine-to-marine gradient and related zonation of macrobenthic communities in the Ob Bay described in previous studies [32,33,36,44]. We could distinguish a brackish-water zone (B1 and B2) in the upper part of the estuary, a transitional zone (T) at the mouth of the Ob Bay, and a marine zone (M) that extends across the shelf towards the sea (Figure 5). The observed macrofauna pattern wasdriven by freshwater inflow and sediment composition, which are known to be the most important abiotic variables for estuarine benthic species [74]. Despite the importance of chemical processes in the water column and thus food availability as a determinant of benthic communities (i.e., the benthopelagic coupling), hydrochemical variables (nutrients) did not explain any patterns in benthic communities, the response of which to nutrient changes in the water column is likely to be indirect and delayed [24]. Salinity has been the most important driver of changes in the species composition and structure of benthic communities throughout the estuary. High variability in salinity represents a physiological barrier for many species adapted to particular salinity range, determined by their physiological osmoregulation mechanisms, and limits their spatial distribution across the estuary [22,75,76]. Zones B1 and B2 undergo similar changes in salinity, with periodic influxes of waters with “critical salinity” (salinity 5–8) [75,76], which leads to the dominance of physiologically adapted oligochaetes and polychaetes, M. arctia, as previously reported [32,33,36,44,45]. Oligochaetes are the characteristic species in estuaries due to their recognized tolerance to low salinity and affinity for high mud and organic content in the sediment [39,71,77,78]. One of the differences between zones B1 and B2 at the species level is the presence of the isopod Saduria sabini only in the zone B2, probably due to the hydrodynamic regime and the influx of saline water. The eastern part of the Ob Bay, occupied by the B2 community, is distinguished from the main river flow by bottom topography and hydrodynamics, primarily because of the presence of the so-called Obsky Sill at 72°15′ N–72°25′ N with a depth of about 10 m [42]. Along this sill, there are troughs through which saline water enters this zone. All three species of Saduria found in our study tolerate a wide range of salinities and are euryhaline osmoregulators, although with varying degrees of effectiveness. The most marine species, S. sabini, has the lowest tolerance, whereas S. sibirica is intermediate, and S. entomon favors brackish water [79].

We also observed a difference in the distribution of the two bivalve species of the genus Portlandia: P. aestuariorum and P. arctica; the first one was observed only in the B2 zone, while the second inhabited the transitional zone T. Portlandia species are widespread in shallow inshore areas with muddy sediments both in the Kara Sea [80,81,82] and in the other Arctic seas [37,83,84,85,86,87], as well as in the areas associated with glacial sedimentation [88,89,90]. Portlandia aestuariorum indicates the areas of the lower part of estuaries exposed briefly to strong freshening [91] and can survive an exposure tosalinities less than 1 [92]. Portlandia arctica, on the other hand, generally requires a salinity of more than 25 [93], suggesting stratification between low-saline surface water and denser bottom water. The distribution of these species indirectly confirms the presence of the northernmost seasonal position of the critical salinity zone between brackish-water zones B1, B2, and the transitional zone T, which is characterized by the constant presence of a frontal zone with high rates of sedimentation of terrestrial particulate matter.

In both the transition (T) and the marine (M) zones, there was a small range in the variability of the near-bottom salinity (28–33), and the main feature of the M zone was a constant negative bottom water temperature and a low percentage of fine fraction in sediments as compared to the zone T. Similar marine communities in this part of the Kara Sea were described as early as the middle of the twentieth century [80]. Stenogaline species, which do not enter the transition zone due to their osmoregulatory capabilities, inhabit only the marine area.

4.2. Problem of Scale

The concept of scale and the importance of heterogeneity in space and time are now widely recognized as key issues in ecology [1,5,7,9,15,94]. In our study area, the spatial distribution of macrobenthic communities corresponded well between different years, which might indicate that on a regional scale (the entire estuary), the zonation of macrobenthic communities within the estuary remained stable over time and is consistent with data obtained for the Ob Bay in 2007–2010 [36]. Moreover, communities with the same dominants were observed earlier in 1991–1996. Benthic communities dominated by Marenzelleria arctia, Pontoporeia affinis, Ampharete vega, Saduria entomon, and oligochaetes were found in the brackish water zone (referred to in this study as the B1 zone) [33], communities with the dominance of Portlandia arctica were described in the area consistent with our transition (T) zone [33,44], and communities dominated by Astarte borealis, Macoma calcarean, and Maldane sarsi were found in the shelf area that corresponds to our marine zone (M) [33]. This long-term stability may be due to the strong vertical stratification of the water layer, which provides a constant two-layer structure with more stable salinity and temperature conditions in the near-bottom layer.

Similarly, stable-over-time zonation of macrobenthic communities was observed in the shelf areas with freshwater runoff (e.g., the shelf of the Laptev Sea and East Siberian Sea ([37,38] and reference therein) and the areas with the low level of river runoff (e.g., Baydaratskaya Bay in the Kara Sea before anthropogenic impact [95]) and might be an overall characteristic for the Eurasian Arctic [96]. Therefore, it might be assumed that macrobenthic shelf communities are well-adapted to natural environmental fluctuations within their habitats, which can be approximated by the combination of bottom water masses and sediment type. Nonetheless, such a stable spatial pattern might be disturbed by either anthropogenic impacts or invasive species [97,98,99].

The observed small-scale variability, which increases with the number of species in the local community, is most likely explained by spatial heterogeneity, stochastic population dynamics, and small-scale resource partitioning by dominant species. For instance, the cyclic changes in population dynamics associated with recruitment success were observed for the common arctic marine bivalve Astarte borealis [100], and similar processes could affect within-zone variability in the marine part of the studied area. Among other factors, biotic interactions and density-dependent regulation are suggested as possible mechanisms of small-scale variability [101,102]. This highlights the need for further research into other factors influencing benthic fauna and for sufficient sample replication to have enough power to detect possible temporal changes. With sufficient replication, hierarchical analysis of spatial and temporal variability can provide information on the range and relative importance of processes most crucial for organisms (environmental factors, population dynamics, and biotic interactions) [12,15,103,104]. Our study shows that in order to assess, monitor, and predict the state of the environment, it is necessary to determine the appropriate scale of natural variability and distinguish between its general and site-specific components. Regional-scale studies alone will not be sufficient to monitor benthic communities and smaller spatial scales are required to make appropriate conclusions about spatial or temporal variation in benthic communities and their relation to environmental gradients [15,94,105]. Establishing monitoring stations with an increased number of samples is recommended to exclude the possible influence of small-scale variability on the long-term trends.