Role of Macrofaunal Communities in the Vistula River Plume, the Baltic Sea—Bioturbation and Bioirrigation Potential

Simple Summary Coastal areas, especially river plumes, are very diverse and dynamic zones where numerous geological, chemical and biological processes take place. This is because fresh water from the river with all substances, including pollutants from the land, mixes with salt water from the sea, creating specific living conditions for the organisms that inhabit the area. These organisms, e.g., macroscopic invertebrates such as mussels or worms, live in the sediment where their movement and feeding activities cause the sediment to mix and allow water to flow through it—these activities are called bioturbation and bioirrigation. Our research aimed to investigate how the structure and functioning of benthic marine ecosystems change with distance from the river mouth. We found that coastal areas are very diverse and host a wide range of organisms that bioturbate and bioirrigate and support sediment transformations relatively deep (up to 15 cm) into the sediment. Farther away from the river mouth, organisms were very scarce and occurred only on the sediment surface and did not burrow into the sediment, so bioturbation and bioirrigation did not take place. The coastal zone is like a hotspot where ecosystem processes and services are intensively reflected, and this is especially important when deeper areas are not functioning properly, as in the Baltic Sea. For this reason, we should consider how we can support the protection and recovery of marine ecosystems. Abstract Macrozoobenthos plays a key role in the transformation of inputs from rivers to the sea, such as nutrients, organic matter, or pollutants, and influences biogeochemical processes in the sediments through bioturbation and bioirrigation activity. The purpose of our study was to determine the structure of benthic communities, their bioturbation (BPC) and bioirrigation potential (IPC), and the vertical distribution of macrofauna in the Gulf of Gdańsk. The study revealed changes in the structure of benthic communities and, consequently, in the bioturbation and bioirrigation potential in the study area. Despite the presence of diverse and rich communities in the coastal zone, BPC and IPC values, although high, were formed by a few species. Both indices were formed mainly by the clam Macoma balthica and polychaetes, although the proportion of polychaetes in IPC was higher than in BPC. In the deepest zones, the communities became poorer until they eventually disappeared, along with all macrofaunal functions. Both indices changed similarly with distance from the Vistula River mouth, and there was a very strong correlation between them. We also demonstrated that the highest diversity of the macrofauna was observed in the upper first cm of the sediment, but the highest biomass was observed in deeper layers—at a depth of up to 6 cm, and single individuals occurred even below 10 cm.


Introduction
Coastal zones provide a variety of benefits derived by humans from ecosystem functions and processes. These include nutrient regulation or waste treatment functions, where biota play an important role in storage, recycling or removal of nutrients and compounds [1]. All of these functions help maintain healthy and productive marine ecosystems. Coastal in both scientific research and environmental monitoring. According to Queirós and colleagues [38], the bioturbation potential index also has limitations, and knowing this can contribute to more informed use of the index as an indicator of benthic function.
A few studies on the role of macrofauna carried out in the Gulf of Gdańsk have addressed the bioturbation potential index (BP C ) or nutrient fluxes between water and sediments [39][40][41][42]. Studies on the functioning of marine ecosystems in the Gulf of Gdańsk in recent years have also considered the influence of organic matter on the structure and functioning of trophic networks [43] and how organic matter is transformed by organisms [44]. So far, no research has been carried out in the Gulf of Gdańsk on bioirrigation processes. There are also few published studies on the distribution of organisms in the sediment. They mostly contain information on the depth of occurrence of individual macrofauna and meiofauna taxa [45][46][47][48], but only single studies addressed entire benthic communities [40,41,49].
The objective of this study was to determine the structure of benthic fauna as well as the bioturbation and bioirrigation potential of macrofauna in the sediments of the Vistula plume area, in the Gulf of Gdańsk. Furthermore, it was determined quantitatively how this impact of benthic communities varies depending on the proximity of the Vistula River mouth, as well as which species are the most responsible for sediment matrix reworking in the area. In addition, we have made an attempt to investigate the vertical distribution of macrofauna taxa, detailing their maximal and typical depth of occurrence in the sediment.
The results presented in this paper will help to demonstrate the zones where, due to the presence of animals in the sediments and their activity, nutrients, organic matter and pollutants carried into the Gulf of Gdańsk by the Vistula River are processed. They will also provide knowledge of the vertical distribution of species in the sediments necessary, among other things, for indices of functionality to assess the functioning of the seafloor and basin. Determining the role of macrofauna will also provide arguments for the protection and proper management of marine areas in estuaries.

Sampling
Bottom water, sediment and fauna were collected during two cruises in the Vistula River plume area and along an offshore depth transect in the Gulf of Gdańsk, the Baltic Sea ( Figure 1). Samples from 11 sites were collected in July 2014 from the deck of RV Elisabeth Mann-Borgese. In March 2016, three more sites were sampled (VE04, VE06, VE07) during a cruise aboard RV Alkor. Bottom water temperature, salinity and dissolved oxygen (DO) concentration were measured at all sites approximately 0.5 m above the sediment using a Seabird CTD-system with an oxygen SBE43 sensor.
For sediment and macrofauna analysis sediment cores (inner diam. 10 cm) were collected with a multicorer and subsamples of coarse-grained sands were collected from a Haps corer. At each site the upper 10 cm sediment sample for sediment parameters was frozen and prior to all analysis the sediment was dried and homogenized. The organic matter content of the dry sediments was measured as the percentage loss on ignition (LOI) after dry combustion for 8 h at 450 • C and for 5 h at 550 • C for samples collected in March 2016. For grain size analysis, samples were sieved using a shaker and a set of standard test sieves with mesh diameters of 2, 1, 0.5, 0.25, 0.125 and 0.063 mm [50]. Based on a percentage of each class in the total sample mass, sediments were classified by the Udden-Wentworth grain-size scale (after Wentworth [51]).

Macrofauna
For benthic fauna analysis, 3 to 5 replicates were collected at each site, with the exception of station VE49, where only 2 replicates could be collected. Sediment cores were divided into layers: 0-1 cm, 1-3 cm, 3-6 cm, 6-10 cm, 10-15 cm and >15 cm depth. We sifted all layers separately through a 1 mm sieve to separate the macrofauna from the sediment and preserved with 4% formaldehyde until analysis (stored for at least 3 months).
Biology 2023, 12, 147 4 of 22 In the laboratory, the fauna was sorted and taxa, with the exception of Oligochaeta and Marenzelleria spp. polychaetes, were identified to the species level. Taxa were counted and weighed to determine their abundance and biomass (wet mass) per square meter.

Bioturbation Potential (BP C ) and Irrigation Potential (IP C )
To calculate the bioturbation and bioirrigation potential, wet mass (WW) was converted to ash free dry mass (AFDW). The conversion was based on literature data; for bivalves, the coefficients were used for individuals with shells [52][53][54][55]. The Bioturbation Potential Community Index (BP C ) at individual sites was calculated by summing the bioturbation potentials (BP i ) calculated for individual taxa [30,36].
where for taxon i: B i is biomass (in ash free dry mass g·m −2 ) and A i is abundance (ind.·m −2 ) at each sample, while M i , mobility, and R i , sediment reworking, are categorical scores assigned to each species (Table A1). The Irrigation Potential Community Index (IP C ) at individual sites was calculated by summing the irrigation potentials (IP i ) calculated for individual taxa [56].
where for taxon i: B i is biomass (in ash free dry mass g·m −2 ) and A i is abundance (ind.·m −2 ) at each sample, while feeding type (FT i ), burrow type (BT i ) and depth (ID i ) are scores for the trait categories assigned to each species (Table A1). Exponent 0.5 used in BP C emphasizes the importance of organisms with high density and relatively low biomass, while exponent 0.75 used in IP C emphasizes the activity of organisms with larger sizes but lower densities [33].

Vertical Distrbution of Macrofauna in Sediment
The analysis of the vertical distribution of macrozoobenthos in the sediment to determine the maximum burial depth of each taxon and the entire community was performed for both the abundance and biomass of organisms from 51 cores. To present the vertical distribution, the benthic macrofauna abundance and biomass measured in separate sediment layers were recalculated per 1 dm 3 volume. The burial depth data were averaged for all cores in which a given taxon occurred. The percentage of individual taxa abundance and biomass (90%) in the studied layers was indicated to determine the typical depth of occurrence.

Data Analysis
Principal Component Analysis (PCA) was carried out to determine the relationship between physicochemical conditions in bottom water and surface sediments, and the variability between the sites. A matrix with normalized data on bottom water temperature, salinity, dissolved oxygen concentration and organic matter content in surface sediments was used in statistical analysis. Environmental parameters were strongly correlated with the depth of the basin i.e., salinity (Pearson's r = 0.95), DO (r = −0.79), type of sediment (r = 0.86) and LOI (0.65).
Prior to biological data analysis, the biomass at each sampling site was averaged and square root transformed. Cluster (Bray-Curtis similarity) and SIMPROF analysis was used to determine the similarity of macrofauna samples. The SIMPER procedure was applied to identify species responsible for similarities/differences in macrozoobenthic communities between the analyzed sites [57]. Biota and Environment matching analysis (BEST BIO-ENV) was performed to determine the effects of temperature, salinity, DO concentration in bottom water and organic matter content in surface sediment on the formation of benthic fauna communities. Distance-based linear models (distLM) were used to examine the effects of environmental variables on biomass, maximum burrowing depth, BP C and IP C [58]. First, the relationships between the variables were examined and oxygen concentration was excluded from the analysis as being strongly correlated with salinity (Pearson's r = −0.84). The following three environmental variables were selected: temperature, salinity and LOI and log(x + 1) transformation was used before analysis. Stepwise selection and the AICc stopping criterion were used in distLM to investigate the role of environmental variables in predicting biological traits of macrozoobenthos. Resemblance matrices were based on the Euclidean distance similarities between the sites. The results of marginal tests indicate the proportion of the variation the predictor accounts for on its own, while the results from the sequential test indicate the proportion added by the predictor to the cumulative total proportion explained. The statistical analyses were computed in PRIMER v6 & PERMANOVA +. Maps with results were prepared using Arc GIS Pro 2.9.0, ESRI Inc., Redlands, California, the United States of America.
Data from individual cores were used to analyze the relationship between the biological parameters. The relationship between bioturbation and bioirrigation potential indices (calculated using WW and AFDW, and two different exponents in the case of IP C ) and the number of taxa, abundance, biomass and maximum burrowing depth were determined by Spearman's rank correlation test. In addition, IP C values obtained when considering the maximum burrowing depth of macrofaunal individuals in the sediment observed in this study were also compared with those assumed based on the literature and expert knowledge. Prior to the statistical analysis, the normality of the data was tested (Shapiro-Wilk test p < 0.05).

Environmental Conditions
Bottom water temperature at the surveyed sites was relatively uniform (below 6.2 • C), except for sites VE03, VE05, and VE18, which were surveyed in the summer season, above thermocline (Table 1). Bottom water salinity was generally higher in the deeper parts and reached 12.7 in the Gdańsk Deep (site TF0233). The opposite situation was observed for dissolved oxygen concentrations in bottom water. Oxygen conditions were above 4.68 mL·dm −3 at the shallow sites, but oxygen deficiency was observed in the deepest part-below 3.41 mL·dm −3 , and the two deepest sites (VE43 and TF0233) showed hypoxia (DO < 2 mL·dm −3 ). Sediment variability was fairly typical for the coastal areas. The shallow sites were characterized by the presence of medium and fine-grained sands, while deeper sites were dominated by clay and silt. PCA analysis was conducted to determine the effect of four physicochemical parameters on the variability between the sites ( Figure 2). The first principal component explains 59.6% (eigenvalue 2.39), and together with the second principal component (eigenvalue 1.08) a total of 86.6% of the total variation (Table 2). Salinity with a coefficient of −0.612 has the largest contribution to the distribution along the PC1 axis. The distribution along the PC2 axis was most significantly affected by bottom water temperature (coefficient 0.939).

Macrofauna
The study revealed the presence of a total of 23 macrofaunal taxa in the Gulf of Gdańsk. Taxa with the highest frequency in the Vistula estuary (above 70%) were the bivalve Macoma balthica, Oligochaeta, the polychaetes Bylgides sarsi, Marenzelleria spp., Pygospio elegans, as well as the crustacean Corophium volutator and the gastropod Peringia ulvae (data not shown). The biodiversity of benthic organisms decreased with depthfrom 16 taxa at site VE05 to no organisms in the Gdańsk Deep. The main factors determining the structure of macrozoobenthos biomass were salinity and oxygen concentration in the bottom water (BIOENV, r = 0.74).
Based on cluster and SIMPROF analysis, three groups of sites were distinguished with respect to the biomass of the identified macrofauna taxa (Figure 3). In both group 1 and group 2, M. balthica was the most dominant species in the biomass and significantly contributed to the similarity of biomass in both groups (contribution to the total biomass of 67% and 79%, respectively). In addition, species that contributed to the similarity in group 1 were Hediste diversicolor (11%), P. ulvae (9%) and Mya arenaria (8%). Other taxa that accounted for the similarity between sites in group 2, in addition to M. balthica, were Marenzelleria spp. (7%) and Halicryptus spinulosus (5%). In addition, group 3 comprised the deepest sites, where only polychaetes represented by the species B. sarsi were observed.

Macrofauna
The study revealed the presence of a total of 23 macrofaunal taxa in the Gulf of Gdańsk. Taxa with the highest frequency in the Vistula estuary (above 70%) were the bivalve Macoma balthica, Oligochaeta, the polychaetes Bylgides sarsi, Marenzelleria spp., Pygospio elegans, as well as the crustacean Corophium volutator and the gastropod Peringia ulvae (data not shown). The biodiversity of benthic organisms decreased with depth-from 16 taxa at site VE05 to no organisms in the Gdańsk Deep. The main factors determining the structure of macrozoobenthos biomass were salinity and oxygen concentration in the bottom water (BIOENV, r = 0.74).
Based on cluster and SIMPROF analysis, three groups of sites were distinguished with respect to the biomass of the identified macrofauna taxa ( Figure 3). In both group 1 and group 2, M. balthica was the most dominant species in the biomass and significantly contributed to the similarity of biomass in both groups (contribution to the total biomass of 67% and 79%, respectively). In addition, species that contributed to the similarity in group 1 were Hediste diversicolor (11%), P. ulvae (9%) and Mya arenaria (8%). Other taxa that accounted for the similarity between sites in group 2, in addition to M. balthica, were Marenzelleria spp. (7%) and Halicryptus spinulosus (5%). In addition, group 3 comprised the deepest sites, where only polychaetes represented by the species B. sarsi were observed. The average dissimilarity between group 1 or group 2 and group 3 was >99%. In both cases, M. balthica accounted for the highest proportion of dissimilarity (>54%). The average dissimilarity between group 1 or group 2 and group 3 was >99%. In both cases, M. balthica accounted for the highest proportion of dissimilarity (>54%).  All the biological parameters studied reached the highest values at the shallow sites and site VE46, and their values gradually decreased in subsequent groups with increasing depth (Table 3). Benthic communities at the shallow and intermediate sites were characterized by high taxonomic diversity of macrofauna. The highest values of density and biomass of macrofauna were observed at the shallow sites and decreased with depth. All the biological parameters studied reached the highest values at the shallow sites and site VE46, and their values gradually decreased in subsequent groups with increasing depth (Table 3). Benthic communities at the shallow and intermediate sites were characterized by high taxonomic diversity of macrofauna. The highest values of density and biomass of macrofauna were observed at the shallow sites and decreased with depth. Similarly, the values of the BP C and IP C indices decreased, with the values of both indices being lower by half at the intermediate sites compared to the shallow sites. Table 3. Number of taxa, maximum burial depth of macrofauna, mean values: abundance, BP C and IP C (min.-max), contribution of individual taxa to the formation of these parameters, and in each group of sites provided in Figure 3.

Group 1 Group 2 Group 3
No. of taxa 11 (6-16) 8 (7-10)  The vertical distribution of organisms in each group differed in terms of both abundance and biomass ( Figure 4). In all groups of sites, the largest number (>62%) of organisms was found in the shallowest layer of sediment. This was also the only layer in group 3 containing organisms. However, the maximum biomass of organisms was observed in the deeper sediment layers-as much as 42% of the biomass at the sites from group 1 was found in the 3-6 cm sediment layer, and in group 2, organisms were found in the shallower layers-almost 60% of the biomass was found in the 1-3 cm sediment layer. This is due to the dominance of M. balthica in the infaunal biomass.  The vertical distribution of organisms in each group differed in terms of both abundance and biomass ( Figure 4). In all groups of sites, the largest number (>62%) of organisms was found in the shallowest layer of sediment. This was also the only layer in group 3 containing organisms. However, the maximum biomass of organisms was observed in the deeper sediment layers-as much as 42% of the biomass at the sites from group 1 was found in the 3-6 cm sediment layer, and in group 2, organisms were found in the shallower layers-almost 60% of the biomass was found in the 1-3 cm sediment layer. This is due to the dominance of M. balthica in the infaunal biomass.  The vertical distribution of organisms in each group differed in terms of both abundance and biomass ( Figure 4). In all groups of sites, the largest number (>62%) of organisms was found in the shallowest layer of sediment. This was also the only layer in group 3 containing organisms. However, the maximum biomass of organisms was observed in the deeper sediment layers-as much as 42% of the biomass at the sites from group 1 was found in the 3-6 cm sediment layer, and in group 2, organisms were found in the shallower layers-almost 60% of the biomass was found in the 1-3 cm sediment layer. This is due to the dominance of M. balthica in the infaunal biomass.
The vertical distribution of organisms in each group differed in terms of both abundance and biomass ( Figure 4). In all groups of sites, the largest number (>62%) of organisms was found in the shallowest layer of sediment. This was also the only layer in group 3 containing organisms. However, the maximum biomass of organisms was observed in the deeper sediment layers-as much as 42% of the biomass at the sites from group 1 was found in the 3-6 cm sediment layer, and in group 2, organisms were found in the shallower layers-almost 60% of the biomass was found in the 1-3 cm sediment layer. This is due to the dominance of M. balthica in the infaunal biomass.  The vertical distribution of organisms in each group differed in terms of both abundance and biomass (Figure 4). In all groups of sites, the largest number (>62%) of organisms was found in the shallowest layer of sediment. This was also the only layer in group 3 containing organisms. However, the maximum biomass of organisms was observed in the deeper sediment layers-as much as 42% of the biomass at the sites from group 1 was found in the 3-6 cm sediment layer, and in group 2, organisms were found in the shallower layers-almost 60% of the biomass was found in the 1-3 cm sediment layer. This is due to the dominance of M. balthica in the infaunal biomass.   The vertical distribution of organisms in each group differed in terms of both abundance and biomass (Figure 4). In all groups of sites, the largest number (>62%) of organisms was found in the shallowest layer of sediment. This was also the only layer in group 3 containing organisms. However, the maximum biomass of organisms was observed in the deeper sediment layers-as much as 42% of the biomass at the sites from group 1 was found in the 3-6 cm sediment layer, and in group 2, organisms were found in the shallower layers-almost 60% of the biomass was found in the 1-3 cm sediment layer. This is due to the dominance of M. balthica in the infaunal biomass.   -max), contribution of individual taxa to the formation of these parameters, and in each group of sites provided in Figure 3. The vertical distribution of organisms in each group differed in terms of both abundance and biomass (Figure 4). In all groups of sites, the largest number (>62%) of organisms was found in the shallowest layer of sediment. This was also the only layer in group 3 containing organisms. However, the maximum biomass of organisms was observed in the deeper sediment layers-as much as 42% of the biomass at the sites from group 1 was found in the 3-6 cm sediment layer, and in group 2, organisms were found in the shallower layers-almost 60% of the biomass was found in the 1-3 cm sediment layer. This is due to the dominance of M. balthica in the infaunal biomass.
Similarly, the values of the BPC and IPC indices decreased, with the values of both indices being lower by half at the intermediate sites compared to the shallow sites.  -max), contribution of individual taxa to the formation of these parameters, and in each group of sites provided in Figure 3. The vertical distribution of organisms in each group differed in terms of both abundance and biomass (Figure 4). In all groups of sites, the largest number (>62%) of organisms was found in the shallowest layer of sediment. This was also the only layer in group 3 containing organisms. However, the maximum biomass of organisms was observed in the deeper sediment layers-as much as 42% of the biomass at the sites from group 1 was found in the 3-6 cm sediment layer, and in group 2, organisms were found in the shallower layers-almost 60% of the biomass was found in the 1-3 cm sediment layer. This is due to the dominance of M. balthica in the infaunal biomass.     -max), contribution of individual taxa to the formation of these parameters, and in each group of sites provided in Figure 3. The vertical distribution of organisms in each group differed in terms of both abundance and biomass (Figure 4). In all groups of sites, the largest number (>62%) of organisms was found in the shallowest layer of sediment. This was also the only layer in group 3 containing organisms. However, the maximum biomass of organisms was observed in the deeper sediment layers-as much as 42% of the biomass at the sites from group 1 was found in the 3-6 cm sediment layer, and in group 2, organisms were found in the shallower layers-almost 60% of the biomass was found in the 1-3 cm sediment layer. This is due to the dominance of M. balthica in the infaunal biomass.

Group 1 Group 2 Group 3
The vertical distribution of organisms in each group differed in terms of both abundance and biomass (Figure 4). In all groups of sites, the largest number (>62%) of organisms was found in the shallowest layer of sediment. This was also the only layer in group 3 containing organisms. However, the maximum biomass of organisms was observed in the deeper sediment layers-as much as 42% of the biomass at the sites from group 1 was found in the 3-6 cm sediment layer, and in group 2, organisms were found in the shallower layers-almost 60% of the biomass was found in the 1-3 cm sediment layer. This is due to the dominance of M. balthica in the infaunal biomass.
M. balthica accounted for the largest proportion of biomass at all but the deepest sites ( Figure 5) (Table A2). The biomass was also composed of Marenzelleria spp., P. ulvae and H. diversicolor. Only epifaunal B. sarsi was observed at the deepest sites. The coastal sites were characterized by the occurrence of taxa such as Marenzelleria spp. and H. diversicolor, which burrow to a depth of 15 cm. With the distance from the Vistula River, fewer taxa were observed burrowing deeper into the sediment.  M. balthica accounted for the largest proportion of biomass at all but the deepest sites ( Figure 5) (Table A2). The biomass was also composed of Marenzelleria spp., P. ulvae and H. diversicolor. Only epifaunal B. sarsi was observed at the deepest sites. The coastal sites were characterized by the occurrence of taxa such as Marenzelleria spp. and H. diversicolor, which burrow to a depth of 15 cm. With the distance from the Vistula River, fewer taxa were observed burrowing deeper into the sediment.
Both the bioturbation potential index and the bioirrigation potential index followed the distribution of biomass, with the highest values in the shallow areas and in vicinity of the Vistula River mouth, and lower values in the deep area and no bioturbation activity in the Gdańsk Deep. BPC and IPC at all (except the deepest) sites were mainly formed by M. balthica. At the shallow sites, the polychaetes, M. arenaria and Pontoporeia femorata contributed relatively significantly to the formation of BPC, while at VE18 it was mainly formed by Marenzelleria spp. In the formation of IPC, Marenzelleria spp. contributed more than other taxa at several sites (VE18, VE06, VE09). The highest BPC (5001) and IPC (1958) values were recorded at site VE05. Among environmental parameters, salinity (and highly correlated DO) was the most important predictor, explaining more than 44% of the variability in biomass, burial depth, BPC and IPC (Table 4). Salinity (and highly correlated DO) and temperature, and in the Both the bioturbation potential index and the bioirrigation potential index followed the distribution of biomass, with the highest values in the shallow areas and in vicinity of the Vistula River mouth, and lower values in the deep area and no bioturbation activity in the Gdańsk Deep. BP C and IP C at all (except the deepest) sites were mainly formed by M. balthica. At the shallow sites, the polychaetes, M. arenaria and Pontoporeia femorata contributed relatively significantly to the formation of BP C , while at VE18 it was mainly formed by Marenzelleria spp. In the formation of IP C , Marenzelleria spp. contributed more than other taxa at several sites (VE18, VE06, VE09). The highest BP C (5001) and IP C (1958) values were recorded at site VE05.
Among environmental parameters, salinity (and highly correlated DO) was the most important predictor, explaining more than 44% of the variability in biomass, burial depth, BP C and IP C ( Table 4). Salinity (and highly correlated DO) and temperature, and in the case of BP C also LOI, explained more than 80% of the data variation in BP C and IP C . There are strong positive correlations between bioturbation and bioirrigation potential indices, as well as between them and key characteristics of benthic communities, i.e., the number of taxa, abundance, biomass as well as maximum burrowing depth (Table 5). Similarly strong and significant relationships exist between BP C and IP C calculated from differently presented biomass data (WW and AFDW), as well as when comparing IP C calculated from benthic fauna burial depth data obtained in this study with the index using the literature data.
The examination of the macrofauna in different layers of the cores showed that the sediments are inhabited to a depth of 15 cm ( Figure 6). All the studied taxa, with the exception of H. spinulosus, are observed in the shallowest layer of sediment. For some taxa (Planaria torva, Ecrobia ventrosa, Potamopyrgus antipodarum, Saduria entomon, Mysis mixta and Neomysis integer), this is the only layer of occurrence. Few-especially polychaetes-are observed in the deeper sediment layers, and their dominant abundance and biomass occurs in the 3-6 cm and 6-10 cm layers. The deepest recorded taxa are the polychaetes Marenzelleria spp. and H. diversicolor, the clams M. balthica and Oligochaetes, which were found in the layer up to a maximum of 15 cm deep into the sediment. Table 5. Spearman's correlation coefficients of biological parameters from all sampling sites. Both indices were calculated using different exponents (0.5 or 0.75) and two types of animal biomass-wet formalin mass (WW) or ash free dry mass (AFDW). Significance level for all correlations was p < 0.000001.

Conditions of Bottom Water and Sediments and Their Impact on Macrozoobenthos
Coastal areas, such as estuaries, lagoons and bays, are dynamic environments with gradients of freshwater and seawater flows, representing transition zones between land and sea [3,59,60]. Gradients in physicochemical parameters of bottom water and surface sediments are typical for the Gulf of Gdańsk [41,42,61] [this research]. As the depth of the basin increases, salinity increases, DO decreases, while the proportion of the finest fraction, organic matter content and hydrogen sulfide concentration increases. The area of the Vistula outflow is characterized by the presence of increased amounts of organic matter and nutrients supplied with river runoff [62][63][64][65]. Organic matter and nutrients, as well as contaminants on their way from land to open sea are transformed, retained or removed by biota or moved unchanged to the offshore areas of the Baltic Sea [26,44,[66][67][68].
The conditions prevailing in the bottom water and sediments affect the distribution and species composition of macrozoobenthos. In the case of benthic communities inhabiting the seabed of the Gulf of Gdańsk, the factors that had the greatest impact on biomass structure, macrofauna burial depth and indices of bioturbation and bioirrigation potential were conditions such as salinity and oxygen concentration in the water above the bottom, and factors strongly related to these, such as sediment conditions. It is known that as oxygen conditions in the water above the seabed deteriorate, the concentration of toxic hydrogen sulfide in the sediments increases [41,[69][70][71].

Macrozoobenthos
The present study revealed the presence of 23 taxa of the benthic macrofauna in the study area. The results were similar to those obtained during other macrozoobenthos studies conducted in the Vistula estuary [39,41,72,73]. The greatest diversity of benthic organisms was observed in the coastal zone, where the density was dominated by P. ulvae, a gastropod species typical of the coastal zone in the Baltic Sea, while in the deeper zones the species composition of the benthic community shifted and the abundance was dominated by P. elegans and M. balthica, species also common in the Baltic Sea. The biomass in all but the deepest zones was dominated by M. balthica.
The highest number of taxa was observed at some distance from the Vistula estuary (at a depth of 16-24 m). Relatively few taxa were found at the shallowest site (15 m depth), due to the fact that the estuary is highly dynamic and the material carried by the river forms an unstable and easily eroded substrate, unfavourable to macrozoobenthos development [72][73][74][75]. Although organic matter carried with river runoff constitutes food resources for macrofauna [44], it can also cause benthic organisms to become covered and buried, leading in extreme cases to the complete disappearance of benthic macrofauna in a given area [76]. As the depth of the water body increases, both the taxonomic diversity and the biomass of the macrofauna decreases. At the deepest sites, the macrofauna is either absent or represented by single individuals of the surface-living, semi-pelagic polychaete B. sarsi. The reason for this is the decomposition of large amounts of organic matter accumulating on the bottom and stable stratification in the deeper area, which leads to oxygen deficiency or anoxia at the bottom and occurrence of hydrogen sulfide in the surface sediments [69]. These conditions adversely affect the behaviour, physiological processes, fitness of the benthic fauna, and consequently lead to a loss of functions performed by the benthic fauna [36,77,78]. Such a loss of biodiversity can result in reduced resistance of the environment to stress [79].

Bioturbation and Bioirrigation
The research carried out has shown that while the zoobenthos biomass in the Vistula estuary is completely dominated by M. balthica, the use of bioturbation and bioirrigation potential indices reveals the role of other species, i.e., those whose biomass is not large but it is known from experimental studies that their activity can significantly affect biogeochemical processes [80,81]. The benthic communities described in this study are characterized by their high bioturbation and bioirrigation potential in the coastal region. This is where their impact on various compounds is most likely to be greatest. M. balthica, whose intensive bioturbation and bioirrigation activity is relatively well studied, had the largest contribution to the indices [80,81]. Polychaetes of the genus Marenzelleria also contributed relatively significantly to the bioirrigation potential index. Experimental studies have shown that this species is an extremely effective bioirrigator and bioturbator [6,48,82,83]. In situ experiments in the Vistula plume showed a significant increase in nutrient fluxes from sediments inhabited by macrofauna, with the greatest impact observed in the presence of polychaetes [39]. In previous studies, a comparison between bioturbation and bioirrigation potential indices maps showed a very similar pattern, but also some differences [56,84]. For example, differences on a spatial scale were found in the German Bight, with higher IP C scores in areas where sessile or semi-sessile species (i.e., Lanice conchilega and Notomastus latericeus) were particularly abundant [56]. In the Vistula estuary, such a difference is apparent only for one site (VE18), where higher IP C values compared to BP C are due to the abundance of Marenzelleria spp.
The present study demonstrated a strong positive relationship between the two indices and their strong correlation with both the maximum burrowing depth of macrozoobenthos, the number of taxa, as well as the abundance and total biomass of macrozoobenthos. Interestingly, there was virtually no difference in these relationships regardless of how the calculations were made (i.e., wet or ash free dry mass). In an earlier study conducted in another region of the Baltic Sea, the authors found no relationship between the bioirrigation index and the number of taxa [84]. According to Queirós et al. [38], BP C was found to be a good predictor of bioturbation distance (average distance travelled by a sediment particle). However, it was found unsuitable for determining other attributes of infauna, such as bioturbation activity, bioturbation depth or diffusion transport. In addition, the index also appears to be a better predictor of community-level estimates, rather than those for individual species. Statistical models using experimental results showed that BP C explained a considerable amount of variance in oxic processes, i.a. oxic mineralization, total N mineralization, and nitrification [85]. Few studies have also been conducted to determine the correlation between bioirrigation potential index values and actual bioirrigation. However, the results of these studies are inconclusive and require further research. A study by De Borger et al. [86] showed that IP C correlates more strongly with burrow ventilation depth than with ventilation rate. The correlation between IP C and irrigation rate was not confirmed by Toussaint et al. [85].
The present study did not use the bioirrigation index (BIP C ) proposed by Renz et al. [32], the scoring system of which additionally takes into account the distinction between the advection and diffusion system performance. The use of this index would result in higher values for free living species and species living in burrows as well as facultative deposit/suspension feeders in advective sediments. Furthermore, it would result in even higher values of bioirrigation potential in the coastal zone, where the advection system dominates, and an even higher proportion of M. balthica or polychaetes H. diversicolor and Marenzelleria spp. in the index for this zone. At the deeper sites where diffusive sediments occur, bioirrigation potential would be much lower than in the coastal zones. The system by Renz and co-workers [32] would emphasize the variability of the bioirrigation index in the Gulf of Gdańsk and the gradual loss of this function in the environment with increasing depth of the water body. Both approaches to the determination of the bioirrigation potential index are certainly worth testing in further studies, especially those combining studies of benthic assemblages, including functional indices, with experimental studies of the impact of macrofauna on biogeochemical processes, or measurements of animal activity.
The indices used provide only a simplified approximation of the potential capabilities of benthic communities. Bioturbation and bioirrigation are dynamic and complex activities performed by those organisms. They are determined by a number of factors that affect the biological functions of these animals. BP C was observed to follow the seasonal pattern in seawater temperature, with the highest values in summer and autumn [38]. However, it should be kept in mind that temperature and food availability have the potential to impact bioturbation and bioirrigation intensity, as these factors affect physiological processes of benthic species. Studies conducted on the polychaete Alitta virens showed that sediment reworking processes could be affected by both low and high temperature, with the lowest bioturbation intensity under low temperature [87]. Oxygen depletion may also change the activity of animals in the sediment, thus affecting bioturbation and bioirrigation. Depending on the oxygen concentration and exposure time, these conditions can result in, for example, an increase in burrow ventilation, a decrease in animal activity or no activity at all [88][89][90].

Burrowing Depth
The burrowing depth of organisms provides, among other things, an indication of the depth to which they can affect the conditions and processes in the sediments. In the present study, most of the organisms (>62% of all individuals) inhabited the shallowest layer of sediment (0-1 cm). The maximum biomass of organisms can be found in the deeper layers of sediment-deeper layers (3-6 cm) at the shallowest sites and slightly shallower layers at the intermediate sites (1-3 cm). A similar distribution of organisms deep into the sediment was observed in earlier studies conducted in the Gulf of Gdańsk-the highest abundance of organisms was found in the shallowest layer of the sediment and it decreased with depth [40,41]. In contrast to the abundance, the biomass of organisms in the shallow water zone did not decrease with depth and its distribution was more varied-the highest biomass was usually observed in deeper layers, i.e., up to 6 cm into the sediment [40].
However, even organisms living on the sediment surface can play an extremely important role by being active in disrupting the diffusive boundary layer, which improves the oxygen conditions of the sediment [91]. Few organisms, i.e., bivalves and polychaetes, burrow naturally into the sediment and are rarely present on its surface [83,92], and their typical burrowing depth is 3-10 cm. The maximum depth of occurrence of a given taxon depends on the ability of the organism to contact the sediment surface, for example, the burial depth of M. balthica depends on the length of the clam's siphon, which is often also related to the size and age of the organism [46]. Our research showed the occurrence of M. balthica below a depth of 10 cm, which is also the maximum depth at which the bivalves bioturbate and bioirrigate the sediments. Other deep burrowing species-from the genus Marenzelleria-were found up to a depth of 15 cm, but some scientists indicate that these species can burrow as deep as 35 cm [93]. These deep burrowing organisms, such as polychaetes, form burrows that enable water transport in the sediment and aerobic chemical reactions in the deeper layers, as well as affect nutrient cycling [83,94,95].
The vertical distribution of organisms is determined by environmental factors. Organisms change the depth of their occurrence seasonally [46], e.g., M. balthica has been shown to burrow deepest in winter and remain shallowly buried in the sediment during the summer season. Oxygen deficiencies and hydrogen sulfide cause the animals to move to the sediment surface or they become periodically inactive [92,96,97]. While animals are present in the sediment, their functions may be temporarily impaired.

Conclusions
Our research has shown changes in the structure and functioning of benthic communities with increasing distance from the Vistula River mouth. Coastal zones are characterized by relatively high biodiversity and great burrowing depth of macrofauna, as well as high bioturbation and bioirrigation potential of benthic communities. However, this activity disappears in deep zones with the absence of benthic organisms. The lack of bioturbation and bioirrigation means there is no support for biogeochemical transformation by the macrofauna in the deep zones. In the study area, only a few species drive bioturbation and bioirrigation-the bivalve M. balthica and the polychaetes H. diversicolor and Marenzelleria spp. Other taxa had a marginal impact. Such a strong dominance of single taxa in performing bioturbation and bioirrigation could lead to instability in ecosystem functioning in the case that these organisms were to disappear as a result of an ecological disaster, envi-ronmental degradation or disease. At the same time, these large organisms were the only taxa burrowing deep into the sediment (below 10 cm), and thus the only ones supporting geochemical processes deep in the sediment. To summarize, the coastal zone, unlike the offshore zone, proved to be a hotspot for bioturbation-and bioirrigation-driven processes, which are responsible for the proper functioning of the seafloor and basin. However, very poor functional diversity of the benthic macrofauna in the deepest zones means that we should appreciate and protect coastal zones more efficiently. Funding: This study was funded by the BONUS COCOA project, which was supported by BONUS (Art 185), jointly funded by the EU and the National Center for Research and Development (NCBiR, Poland) and Norway Grant (2020/37/K/NZ8/02782 "Benthic coastal buffers against climatic and eutrophication extremes, BUFFER").

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement: Not applicable.
Acknowledgments: We would like to thank Maren Voss for her excellent cruise leadership, scientific and operational crews and captains of the ships, especially Michael Pötzsch, Anna Borecka, Heather Reader, Natalia Kozak and Monika Lengier for their help during sampling. We also thank Ines Bartl for providing data on sediment parameters from two sites in 2016.