A Comparative Approach to Detect Macrobenthic Response to the Conversion of an Inshore Mariculture Plant into an IMTA System in the Mar Grande of Taranto (Mediterranean Sea, Italy)

Abstract

The expected bioremediation effect, driven by the conversion of an inshore mariculture plant into an Integrated Multi-Trophic Aquaculture (IMTA) system, which could mitigate the fish farm impact, related to the accumulation of organic matter on the seabed, has been studied. The ecological quality status was studied following a Before-After-Control-Impact (BACI) design and variation measured through M-AMBI and compared with the results of univariate and multivariate analyses of variance, to evaluate the sensitivity of the two methodologies. Results from M-AMBI indicated a sharp change in the ecological quality status, just after one year of the conversion of the plant. By contrast, although changes were detected also utilizing univariate and multivariate statistical analysis, the natural temporal variability characterizing the area partially masked evidence of environmental amelioration.

1. Introduction

Finfish mariculture has expanded significantly in recent years, due to rising demand for seafood and the generally stable production from wild capture fisheries [1,2]. Consequently, there has been increased attention on developing innovative techniques to control and minimize the environmental impact of this activity [3,4]. Indeed, intensive mariculture of high-value carnivorous organisms can frequently lead to several environmental and sustainability issues [5,6,7,8], mainly referring to the release of organic and inorganic substances as both particulate waste and dissolved nutrients in the water column [9,10,11].

Ecosystem degradation and algal blooms driven by water eutrophication are thus potential outcomes of inshore plant activities located in enclosed environments, representing a major threat to the tourism sector, which is severely competing for space with aquaculture in Mediterranean coastal waters [12]. The considerable nutrient load spilled from the fish farms can negatively affect the benthic environment due to increased organic enrichment, leading to alterations in sediment chemistry and benthic biodiversity [10,11,13,14,15].

One possible solution to reduce the environmental impact of marine aquaculture, particularly for inshore plants, is the use of Integrated Multi-Trophic Aquaculture (IMTA) [16]. The IMTA technique consists of combining the cultivation of fed aquaculture species with additional commercially relevant organisms (bioremediators) capable of extracting organic and/or inorganic compounds from the seawater, assimilating aquaculture-derived waste to build their own biomass [17,18,19,20]. Hence, this practice could potentially biomitigate the negative environmental impacts of aquaculture, whilst simultaneously providing a possible economical return for the farmers [17,21,22,23,24,25,26].

Benthic assemblages are effective indicators of the health status of coastal marine environments [27,28,29], which have been extensively used to assess the impact of aquaculture [30,31,32,33,34], as well as to detect their recovery after the cessation of the activity [35]. However, the possible effect of IMTA practice has never been measured in the field by monitoring changes in benthic assemblages. Only measurements of organic matter have been performed to test the action of the IMTA, indicating that water quality can be sustained by IMTA activities [36,37].

The present paper refers to an innovative inshore IMTA rearing model, performed at a preindustrial level within the EU Remedia Life project (LIFE16 ENV/IT/000343) in the Gulf of Taranto (Ionian Sea), where a new set of filter feeder bioremediators, such as polychaetes, sponges, mussels coupled with macroalgae, have been reared within a fish farm for the first time at the European level [18]. One of the major novelties of this project was the experimentation of artificial vertical collectors placed in the water column for enhancing the natural settlement of extractive sessile macroinvertebrates, giving a good performance in terms of production [18]. An accurate knowledge of the area was required before implementing an integrated approach in IMTA systems, including the dynamic of the local fouling assemblages [38,39]. Hydrodynamics is also an important factor that allows to optimize mitigation efforts where plant waste accumulates [40].

The area of interest of the project is located in a confined environment (i.e., Mar Grande) that does not allow high dilution and dispersion of the wastes far from the cages, so that organic loads are expected to accumulate in the nearby sediments. Hence, the first action of the project has been a survey to evaluate the environmental conditions below and around the cages before the start of the project. In addition, the assessment of the ecological status allowed to individualize the best area where waste accumulated and in turn, where the bioremediating organisms could have been placed [41]. Indeed, the ex-ante survey highlighted a local impact of the farm activity mostly affecting the area located northwest with respect to the aquaculture plant [41].

Changes in biological assemblages along natural gradients or following the effect of anthropogenic perturbations can be detected by applying proper experimental designs that consider natural spatial-temporal variability and include an appropriate replication of reference conditions. Based on this principle, the BACI (Before-After-Control-Impact) assessment of the status of the macrobenthic assemblages represents one of the most effective statistical tools for analyzing the effect of drivers of change on the natural systems [42,43,44,45]. However, the BACI approach in the detection of anthropogenic changes in assemblages may be insufficient since natural variation between the two sites (impacted vs. control) might be interpreted as human disturbance [46]. For these reasons, asymmetrical designs including one impacted vs. multiple control locations analyzed before and after the impacting event (beyond-BACI designs) were proposed as tools to mitigate the effect of natural variation on the analyzed system [46].

Such kind of studies, following the logic of the hypothetic deductive method, commonly make use of proper analytical tools, such as the Analysis of Variance (ANOVA) or its multivariate counterpart PERMANOVA (Permutational Multivariate Analysis of Variance) to detect changes in assemblages between locations before and after the perturbation occurs [44,47].

On the other hand, biotic indices, which reflect the quality of the environment, are widely used in the marine realm in the context of the European Union Water Framework Directive (WFD), as well as Marine Strategy Framework Directive (MSFD). Synthetic indices, such as the AZTI’s marine biotic index (AMBI) and the multivariate-AMBI (M-AMBI), which integrates the AMBI, the Shannon-Wiener index of diversity and the species richness, have been extensively applied to detect different impact sources worldwide [48,49,50,51,52,53,54] and successfully utilized to evaluate the effects of aquaculture [31,55,56,57,58].

M-AMBI provides a quantitative–qualitative view of the assemblages, as the species are divided into five ecological groups (EGs) in relation to their different degrees of tolerance to a progressive increase in stress: the lower the EG and the higher the sensitivity of the species [48]. Hence, the level of diversity and abundance of invertebrate taxa and the proportion of disturbance-sensitive taxa, allow to obtain a realistic indication of the quality of the benthic habitats, following the criteria required by the monitoring agencies [52].

The main aim of the present study is to assess, for the first time, the bioremediation effect derived by the conversion of the mariculture plant into an IMTA system, analyzing variations in macrobenthic assemblages before and after the intervention, using a BACI approach. Moreover, we aim to compare the assessment of the quality status based on the synthetic index M-AMBI with the results of univariate and multivariate analyses of variance, to evaluate the sensitivity of the two methodologies in assessing changes in assemblages.

2. Materials and Methods

2.1. Study Area

The study area is located on the southwest coast of the Mar Grande of Taranto (40°25′56″ N; 17°14′19″ E) (Southeast Italy, Ionian Sea) (Figure 1). The Mar Grande of Taranto is a semi-enclosed basin reaching a maximum depth of 45 m, with an average annual water temperature of about 18 °C and salinity of about 38‰, almost uniform over the year. The local surface current is directed from the northeast to the southwest at a speed of about 3 cm s⁻¹. At the bottom, the direction of the current is inverted, proceeding from southwest to northeast at a speed of about 1.3 cm s⁻¹.

The investigation was performed at the aquaculture plant “Maricoltura Mar Grande”, which covers a surface of 0.06 km² and is located about 600 m away from the coast. It consists of 15 cages (Ø 22 m), placed at a depth ranging from 7 to 12 m and producing about 100 tons year⁻¹ of European seabass Dicentrarchus labrax (Linnaeus, 1758) or gilthead sea bream Sparus aurata, Linnaeus, 1758.

2.2. Ex-Ante Analysis

In order to investigate the possible effect of the innovative IMTA system on the surrounding environment, an accurate monitoring survey, including both biological and physico-chemical variables, was performed in the area before the beginning of the experimental activities [41]. This ex-ante survey aimed to identify the most suitable site for the placement of the bioremediating system, but also to assess a reliable reference baseline for evaluating the possible environmental changes after the bioremediation activity. The analysis revealed an area of the plant which was more impacted, and which was chosen as the treatment area.

The benthic assemblages were analyzed in four sites (A, B, C, D) located at the corners of the plant (Figure 2) in order to detect the area of the surrounding bottom mostly affected by the mariculture activities, to be used to test the effect of the bioremediation system. The M-AMBI approach was used to analyze the Ecological Quality Ratio (EQR) at each site.

Following the results of the M-AMBI preliminary analysis (Figure 3), site A revealed an EQR clearly lower than the other stations and it was selected as an impacted site to be treated throughout the bioremediation system (Treatment site), while site B, which showed the highest EQR, was selected as a Control site [41].

Therefore, half the plant was converted into an IMTA system in October 2018. Here, three long lines were realized and placed around six cages. Another six cages, separated by an additional line of cages, acted as the control area (Figure 4).

Within the long lines, the space between two consecutive buoys forms a breeding “chamber”, housing modules for the bioremediating organisms (Figure 4). The IMTA rearing system, using the polychaete Sabella spallanzanii (Gmelin, 1791), the sponge Sarcotragus spinosulus Schmidt, 1862, the mollusc Mytilus galloprovincialis Lamarck, 1819 and the macroalgae Chaetomorpha linum (O.F.Müller) Kützing, 1845 and Gracilaria bursa-pastoris (S.G.Gmelin) P.C.Silva, 1952, is described in Giangrande et al. (2020) [18].

The Remedia Life project had a duration of 4 years, and the present study refers to the first 3 years of activity, from 2018 to 2020, including two productive cycles (2018–2019 and 2019–2020).

2.3. Sampling of the Macrozoobenthos

The seabed under the cages was mainly composed of mud. The only hard substrates were represented by concrete blocks and metal chains for fish cage anchoring. Both soft and hard bottoms were investigated, and the results from hard bottom assemblages are available in Arduini et al. (2022) [38]. For all the sampling campaigns, soft bottom samples were collected using an Ekman grab (152 mm³) and three replicates were taken at sites A (Treatment) and B (Control) in February 2018 and July 2018, before the conversion of the plant; and in July 2019, February 2020 and July 2020 after the bioremediating action using a BACI approach. Nine months elapsed from plant conversion (October 2018) to first sampling (July 2019), so that the bioremediating organisms had time to grow and hypothetically act on the system.

Sediment samples were sieved in the field through a 0.5 mm mesh and the individuals retained were preserved in 90% ethanol. Samples were then washed in the laboratory and stored in 70% ethanol. After sorting, the specimens were identified to the lowest possible taxonomic level and counted.

2.4. M-AMBI

M-AMBI values were obtained using the AMBI Software (version 6.0), which can be downloaded from the freeware http://ambi.azti.es (accessed on 1 June 2022) and using the updated May 2022 species list. M-AMBI integrates the Shannon-Wiener H’ diversity index, the number of species (S) and the AMBI biotic index [48], which is calculated with the following formula:

$$\mathrm{AMBI}=\left[(0\times \%\mathrm{EGI}\right)+\left(1.5\times \%\mathrm{EGII}\right)+\left(3\times \%\mathrm{EGIII}\right)+\left(4.5\times \%\mathrm{EGIV}\right)+\left(6\times \%\mathrm{EGV}\right)]/100,$$

(1)

where EGs I-V represent the ecological groups.

The value of the M-AMBI varies between 0 and 1 and corresponds to the Ecological Quality Ratio (EQR), which is expressed as High-Good-Moderate-Poor-Bad conditions for equal size intervals, and which is calculated starting with the following formula:

$$\mathrm{EcoQS}=\mathrm{K}+\mathrm{a}\times \mathrm{AMBI}+\mathrm{b}\times {\mathrm{H}}^{\prime }+\mathrm{c}\times \mathrm{S},$$

(2)

where a, b, c, and K are the discriminant coefficients related to the typology of the water body, and the EcoQS is the Ecological Quality Status. The comparison of monitoring results with the reference conditions allows to derive the EQR [49,52].

The reference conditions and the type-specific EQR for the application of M-AMBI were obtained from the Regional Agency for Environmental Protection report [59] (related to the Apulia region, where the study area is located) and are shown in

Table 1. Reference conditions and the type specific EQR obtained from the 2016–2018 Regional Agency for Environmental Protection report.

Reference Conditions			EQR
AMBI	H’	S	High/Good	Good/Moderate
0.5	4.8	50	0.81	0.61

.

2.5. Univariate and Multivariate Analyses

The design included the factor “Site”, fixed, with two levels (i.e., the Treatment and the Control site); the factor “Before versus After”, fixed, with two levels, crossed to Site; and the factor “Time”, random, with five levels, nested in “Before versus After” and crossed to Site. Although a Beyond-BACI design including appropriate spatial replication (i.e., multiple controls) would have allowed the detection of the causal relationship between human activities and environmental changes, this design still allowed to detect sustained temporal patterns of difference between the control and the putatively bioremediated site, if any of the F-tests for the interaction terms [Site x Before vs. After] or [Site x Time (Before vs. After)] was significant [44,60].

To evaluate eventual differences in the temporal patterns of the assemblage structure, diversity and evenness between the bioremediated and the control location, multivariate and univariate analyses were used. Multivariate permutational analysis of variance (PERMANOVA) [47] used Bray–Curtis similarity matrices [61] of square root transformed abundance of taxa, to reduce the influence of most abundant species, with 9999 permutations of residuals under a reduced model or using Monte Carlo random draws [62] from the asymptotic permutation distribution when too few permutations were available for a given test (i.e., pairwise comparisons a posteriori). Even if Time was a random factor, performing pairwise comparisons for pairs of sampling times (Before or After) across locations, allowed investigating if and how patterns of temporal variability differed between the control and the potentially bioremediated sites. Homogeneity of multivariate dispersions for the interaction terms was verified with PERMDISP analysis [63].

Non-metric multidimensional scaling (nMDS) based on centroids was used as a graphical ordination of the data.

Species total abundance (S), Margalef’s diversity index (d), and Pielou’s evenness index (J’) calculated with the DIVERSE routine on untransformed data were analyzed by means of PERMANOVAs on Euclidean distances, including the same factors described for the multivariate analysis. Cochran’s C-test was used before each analysis to check for homogeneity of variance [64].

The Similarity Percentages analysis (SIMPER) [65] was used to evaluate the species contributing most to dissimilarity between the two sites before and after the conversion of the plant.

All multivariate tests and univariate PERMANOVA were run in the software package PRIMER-E v7 [66] with the PERMANOVA extension [67].

3. Results

3.1. M-AMBI

In total, 3422 organisms were identified at the finer taxonomical level (i.e., identification at the species level was possible for 97% of the individuals) and classified into 235 different taxa (Appendix A).

Before the implementation of the IMTA system, the Treatment site was characterized by a low diversity and the presence of opportunistic species that tolerate high organic loads, such as Capitella capitata. Nine months after the conversion of the plant, the assemblage composition at this site changed: C. capitata population decreased and was replaced by species indicating good environmental conditions such as Spiochaetopterus costarum (Appendix A).

Although species richness (S) varied only seasonally at both stations, the diversity index H’ values increased over time at the Treatment site, with the highest value in February 2020 (

Table 2. Richness (S), diversity (H’) and Ecological Quality Ratio (EQR) of the two sampling sites over time.

February 2018 (T1)				July 2018 (T2)
Site	S	H’	EQR	Site	S	H’	EQR
Treatment	44	2.45	Moderate	Treatment	31	2.42	Moderate
Control	68	4.93	High	Control	68	5.19	High
July 2019 (T3)				February 2020 (T4)
Site	S	H’	EQR	Site	S	H’	EQR
Treatment	36	3.90	Good	Treatment	78	4.83	High
Control	57	4.91	High	Control	94	5.18	High
July 2020 (T5)
Site	S	H’	EQR
Treatment	36	4.19	Good
Control	65	4.61	High

). At the Control site, species diversity remained almost unchanged (Table 2).

In July 2019, M-AMBI values improved at the Treatment site, whose rating changed from “Moderate” to “Good” (Table 2). In February 2020 there was the most evident change in M-AMBI values, with the Treatment site obtaining the EQR rating “High”, such as the Control site, while in July 2020 EQR was similar to July 2019, with the Treatment site rated as “Good” and the Control site as “High” (Table 2).

Overall, M-AMBI values showed a clear improvement of the treated site after the conversion of the plant to IMTA (Figure 5), leading to a positive change in its EQR.

3.2. Univariate and Multivariate Analyses of Variance

PERMANOVA on the structure of benthic assemblages revealed significant differences between the bioremediated and control site in the patterns of temporal variability before and after the conversion of the plant to IMTA system (S × T(B),

Table 3. PERMANOVA on benthic assemblages of Treatment and Control sites at Before (T1, T2) and three After (T3–T5) the bioremediating action. Abundances of taxa were square root transformed.

Source of Variation	df	MS	Pseudo-F	P(perm)	perms	P(MC)	Denominator
Site (=S)	1	9072.2	2.694	0.088	9189	0.012	S × T(B)
Before vs. After (=B)	1	13110	1.840	0.099	10	0.091	T(B)
Time(B) (=T)	3	7125.1	3.978	0.0001	9861	0.0001	Residual
S × B	1	6563.1	1.949	0.159	9236	0.055	S × T(B)
S × T(B)	3	3368	1.880	0.0001	9835	0.002	Residual
Residual	20	1791

). This effect was not due to differences in the heterogeneity of assemblages (PERMDISP on S × T(B): F_9.20 2.6022, P = 0.4296).

The nMDS showed that despite the temporal trajectory had a similar shape for both the Control and the Treatment sites, the latter showed a lower similarity (higher segregation) in the structure of benthic assemblages across time (Figure 6).

The pairwise comparisons for the factor S × T(B) for pairs of levels of factor Site (Control vs. Treatment) revealed significant differences in the structure of the assemblage between the Treatment and the Control site in July (T2, T3, T5) but not in February (T1, T4), both before and after the implementation of the IMTA system (

Table 4. Pairwise PERMANOVA on benthic assemblages of Treatment and Control sites at Before (T1, T2) and three After (T3–T5) the bioremediating action.

Pairwise of S × T(B) for pairs of levels of factor Site (Control vs. Treatment)
	Time	t	P(perm)	perms	P(MC)
Before	T1—February	1.5256	0.1032	10	0.1060
	T2—July	1.7407	0.0980	10	0.0526
After	T3—July	1.7895	0.0995	10	0.0537
	T4—February	1.4622	0.1003	10	0.1147
	T5—July	1.8557	0.0969	10	0.0509
Pairwise of S × T(B) for pairs of levels of factor Time (Before vs. After)
Time	Site	t	P(perm)	perms	P(MC)
Before T1 vs. T2	Treatment	1.1055	0.2890	10	0.3348
	Control	1.4358	0.1016	10	0.1280
After T3 vs. T4	Treatment	2.2320	0.1046	10	0.0185
	Control	1.9917	0.1034	10	0.0239
After T3 vs. T5	Treatment	1.7216	0.0984	10	0.0626
	Control	0.0995	0.0995	10	0.0248
After T4 vs. T5	Treatment	1.7409	0.0997	10	0.0595
	Control	2.2759	0.0994	10	0.0187

).

The pairwise comparisons for the factor S × T(B) performed for checking differences between sampling times for pairs of levels of factor Site (Control vs. Treatment site) revealed no difference between T1 and T2 (Before) both for the assemblage of the Treatment site and the Control one. After the conversion of the IMTA system (T3, T4, T5) the structure of the assemblage at the Control site significantly varied among sampling times, whilst the assemblage of the Treatment site significantly varied between T3 and T4 but not between T3 and T5, nor between T4 and T5, contrary to what observed by to nMDS representation of the sites across time (Figure 6). The pattern derived from the pairwise comparison suggests the presence of high natural variability in the Control site and the progressive homogenization of the assemblages of the Treatment site after bioremediation. In general, differences in the structure of the assemblage between Control and Treatment sites appear more evident in July compared to February.

The ANOVA carried out on species richness (S) and Margalef’s diversity index (d), showed significant differences between sites and among times, but no differences between the Control and Treatment site among samplings before and after the implementation of the IMTA system were detected (

Table 5. PERMANOVA on species richness (S) of benthic assemblages of Treatment and Control sites at Before (T1, T2) and three After (T3–T5) the bioremediating action. Cochran’s C-test was used before each analysis to check for homogeneity of variance.

Source of Variation	df	MS	Pseudo-F	P(perm)	perms	P(MC)	Denominator
Site (=S)	1	1596.100	32.402	0.0150	9117	0.0099	S × T(B)
Before vs. After (=B)	1	314.690	0.527	0.9010	7	0.5230	T(B)
Time(B) (=T)	3	596.930	11.531	0.0002	9944	0.0003	Residual
S × B	1	15.022	0.305	0.6150	8991	0.6159	S × T(B)
S × T(B)	3	49.259	0.952	0.4320	9941	0.4310	Residual
Residual	20	51.767
Cochran test	C = 0.259		P = 0.6695

and

Table 6. PERMANOVA on Margalef’s diversity index (d) of benthic assemblages of Treatment and Control sites at Before (T1, T2) and three After (T3–T5) the bioremediating action. Cochran’s C-test was used before each analysis to check for homogeneity of variance.

Source of Variation	df	MS	Pseudo-F	P(perm)	perms	P(MC)	Denominator
Site (=S)	1	63.960	61.650	0.007	9164	0.0042	S × T(B)
Before vs. After (=B)	1	14.234	1.048	0.509	10	0.3797	T(B)
Time(B) (=T)	3	13.581	10.794	0.0003	9951	0.0001	Residual
S × B	1	2.466	2.377	0.222	9115	0.2181	S × T(B)
S × T(B)	3	1.037	0.824	0.495	9954	0.4909	Residual
Residual	20	1.258
Cochran test	C = 0.373		P = 0.149

). On the contrary, significant differences in the S × B interaction term were detected for the evenness (J’) of benthic assemblages (

Table 7. PERMANOVA on Pielou’s evenness index (J’) of benthic assemblages of Treatment and Control sites at Before (T1, T2) and three After (T3–T5) the bioremediating action. Cochran’s C-test was used before each analysis to check for homogeneity of variance.

Source of Variation	df	MS	Pseudo-F	P(perm)	perms	P(MC)	Denominator
Site (=S)	1	0.138	25.526	0.0215	9081	0.0155	S × T(B)
Before vs. After (=B)	1	0.097	30.820	0.0999	10	0.0125	T(B)
Time(B) (=T)	3	0.003	0.296	0.8410	9966	0.8324	Residual
S × B	1	0.128	23.699	0.0181	9115	0.0177	S × T(B)
S × T(B)	3	0.005	0.507	0.6860	9956	0.6771	Residual
Residual	20	0.011
Cochran test	C = 0.373		P = 0.149

).

Pairwise comparisons showed that the evenness of bioremediated assemblages changed after the implementation of the IMTA system, becoming more similar to that of the Control site (

Table 8. Pairwise PERMANOVA on Pielou’s evenness index (J’) of benthic assemblages of Treatment and Control sites at Before (T1, T2) and three After (T3–T5) the bioremediating action.

Pairwise of S × T(B) for pairs of levels of factor Site (Control vs. Treatment)
		T	P(perm)	perms	P(MC)
Before	C vs. Tr	48.2690	0.2616	3	0.0139
After	C vs. Tr	0.1190	1	38	0.9186
Pairwise of S × T(B)for pairs of levels of factor Time (Before vs. After)
Site	B	T	P(perm)	perms	P(MC)
Treatment	Bef vs. Aft	6.7528	0.0971	10	0.0062
Control	Bef vs. Aft	0.5353	0.6020	10	0.6316

).

The similarity percentages analysis (SIMPER) showed the overall average dissimilarity between the two sites ranging from 89.25, before, and 83.95, after the conversion, and highlighted the species cumulatively contributing up to 50% of the total dissimilarity (

Table 9. List of discriminating species contributing more than 50% of the cumulative dissimilarity by SIMPER analysis, between the Treatment (A) and the Control (B) sites, before (a) and after (b) the conversion of the plant.

(a)		A (before)	B (before)
EG	Species	Av.Abund	Av.Abund	Av.Diss	Diss/SD	Contrib%	Cum.%
V	Capitella capitata	90.33	3.67	27.93	1.34	31.29	31.29
I	Amphipholis squamata	0.50	13	6.32	1.05	7.09	38.38
IV	Cirrophorus nikebianchii	8.17	5.17	3.69	1.10	4.14	42.52
II	Amphiura chiajei	0.17	6.67	3.27	0.71	3.67	46.18
IV	Cirriformia tentaculata	0.17	7.17	3.02	0.47	3.38	49.57
(b)		A (after)	B (after)
EG	Species	Av.Abund	Av.Abund	Av.Diss	Diss/SD	Contrib%	Cum.%
V	Capitella minima	14.33	2.56	5.88	1	7.01	7.01
IV	Cirrophorus nikebianchii	0.22	9.44	4.87	0.93	5.80	12.81
II	Prionospio depauperata	3	11.56	4.63	0.99	5.52	18.33
II	Aricidea (Acmira) catherinae	1.67	8.89	4.58	0.69	5.46	23.79
II	Pseudolirius kroyeri	3.67	5.89	3.91	0.87	4.66	28.45
III	Abra alba	6.67	6.11	3.74	1.09	4.46	32.90
IV	Varicorbula gibba	6	3	3.47	0.89	4.14	37.04
I	Pisidia bluteli	2.22	6.56	3.20	0.57	3.81	40.85
IV	Caulleriella alata	1.22	5.78	2.34	0.77	2.79	43.64
I	Amphipholis squamata	4.44	1.89	2.01	0.66	2.40	46.04
II	Amphiura chiajei	0.89	3	1.97	0.68	2.34	48.38

). At the beginning, C. capitata was responsible for 31.29% of the dissimilarity, being very abundant in the Treatment site compared to the Control one. Only six species reached 50% of dissimilarity between the two sites before the IMTA conversion, while the number of species needed to reach 50% of dissimilarity doubled after the conversion.

4. Discussion

Aquaculture, especially mariculture, represents today the fastest-growing food production activity [1]. Notwithstanding the benefits that marine aquaculture has brought to society, the negative impact that fish farming can have on the environment is well documented [3,35]. By contrast, if carefully designed, mariculture may be an effective strategy to achieve positive ecological, economic, and societal impacts [38,40,68,69].

The IMTA can be a valuable tool for building a sustainable aquaculture industry. It is a very flexible concept that gives additional value to plants, relying on the appropriate choice of organisms with complementary functions within the ecosystem and with relevant economic value. The benefits of an integrated approach also include the recycling of waste by producing valuable biomass as a by-product, with consequent minimal impacts on the surroundings, and the opportunity to effectively farm fish in a healthy environment [22,23].

To evaluate the efficacy of the reared species as bioremediators in improving the environmental quality, particular attention was devoted to the changes in the structure of benthic soft bottom assemblages. During the ex-ante analysis conducted in 2018, the assemblage of the treatment site was dominated by the presence of the species C. capitata, which is generally abundant in organic enriched, low hydrodynamic environments, emphasizing the high organic content present in this area compared to the other sites, which by contrast showed a good health status [41].

The ecological changes that occurred in the treatment site after the conversion of the plant into an IMTA system were emphasized by the M-AMBI, which suggested an improvement in environmental conditions. Indeed, results obtained from the comparison of macrozoobenthic assemblages in terms of species composition, before and after the conversion of the plant, showed that the site where the bioremediating system was placed improved its ecological quality status after only one year, maintaining this condition through time, according to the EQR classification. Many of the species previously found in the treatment site were no longer collected after conversion, while others such as the polychaetes S. costarum, Aricidea (Acmira) catherinae, and Naineris setosa, belonging to lower EG (less tolerant to organic pollution), appeared to increase their abundance.

After the conversion, the Treatment site changed from “Moderate” rating in 2018 to “Good” in July 2020, reaching, during February 2020, the EQR classification of “High”, possibly due to the milder impact of the plant during winter months (the amount of feed spilled is less in winter than in the summer). Overall, M-AMBI showed an improvement in environmental conditions in the Treatment site, while the Control site persisted in a high-quality state with no ecological state variations through time.

This result was corroborated by a previous paper, conducted on the hard bottom under the cages, revealing an increase in biodiversity in time in the Treatment site [38], as well as an amelioration of some physico-chemical and microbiological variables recorded after the conversion of the plant [70].

The univariate (ANOVA) and multivariate (PERMANOVA) analyses revealed, however, a more complex scenario of spatial and temporal variability of benthic assemblages. Even before the implementation of the IMTA system, assemblages at the Control and at the putatively bioremediated sites did not significantly differ in February, while the difference between them was detected in July. Hence, PERMANOVA revealed that the change in environmental quality highlighted by the M-AMBI index (“Moderate” vs. “High”) derives from a more complex dynamic of change in assemblages.

The structure of the assemblages of the two sites was consistently different at samplings in July (T3 and T5) after the conversion of the IMTA system, but not in February 2020 (T4). Differences in the assemblages in the two sampling periods can be related to the seasonal variations in the physico-chemical properties of the water column, which appear to regulate the benthic composition directly or indirectly by influencing food availability, bottom-water oxygenation, and larval dispersion [71,72,73]. Moreover, the amount of feed released in the breeding cages is significantly lower in the winter period, suggesting greater similarity between the sites in February, as the sediment at the Treatment station is not affected by a large accumulation of organic matter as in the summer months.

While the M-AMBI index consistently characterized the Control as a “high quality” site from an ecological point of view, PERMANOVA emphasized the natural variability of benthic assemblages across time in the control site, that is one of the potential issues occurring when the BACI design with one control vs. one impacted location is used instead of a beyond-BACI design including multiple controls [46]. Specifically, higher temporal variability in the structure of benthic assemblages characterized the control site after the conversion of the IMTA system compared to the treatment site, possibly indicating a side effect of the IMTA system, which reduces the natural variability of the system in addition to mitigating the effect of aquaculture cages.

This result is in contrast with the M-AMBI index, which outlined an opposite pattern of temporal variability of environmental quality, describing the Control site as persistently in a “high quality state”, and the Treatment site ranging from “good” to “high quality” at different sampling times after the implementation of the IMTA system. This could be explained by the AMBI method, weighing species based on their ecology, based on the percentage of species belonging to five different EGs.

Probably univariate and multivariate analyses of variance were not able to highlight this change due to the great sensitivity to variations in the species composition and abundance of the assemblage. Benthic assemblages differing in species composition can share similar percentages of species for each EG but the identity of the species of a given EG may change time by time without affecting the EQR. In the case of the treatment site, the decrease in the number of species belonging to EG V (i.e., C. capitata) and the appearance of species belonging to EG I (i.e., S. costarum) have led to changes in the EQR, without increasing the assemblage variability across time as much as PERMANOVA detects for the Control site. Here, temporal changes in the assemblage (i.e., species identity and abundance), albeit being higher than in the Treatment sites, did not affect the relative proportion of EGs, thus resulting in a stable EQR classification among times.

Although changes in the structure of benthic assemblages were not a clear-cut effect of the bioremediating action and rather highlighted natural patterns of temporal variability, among univariate indexes, evenness was particularly diagnostic in revealing differences between control and treatment sites before and after the implementation of the IMTA system. Evenness of the assemblage increased after the implementation of the IMTA system and reached values comparable to those observed at the control site, which did not vary before or after the implementation of the IMTA system.

Species richness and evenness respond in different ways to human impact [74,75] and evenness often responds more rapidly to altered environmental constraints than species richness [76]. Therefore, although the overall richness or diversity of the system has not significantly changed due to the bioremediation effect, the distribution in the relative abundances of the species has changed. Since the dominant species in the treatment site, before the implementation of the IMTA system, were mainly opportunistic species tolerant to organic pollution, an increase in the evenness mirrors the improvement in environmental conditions, as suggested by the M-AMBI output, which relies on the proportion of disturbance-sensitive taxa.

Indeed, SIMPER analysis showed that few species were responsible for the dissimilarity (89.25) between the two sites before the IMTA implementation, with C. capitata contributing 31.29% of the total. High dissimilarity between the two sites (83.95) was also found after the implementation of the plant, but no species contributed to more than 10% to the dissimilarity. The benthic assemblage present in the treatment site has changed from the dominance of one species (belonging to EG V) to a more diversified one, reducing the dissimilarity with the control site. This result is consistent with M-AMBI outputs, which showed a progressive transformation of the treatment site assemblage, until it became more similar to the control one, in terms of species composition.

5. Conclusions

The bioremediation effect, derived by the conversion of a mariculture plant into an IMTA system was assessed following a BACI design. Changes in the structure of benthic assemblages, before and after plant conversion, were analyzed through M-AMBI and compared with the results of univariate and multivariate analyses of variance, to evaluate the sensitivity of the two methodologies. The following conclusions can be highlighted:

According to the M-AMBI index, the IMTA system proposed by the Remedia Life project, appears to be able to reduce the impacts related to the accumulation of organic matter in the sediment due to mariculture activities. Indeed, the treatment site improved its EQR value after the implementation of the IMTA system, while the control site maintained its starting EQR value over time.

PERMANOVA demonstrated a higher sensitivity to temporal variability in the structure of the benthic assemblages compared to the M-AMBI index, suggesting that it should be more widely applied for the detection of bioremediation effects, in addition to that of anthropogenic impacts. Yet, the selection of a unique control site as opposed to a set of reference conditions, and the replication of samplings at regular as opposed to randomly chosen time intervals might have partially prevented the detection of bioremediation effects, due to an underestimation of the magnitude of bioremediation. In contrast, a Beyond-BACI design would have allowed the detection of the causal relationship between human activities and environmental changes.

Therefore, additional samples in time and space are needed to better investigate the changes occurring in the assemblages from a statistical point of view. Further studies conducted on other plants with different environmental conditions are also needed to understand whether the bioremediation pattern can be generalized to a broad spectrum of environmental conditions.