## 2 Materials and Methods
### Statistical Analysis  
The structure of the macrobenthic community was analyzed in terms of species richness (S), Shannon-Wiener diversity (H\({}^{\prime}\)), Margalef index (d), total abundance (A), and biomass. Macrobenthic community composition was calculated using similarity percentage analysis (SIMPER) and analysis of similarity (ANOSIM) in PRIMER 7 and PERMANOVA+ software package, version 7.0.17. Multivariate analysis (Bray-Curtis similarity) was performed with fourth square-root transformed data using PRIMER 7 and PERMANOVA+ software package, version 7.0.17. Abundance data (transformed with log (x + 1) function) were used to create a similarity matrix based on the Bray-Curtis similarity index. The matrices were used to perform hierarchical cluster analysis (CLUSTER) and nonmetric multidimensional scaling ordination (MDS), in order to evaluate similarities between samples (grouped by stations). To evaluate the contribution of the species that determined the differences between the groups created by the cluster analysis, the similarity percentages were calculated (SIMPER analysis). We used a one-way analysis of similarity (ANOSIM) to test the null hypothesis of no indices differences among all sampling sites. The statistical coefficient Global R-value highlights variable separations between the sites with R-value between 0 to 1. If the ANOSIM R statistic is close to its maximum range of 1 that shows that there is a clear separation of the habitats.  
These data were processed with Ocean Data View, version 5.4.0 [46]. The AZTI Marine Biotic Index, AMBI [47], and the multivariate AMBI, M-AMBI [48] were calculated using the freeware program available on www.azti.es (accessed on 21 March 2023). Based on the species list provided in AMBI v6.0, we divided the macrobenths collected in the study areas into five ecological groups (EG) ranked according to their sensitivity to a progressive gradient of stress (EG I-sensitive, EG II-indifferent, EG III-tolerant, IV second-order opportunistic, and EG V-first-order opportunistic). The percentage of macrozoobenthic species belonging to the different ecological groups was then calculated for each station. AMBI results vary from 0 (high environmental quality) to seven (extremely polluted environment). M-AMBI was then derived based on a factor analysis of AMBI, species richness and H\({}^{\prime}\)[48]. M-AMBI\({}^{\prime}\)(n) was calculated as the arithmetic mean of the min-max normalized AMBI, H and S [49]. It was proposed as one of the indicators for assessing the good environmental status (GES) of marine habitats in Romanian marine waters. The threshold values of the AMBI, S, M-AMBI\({}^{\prime}\)(n) indices were used according to [PERSON] et al. (2020) [9]. Therefore, GES for the three habitats was established according to the thresholds of M-AMBI\({}^{\prime}\)(n), 0.68--circalittoral mud dominated by _Melinna palmata_, 0.71--circalittoral mud with _Abra nitida_, _Pitar rudis_, _Spisula subtruncata_, _Acanthocardia pauccostata_, _Nephthys hombergii_, and _Mytilus galloprovincialis beds_ and 0.60--circalittoral mud and sandy mud with _Upogebia pusilla_.  
In this paper, no relevant statistical correlations could be established between the biological, geophysical, and fishing data, as the information extracted were retrieved from several non-targeted projects and research, which did not overlap exactly.