Macrozoobenthos of the Pechora Bay in 2020–2021 indicates a likely change of common bivalve molluscs in the Arctic estuary

The Pechora Bay is a hydrologically and ecologically important area of the Barents Sea but there are still gaps in our knowledge of biodiversity of the area, including macrozoobenthos. In the first half of the twentieth century, the Pechora Bay was noted as a type locality for the bivalve mollusc Portlandia aestuariorum. Only a few surveys of macrozoobenthos have been conducted since and the last work from this area indicated the absence of P. aestuariorum. In this paper, we described macrozoobenthos and hydrological conditions of the bay based on the data collected in field campaigns in 2020–2021. All estuarine stations corresponded to a monodominant community of Macoma balthica poor in biomass (27.29 ± 20.82 g m−2) and species richness (33 species of macrozoobenthos recorded from 10 stations). The seaward most station was occupied by a marine assemblage of macrozoobenthos dominated by polychaetes Nephtys longosetosa and Cistenides hyperborea. Macrozoobenthos of the bay forms an ecocline from estuarine to marine species along the gradient of salinity. No significant differences in macrobenthic abundance, biomass and species richness were found between the 2 years of sampling and described fauna corresponds well to communities sampled in 1990s. Portlandia aestuariorum was absent in all our samples, which supports the hypothesis of disappearance of this previously common bivalve mollusc from the area, but the question of the driver of the change in macrozoobenthos remains open and requires further investigation.


Introduction
The Pechora Sea in the south-east basin of the Barents Sea has been attracting noticeable attention in the recent year due to its ecological and economical significance-the Pechora Sea has been identified as an Ecologically and Biologically Significant Marine Area by the Convention on Biological Diversity (CBD 2022) and by the Arctic Council (AMAP/ CAFF/SDWG 2013) because of its importance as moulting and staging areas for marine birds and summer haul-outs of Atlantic walruses. At the same time, the Pechora Sea remains the location of the only offshore oil-producing project in the Russian Arctic, the Prirazlomnoye oil field.
The special issue on the ecology of the Pechora Sea was published in Polar Biology in 2019 comprising research articles looking at different components of the Pechora Sea ecosystems (Sukhotin et al. 2019). The special issue had a prominent focus on the environmental changes observed in the region and effects of climate change on marine ecosystems. The authors highlighted that overall, the state of knowledge on the Pechora Sea ecosystems remains fragmentary and providing baseline data on current state of the Pechora Sea ecosystems is of critical importance as it will allow us to determine changes and possible ecosystem shifts in the future (Sukhotin et al. 2019).
The Pechora Bay is of particular interest as it provides up to 85% of riverine input into the Barents Sea (Dobrovolsky and Zalogin 1982). The hydrological conditions of the bay are determined by the interaction of the freshwater runoff from the Pechora River and the saline Barents Sea waters entering northern part of the bay through the Gulyayevskiye Koshki archipelago (Byshev et al. 2001(Byshev et al. , 2003Nikiforov et.al. 2005). The maximum salinity in the Pechora Bay reaches 33-35 PSU during the winter baseflow when the riverine input is minimal. In spring, salinity declines as the continental runoff accumulates under the fast ice (Polonsky et al. 2007). The summer salinity depends on the river discharge volume and the spread of the low-saline waters and can be as low as 0-1 PSU in the surface layer and 10-20 PSU in the near-bottom layer. Temperature values during the freshnet conditions can vary from 0.8 to 12 °C in the bottom layer. The warmest waters of the Pechora Bay are located near the mouth of the Pechora River, cooling down as they move to the northeast . The southern part of the bay is more influenced by the Pechora River discharge, whereas the northernmost part is more influenced by the Barents Sea waters and is typically more saline (30)(31)(32)(33). Wind conditions can affect the water mixing, bringing discharge waters to the northern area of the bay and decreasing salinity to 10-15 PSU in the surface layer (Polonsky et al. 2007). In general, the thermohaline characteristics of the bay waters are subject to significant temporal and spatial variability (Adrov and Denisenko 1996).
The biological conditions of the bay, including macrozoobenthos (benthic invertebrates > 5 mm), have been studied sporadically. Some stations in the Pechora Bay have been sampled during the big benthic surveys throughout the twentieth century with the most detailed surveys in the area in 1990s (Dahle et al. 1998;Denisenko et al. 2003). In 2019, a study of macrozoobenthos of the Pechora Bay was published by Denisenko et al (2019) based on the samples collected from 22 stations in the central part of the bay in 1995. Also in 2019, a paper by Gebruk et al. (2019) described the shallow-water community of the bay based on the eight stations in the Russky Zavorot Peninsula collected in 2016.
In one of the earliest studies in the Pechora Bay in 1920, bivalve mollusc Portlandia aestuariorum was first described as a subspecies of Portlandia arctica (Yoldia arctica subsp. aestuariorum natio petschorae Mosevich 1928)-the subspecies was abundant in the area, and the bay should be claimed its type locality (Mosevich 1928). Later studies in the 1990s also indicated the presence of Portlandia as one of the subdominants in the macrobenthic community in the Pechora Bay (Denisenko et al. 2003. The intention of this study was to survey macrozoobenthos in the bay and collect several specimens of Portlandia for genetic studies to resolve the controversy around the status of the subspecies. During the first survey in 2020, no Portlandia were found at the position given by Mosevich, or other stations in vicinity (Dgebuadze et al. 2021). The lack of a common species raised a question of possible causes and consequences of such a quick change. To address this question, repeated sampling of macrozoobenthos in the area was conducted in this study together with more detailed analysis of hydrological information. The changes in freshwater supply were hypothesised as a possible driver of changes in estuarine zone. Interannual variability of the communities and populations in a dynamic environment was also predicted as another likely reason of observed changes. In the present study, the samples from 2020 (partly described in Dgebuadze et al. 2021) were combined with the newly collected samples in 2021 and hydrological data.
Thus, the aims of the present study were (1) to characterise hydrological conditions in the bay to provide context for biological processes; (2) to review the present state of macrozoobenthos in the Pechora Bay based on the field data collected in 2020-2021; (3) to compare macrozoobenthos of the bay with the earlier studies based on the samples collected in the 1990s.

Study area and sampling
The Pechora Bay is a large shallow-water estuary in the Pechora Sea, formed by the Pechora River. From the north it is limited by the Russky Zavorot Peninsula, and a group of small islands Gulyayevskiye Koshki (Fig. 1).
Samples of macrozoobenthos were taken for this study from aboard RV Kartesh in 2020-2021 with benthic grabs. In 2020, one sample was taken from each station with a Van-Veen grab with 0.2 m 2 capture area, and in 2021 a smaller grab with 0.1 m 2 capture area was used and four replicates were taken from each station. Ten different locations were sampled in total, seven in 2020, and seven in 2021, with four stations repeated in both years (Table 1).
Some of the sampling stations repeated the stations sampled by co-authors in 1995 (Denisenko et al. 2019), and station D0 repeated the sampling location reported in 1924 as type locality of P. aestuariorum (Mosevich 1928). Figure 1 illustrates the locations of the sampling stations for this study and historical sampling locations from Denisenko et al. (2019).
Environmental data recorded for each station with a CTD profiler YSI CastAway included water depth, temperature and salinity. Discharge data used in this study were obtained at the downstream-most gauge station at the Pechora River in Ust-Tsilma by Roshydromet specialists. Station is located about 300 km from the mouth of the Pechora, and not affected by surge effects (Polonsky 2012). Processed discharge data were prepared by Arctic Great Rivers Observatory (Shiklomanov et al. 2021).

Sample handling and data analysis
Bottom sediments from the benthic grabs were washed with pumped seawater over the 0.5-mm nylon mesh. In 2020, one sample per station was taken and preserved with 96% ethanol on board (for possible future genetic analysis). In 2021, three samples were preserved with 4% formalin solution and one with 96% ethanol at each station on board. Then all samples were washed and transferred into 70% ethanol solution for storage in the ashore laboratory.
Macrobenthic invertebrates were identified with the maximum level of certainty through optical microscopy (Olympus SZ2) using regional taxonomic keys (Zhirkov 2001;Naumov 2006;Buzhinskaja 2009;Vassilenko and Petryashov 2009). All taxonomic names were standardised using the World Register of Marine Species (WoRMS). All specimens have been counted and weighted (wet biomass) on Ohaus Adventurer scales with reported accuracy to 0.01 g. Bivalve molluscs and gastropods were weighed in shells. For all specimens of Macoma balthica (130 specimens in 2020 and 370 in 2021), morphometric measurements were taken with accuracy to 0.01 mm using Vernier calliper or microscope eyepiece with ruler and included the following: SL-Shell Length; SW-Shell Width, SH-Shell  Height (Depth) for further analysis on population structure. Here, all the specimens were briefly divided in two groups: less 5 mm SL ("juveniles") and > 5 mm ("adults"). Mean values ± standard deviations were measured for biomass (g m −2 ) and abundance (ind. m −2 ) to characterise macrozoobenthos. In addition, "relative production" parameter was calculated as an relationship between biomass and abundance based on the average exponent of annual production on body size for macrobenthic invertebrates following Clarke and Warwick (2001), as a simplified variation of initial equation (Kucheruk and Savilova 1985;Azovsky et al. 2000): where B refers to biomass, A refers to abundance, 0.73 refers to the average exponent of annual production on body size for macrobenthic invertebrates.
The "relative production" index is introduced here as an approximated relationship between biomass and abundance to address input both from abundant but low in biomass species and larger organisms that dominate biomass but occur in samples less frequently.
Statistical calculations were performed using free software PAST version 4.09 (Hammer and Harper 2008). Shannon diversity index was calculated to characterise the diversity of macrozoobenthos (Hammer and Harper 2008). To assess predicted species richness (S ) in the research area, a species accumulation curve was built with the Chao-2 type estimator expressed as the following: where H refers to samples, S obs refers to the total number of observed species and s 1 refers to the number of species found in exactly one sample.
The non-metric multidimensional scaling (nMDS) based on the Bray-Curtis similarity index and the hierarchical cluster analysis based on an unweighted pair group method with arithmetic mean (UPGMA) algorithm with Bray-Curtis similarity index were used to distinguish macrobenthic communities. SIMPER (similarity of percentages) analyses were then carried out to assess contributions of individual species to observed differences. A pairwise analysis of similarities (ANOSIM) analysis was performed, p values of each pair were given to ascertain statistical significance with α set to 0.05, and sequential Bonferroni corrections were applied.
Abundance/biomass comparison (ABC) plots are calculated in Primer v 7 to determine levels of disturbance of macrozoobenthic communities. Stress conditions of macrozoobenthos are further expressed through W parameter following Clarke and Warwick (2001)-W parameter presents species assemblages in r-K continuum and is often used as a measure of the environmental stress. W takes values in the range (− 1, 1), where + 1 corresponds to even abundance across species but biomass dominated by a single species (admittedly "normal conditions") and − 1 stands for the opposite and corresponds to the stress conditions, for more detailed formula explanation please refer to (Clarke and Warwick 2001). Constrained seriation algorithm (Hammer and Harper 2008) based on the presence-absence matrix for all years was used to order the species according to when they appeared or disappeared in the samples. To standardise observations only species that occurred in more than one sample were considered.
Canonical correspondence analysis (CCA) was used to assess input from environmental variables into the distribution of macrozoobenthos (Hammer and Harper 2008). Biomass data were used for site/species matrix; temperature, salinity and depth were used as environmental variables (temperature, depth and salinity).
Maps were generated using ArcMap v10.8.1. by the standard geoprocessing tools with the reference coordinate system UTM/WGS84 Zone 40N.

The Pechora River discharge and vertical salinity structure
The annual river discharge of the Pechora River was analysed for 1995 (corresponding to the period of macrozoobenthic Fig. 2 Distribution of Pechora River discharge at the gauge station located in Ust-Tsilma in summer period of 1995, 2020 and 2021. Red and blue stars show dates of sampling for 1995 and 2020-2021 years, respectively. River discharge data were provided by Arctic Great Rivers Observatory (Shiklomanov et al. 2021) sampling by Denisenko et al. (2019), 2020 and 2021 (corresponding to the period of sampling for macrozoobenthos in the present study). According to the multi-annual analysis by Magritsky et al. (2017) for 1935-2013, the mean average discharge volume of the Pechora River is estimated as 147-150 km 3 . During the spring thaw in May-June, peak discharge can reach as much as 5000-40,000 m 3 s −1 . Figure 2 shows distribution of the Pechora River discharge in the summer period of 1995, 2020 and 2021. Maximum discharge of 25,000 m 3 s −1 was observed in 2020, minimum discharge of 16,000 m s −1 (at peak) was observed in 2021. River discharge in 1995 can be characterised as lower than in 2020, but higher than 2021. Distribution of discharge in 1995 clearly demonstrates the prolonged spring thaw characterised by the second discharge peak, one month after the main peak. The total annual flow of the Pechora River in 2020 (151 km 3 ) is 30 km 3 more than in 2021 (120 km 3 ). Total river discharge values in 1995 (157 km 3 ) were similar to 2020.
To assess the spatial and temporal variability of water temperature and salinity in the Pechora Bay, in-situ measurements were taken from aboard RV Kartesh with CTD profilers at the same stations where the macrozoobenthic samples were collected. The measurements in 2020 and 2021 were carried out during the same time of year and within the same research area. Temperature and salinity values in the surface and bottom layers are shown in Table 1.
Temperature conditions in the Pechora Bay in 2020 vary from 10.9 °C in the near-bottom layer to 14.3 °C in the surface layer. Values in 2021 are quite similar-from 9.8 to 14.9 °C. A slight annual difference in temperature values is associated with the same dates of sampling. Figure 3 illustrates vertical salinity structure in the Pechora Bay in 2020 (a) (19-20 July) and 2021 (b) (19-20 July). In 2020, salinity values varied from 3 PSU for the subsurface waters at the river mouth (station D11), which is strongly affected by the riverine influence, to 24-26 PSU for the near-bottom waters of the seaward station D2, where saline Barents Sea waters are having a strong influence on the near-bottom water mass (Byshev et al. 2003;Rogozhin et al. 2023). The Pechora Bay waters were fully mixed from surface to bottom, evident from the vertical salinity profiles (Fig. 3) for the stations D5-D0. Seaward from the station D7, salinity rapidly increased from 18 to more than 25 PSU at the depth of 6-8 m.
In 2021, minimal salinity values were 18-22 PSU near the river mouth (station D15), and the maximum values were over 30 PSU in the seaward stations. The frontal zone, distinguishing the riverine plume and sea water, was near the station D5-relatively narrow fully mixed from surface to bottom zone of approximately 20 km with sharp salinity gradient-from 19 to 26 PSU.

Species richness, biomass and abundance
A total of 26,230 specimens corresponding to 33 species of benthic invertebrates were identified in the grab samples from the 10 sampling stations in the Pechora Sea in 2020-2021 at the depth range between 5.3 and 17.5 m. The most species rich groups were Polychaeta (14 species) and Malacostraca (11 species) followed by Bivalvia and Gastropoda with 2 species each, and singular species of Anthozoa, Clitellata, Nemertea and Priapulida. The mean species richness was 10 ± 3 species per station (from 4-6 to 12-15 species per station along the gradient of salinity). The full list of species with mean occurrences across the sampling stations is presented in Table 2.
The observed number of species (33) was close to that predicted by Chao-2 estimator (35.4 ± 12.2 species) (Fig. 4), therefore, the observed number of species in the samples represents approximately 95% of the total predicted number of species. The value of Shannon diversity index across the sampling area was relatively low (H′ = 1.61 ± 0.3 mean for all stations), corresponding to a community with a small number of taxa with a few individuals and a strong dominance of one species.
Macoma balthica dominated macrozoobenthic biomass at all stations both years, with the exception of the most seaward station DS sampled in 2021 (Fig. 6). Stations with higher biomass alternate with more barren stations with no clear spatial distribution patterns (Fig. 6).
In 2021, a peak of juvenile recruitment of M. balthica was observed with 262 of 370 specimens corresponding to the post-settlement juveniles with shell length under 5 mm. As is evident from the map in Fig. 7, the juvenile molluscs (< 5 mm shell size) were found across the whole research area, including more seaward stations (D5) and more estuarine stations (D11, D15, D16). Maximum number of juveniles (87) were found in the northeast of the bay at station D5.
To assess dynamics of species composition between the yeas, constrained seriation algorithm was used for the presence-absence matrix. Seriation allowed to sort species in a way that showed which species appeared in the research area and which species disappeared (Fig. 9). The occurrence and distribution of the subdominants within the M. balthica community demonstrates a gradual shift of macrozoobenthos from marine species (Ophelia limacina, Scoloplos armiger, N. longosetosa) to brackish (Saduria entomon, M. affinis, P. femorata) along the gradient of salinity without strictly defined assemblages (Fig. 9).
In 2020, the spatial profile of W along the estuarine gradient shows the overall decreasing trend from values corresponding to "stressed" (− 0.1) to "normal" (0.6) conditions along the bay with the most "stressed" station closer to the river mouth (Fig. 10). ABC curves are presented for the stations with minimum and maximum W values in 2020 and 2021 (stations D2 and D13 in 2020; stations D11 and D2 in 2021). Notably, in 2021 W values fluctuated without a clear trend-all but one values are close to mean (0.25) (Fig. 10).

Interannual variability
To assess differences in macrozoobenthos between the two years of sampling, SIMPER analysis was performed. Overall observed average dissimilarity values and contributions to observed dissimilarity are presented in Table 4. Statistical significance of the observed difference was assessed with ANOSIM analysis (Table 4)-notably the observed dissimilarity was not statistically significant for any of the three analysed parameters (biomass, abundance and "relative production") (p > 0.05). Therefore, no statistical differences in macrozoobenthos were observed between the 2 years.
CCA for pooled biomass data of two consequent years confirmed the existence of a gradient in species distribution and no prominent difference between the years (Fig. 11).   The analysis included most of the species (excluding singletons) and three environmental variables for the repeated stations. The first axis that coincides with salinity (S) and depth explains 88% of variation, the second axis explains 12%. The position of the stations on the gradient did not change from year to year. To strengthen the comparison with historic data collected in 1990s, stations sampled in this study were mapped onto the five macrozoobenthic assemblages distinguished by Denisenko et al (2019) in the Pechora Bay. Corresponding station locations were identified for Type 1 (O. limacina-community), Type 2 (M. balthica-community) and Type 3 (M. affinis-Oligochaeta Gen.sp. community). Abundance, biomass and occurrence were compared between the types and between the years (1995 and 2020-2021), full data are presented in "Appendix 2". ANOSIM did not identify statistically significant differences between 1995 and 2020-2021 (p = 0.9, M. balthica contributing to 42% of observed difference according to SIMPER). At the same time Bray-Curtis similarity index showed 21-36% similarity for abundance data between 1995 and 2020-2021 (Type 1 = 0.36, Type 2 = 0.25 and Type 3 = 0.21) and 16-35% for biomass data (Type 1 = 0.16, Type 2 = 0.35 and Type 3 = 0.27).  Table 4 Observed differences in biomass, abundance and "relative production" of macrozoobenthos of the Pechora Bay between 2020 and 2021 according to SIMPER analysis-overall dissimilarity and contribution to dissimilarity (SIMPER); significance of observed difference (ANOSIM) Biomass (wet) (g m −2 ) Abundance (ind. m −2 ) "Relative production" index Overall average dissimilarity (%) 54 70 57 Largest contribution (%) Macoma balthica, 77% Pontoporeia femorata, 23% Macoma balthica, 67% ANOSIM significance of observed difference p = 0.9; R = − 0.08 p = 0.4; R = − 0.01 p = 0.8; R = − 0.07 1 3

Hydrological conditions in the Pechora Bay
In general, long-term dynamics of the Pechora River discharge are characterised by two altering phases (Georgiadi and Groisman 2022): decreased flow (1936-1955; 1967-1980) and increased flow (1956-1966;1981-present) with difference of up to 500 m 3 s −1 (8-10 km 3 year −1 ) during freshet period, corresponding to less than 10% of total annual discharge. This is noticeably lower than interannual variability with the difference in maximum and minimum discharge in spring flood period between the consecutive years sometimes reaching 30-40 km 3 according to the analysis of multi-annual discharge data (Magritsky et al. 2017;Rogozhin et al. 2023). Overall, no significant changes in hydrological regime of the Pechora River have been reported in literature (Magritsky et al. 2018;Frolova et al. 2022). Estuarine mixing in 2020 was typical for shallow-water tidal plumes (e.g. the Khatanga plume), as opposed to water stratification in the estuaries with weak tidal forcing (e.g. the Yenisei plume) (Osadchiev et al. 2020). In the presence of a significant tidal component, distance and depth of distribution of freshened waters in the water area of the Pechora Bay increases (i.e. there is practically no stratification). This means that initial freshwater volume forms a deep, but diluted river plume (Osadchiev et al. 2020). The tidal currents in the research area are low speed (Nikiforov et al. 2005), and shown to have little impact in the southern part of the Pechora Bay (Polonsky et al. 2007). At the same time, the tide height in the Pechora Bay can reach 1-1.5 m, which can affect the mixing of waters . It is also possible that the observed water mixing in the Pechora Bay was due to the wind influence.
Interannual variability of external factors in the Pechora Sea region (river runoff, wind conditions and others) has a noticeable impact on the salinity values across the Pechora Bay (Rogozhin et al. 2023), which is why the salinity of both surface and bottom water in the research area in 2020 was lower than in 2021. In addition, in 2021 the water column was stratified across the plume area, and the frontal zone dividing Pechora Bay water from the more saline Barents Sea water was located 30-40 km closer to the river mouth than in 2020. Salinity values in 2021 were also different from 1995, according to Denisenko et al. (2019)-during the freshet period in 1995 the salinity values were significantly lower than 2021 and comparable with data from 2020, which is likely linked to the prolonged spring thaw in 1995. Thus, the frontal zone in 1995 and 2020 was outside of the Pechora Bay eastward from the Gulyayevskiye Koshki Archipelago. It is worth noting that the hydrological conditions of the Pechora Bay outside the survey period may be significantly different from those described in this section due to influence of tidal forcing, variability of river runoff and wind forcing. The review of historical hydrological data together with recent observations do not confirm the hypothesis of significant changes in hydrological regime as a likely cause of the disappearance of P. aestuariorum. The differences in freshwater supply are within the typical range of interannual variability known for the area from literature.

Macrozoobenthos of the Pechora Bay
In general, macrozoobenthos of the estuarine part of the Pechora Bay is formed by a small number of species with a constant predominance of bivalves M. balthica (by biomass). Low species diversity and relatively high abundance of euryhaline species is a common trait in estuarine ecosystems (McLusky et al. 1993;Elliott and Whitfield 2011;Denisenko et al. 2019). Strong dominance of one or a small number of species (2-4) is also a common in mesohaline environments (Vedenin et al. 2015;Gebruk et al. 2019). Despite the slightly different sampling locations in 2020 and 2021 (only four stations were the same in both study years), no statistical differences were found in the species composition between the 2 years.  (St D9, D11, D12 in 2020 andD11, D15, D16 in 2021). No strong borders occur between these variants and their delimitation is rather speculative as it is based on independent distribution of the restricted set of species. Assessment of the disturbance of macrozoobenthos showed lower W values corresponding to more disturbed conditions closer to the river mouth in 2020 (Fig. 10), which corresponds to the analysis conducted by Denisenko et al. (2019) and can be explained by the influence of the river outflow, including sedimentation of inorganic matter. However, in 2021 no clear trends were observed with all but one W values close to mean.
A stable transition from estuarine to marine community was noted in both years along the gradient of salinity. Salinity is widely accepted as a primarily environmental driver in many estuarine ecosystems (Elliott and Whitfield 2011;Vedenin et al. 2015). Notably, some researchers reported that in the Pechora Bay, the primary driver specifically for the macrozoobenthic species richness is the sediment structure . The observed in this study distribution of macrozoobenthos corresponds to an ecocline of species from marine benthic community in the seaward part of the bay to depleted set of brackish-water species in the inner part. There is a substantial body of scientific debate around the nature of ecological boundaries in estuarine ecosystems (reviewed in Attrill and Rundle 2002), but an ecocline model (i.e. a continuum of assemblages along the gradient of an environmental factor) has been described for many estuaries including the Yenisei Estuary and the Gulf of Ob in the adjacent Kara Sea (Vedenin et al. 2015;Udalov et al. 2021), the Thames Estuary (Attrill and Rundle 2002), the Forth Estuary (McLusky et al. 1993) and other aquatic ecosystems characterised by gradient of environmental conditions-e.g. the Spitsbergen Fjords (Włodarska-Kowalczuk et al. 2012).

Temporal variability of macrozoobenthos
Overall macrozoobenthos described in this study was similar to community described by Denisenko et al. (2019) from the central part of the Pechora Bay based on samples collected in 1995. In total, 22 species occurred in both studies, 10 species were unique to 2020-2021 dataset, and 77 species were only found in 1995 but not 2020-2021, which is likely due to larger number of stations sampled in 1995. Additional species occurred in 1995 were mostly marine species in the outer areas of the estuary closer to the open water. No evidence of changes in hydrologic conditions (i.e. position of frontal zone or freshwater supply) causing the reduction in stenohaline species have been found. More detailed study of the hydrological structure of the Pechora Bay over a long-term period is needed to gain a better understanding of the temporal variability of environmental conditions in the bay. The species list was also compared to species listed by Gebruk et al. (2019) from the near-shore shallows of the Pechora Bay based on samples collected in 2016 from the coasts of Russky Zavorot Peninsula. Five species occurred in 2016 but not in other years, and 8 species were found in all years (1995; 2016; 2020-2021). Among them Cyrtodaria kurriana, the third dominant bivalve in the 1995-survey, was rather abundant in the shallows, but did not appear in the deeper stations in 2020-2021. The species that were unique to 2020-2021 included the following-C. hyperborea, Prionospio cirrifera, Lamprops fuscatus, Priscillina armata, Cylichnoides occultus and Margarites helicinus. "Appendix 1" provides a full list of species found in 1995, 2016 and 2020-2021.
As demonstrated on Figs. 9 and 10, the nature of ecological boundaries in the Pechora Bay corresponds to an ecocline with gradual substitution of species and without clearly defined community boundaries. In general, the majority of macrobenthos that occupied central and northern parts can be interpreted as a variation of the M. balthica community, described by Denisenko et al. (2019), also corresponding to a reduced in richness and biomass form of M. balthica community described by Gebruk et al. (2019) at the periphery of its distribution in Russky Zavorot Peninsula. It also corresponds to faunal association typical for low salinity environments, described by Dahle et al. (1998) in the Pechora Bay based on one sampling station as a reduced in abundance and richness community of P. femorata, M. balthica and Halicryptus spinulosus. The lack of change between macrozoobenthos in 1990s and 2020-2021 ("Appendix 2") can be explained by the fact that high spatial heterogeneity of estuarine macrozoobenthos combined with noticeable fluctuations in consecutive years show stronger differences than temporal changes observed on a relatively short time-scale.
The low level of integrity of the described macrobenthic community is evident from the prominent differences in the abundance and population structure of the dominated species-M. balthica, in consequent years. The population structure of this species is likely driven by interannual dynamics and does not correlate with the abundance of other species and vice versa. Regular benthic surveys are needed in the future to assess the state of the M. balthica population and abundance. The profile of ABC curves also suggests low integrity of assemblages (Fig. 10)-rapid change from clear transition from "stressed" to "normal" conditions in 2020 followed by fluctuation of values in 2021 can be explained by low resolution of the method, and also by rapid interannual rearrangement of species due to low integrity of macrozoobenthic community (Clarke and Warwick 2001). The low integrity of benthic community suggests the lack of strong biological interactions between the species within the community-therefore it is unlikely that disappearance of Portlandia in the area between 1990 and 2020s is due to competitive exclusion. The reasons for the change in common bivalve molluscs in the area are not apparent from available hydrological and biological data.

Conclusions
Overall, the observed macrobenthos in the Pechora Bay represents an ecocline of species with high level of independence in population dynamics and strong affinity with the gradient of salinity. The observed macrozoobenthos corresponds in principle with the community reported in the bay in 1990s (Dahle et al. 1998;Denisenko et al. 2019)lower values were noted for species richness, biomass and abundance, but were not statistically significant. The lack of statistically significant differences between data from 1990s and 2020-2021 can be explained by a potentially high stability of estuarine macrozoobenthos, but also by the high spatial heterogeneity and interannual variability of macrozoobenthos in the area-to gain a better understanding of the temporal dynamics of macrozoobenthos in the context of ongoing climate change in the region, long-term timeseries data are essential. No significant changes were found in freshwater input during last decades (Magritsky et al. 2017;Rogozhin et al. 2023) to explain the lack of P. aestuariorum in the area, thus the nature of the disappearance of the previously common bivalve mollusc (Mosevich 1928) remains uncertain and requires further investigation.