*3.2. Macrozoobenthic Assemblages*

In the present study, we found 43 taxa which marked a significant increase in benthic species richness and diversity in the Yundang Lagoon compared to previous years and decades (Table 4). In particular, according to the available published and gray literature, no species were present in the 1980s, a single polychaete, *Neanthes succinea*, was found in the 1990s, and a few more species (n = 9) were reported during the 2000s, with two polychaetes dominating (*Neanthes glandicincta* and *Capitella capitata*), and the first report of the invasive bivalve *Mytilopsis sallei* (Table 4). In our study, *M. sallei* was the most abundant benthic species, accounting for 75% of the total individuals and 88% of the total biomass. Several other less abundant species included the gastropods *Stenothyra glabra* and *Rissolina plicatula*, the bivalve *Pseudopythina tsurumaru*, the crustaceans *Corophium* sp. and *Corophium uenoi*, and the polychaete *Cossurella dimorpha* (Table 4). All these species, except *C. dimorpha*, were found at site A. *M. sallei* was dominant at all stations and replicates of sites A and C, very abundant at a single station of site F, and was not found at sites B, D, and E. Another invasive species, the ascidian *Styela plicata*, was found only in some replicates of site C, but followed *M. sallei* in terms of its contribution to the total biomass (8.4%). *S. glabra, P. tsurumaru* and *N. glandicincta* occurred at all sites except B and F, while *C. dimorpha* was dominant in B and absent in A and C. Finally, different subsets of species were exclusive to a single site: nine species in A, four in B, four in C, nine in D and only two in E (Table 4).

The PERMANOVA conducted on the abundance of all species revealed significant differences among sites. Site A was distinct from all the others, while site C differed from B, D, and E, and site F from D and E (Table 5). These differences were similar to those found for *M. sallei* which showed the highest abundance and biomass at site A. Significant differences also emerged from the structural biotic parameters, including total abundance, species richness, diversity, and biomass (Table 6, Figure 6). In the pairwise comparisons, site A was significantly different from most of the other sites in terms of total abundance and species richness, and less so in terms of diversity (- from F) and biomass (- from B, D, and E) (Table 7). Results sugges<sup>t</sup> that the dominance of *M. sallei* appeared to be one of the discriminating factors responsible for the differentiation of sites and their species richness. In addition, site A had the highest species richness, followed by C, D, and E, and was significantly different in terms of diversity from site F which had the least number of species (Table 7, Figure 6).


*Water* **2019**, *11*, 1692

**Table 4.** Changes in species number and composition of macrozoobenthos collected from the bottom sediment in the Yundang Lagoon in the years before the


**Table 5.** Results of PERMANOVA (based on Bray–Curtis similarity index) testing di fferences among sites (A, B, C, D, E, F) for the abundance of all the taxa. Comparisons of post−hoc test between all pairs of sites are also given. Significant comparisons at *p* < 0.05 are marked by \* and at *p* < 0.01 by \*\*. The Bonferroni correction was used. Permutation number = 9999.

**Table 6.** Results of PERMANOVA (based on Euclidean distance) testing di fferences among sites for abundance, number of species, Shannon diversity (H'), and biomass. Permutation number = 9999. All values are significant at *p* < 0.0005.


**Figure 6.** Mean values (n = 9; ±SE standard error) of the total number of individuals (N, ind. 600 cm<sup>−</sup>2), total number of species (S), Shannon index (H'), and biomass (g.wet weight 600 cm<sup>−</sup>2) at each site.


**Table 7.** Results of post−hoc test between all pairs of sites. (A–F). Significant comparisons at *p* < 0.05 are marked by \* and at *p* < 0.01 by \*\*. The Bonferroni correction was used.

Di fferent aspects of diversity, i.e., species richness and the even distribution of individuals over the species, di fferently a ffected the benthic distribution patterns of the study sites, as shown by the cumulative dominance curves (Figure 7). The most even distribution in the number of individuals among species was exhibited by site B (low species richness and low diversity), D (high species richness and highest diversity), and E (intermediate between B and D in terms of species richness and diversity). Instead, the highest abundances of *M. sallei* mostly a ffected the dominance curve of sites A and C notwithstanding their high species richness (highest at site A). Finally, the short and flattened curve for site F reflected both the lowest number of species (n = 5) observed and the high number of *M. sallei* individuals which occurred at one of its three stations.

**Figure 7.** Cumulative dominance curves on macrobenthic abundance data from each sampling site.

The separation of sites based on the macrobenthos structure was revealed by the MDS (Figure 8). The low level of stress (=0.05) confirmed that there was a good representation of the community pattern. Site A was most markedly distinguished from all the other sites, with the exception of one station of site F dominated by *M. sallei*, while both sites A and C were clearly separated from the other sites. The significance of di fferences among the sites was confirmed by the PERMANOVA (Table 5). Canonical correspondence analysis (CCA) revealed the relationship between the selected environmental variables and the ordination model produced (Figure 9). The first two axes accounted for 69.1% of the explained variability (51.5% and 17.6% for the *x*- and *y*-axes, respectively). The environmental variables were plotted as correlations with site scores and the strongest correlations were found for sand (negative) and for SS, Md, nitrate, and RP (positive); TOC was less correlated with the axes. Overall, the environmental variables mostly a ffecting the separation of sites were those connected with the trophic and hydrodynamic features of the inner lagoon sites (D, E), with sand vs. Md and SS in opposite directions to indicate, as proxies, higher and lower hydrodynamics, respectively (Figure 9).

**Figure 8.** Non-metric multidimensional scaling (nMDS) ordination model of log transformed abundance data (stress: 0.05). Each symbol represents a sampling station at each site obtained as the mean of three replicate samples.

džŝƐ ϭ

**Figure 9.** Plot of Canonical correspondence analysis (CCA) performed on species abundance and environmental variables most correlated with two axes and indicative of trophic (i.e., ammonia, nitrate, reactive phosphorus, and TOC) and hydrodynamic (i.e., suspended solids, sand, and size-particle) conditions.
