Distribution Patterns of Benthic Foraminifera in Fish Farming Areas (Corsica, France): Implications for the Implementation of Biotic Indices in Biomonitoring Studies

Abstract

Corsican marine aquaculture is one of the highest contributors of fish production in France, which may result in environmental perturbations caused by organic matter (OM) accumulation under fish farms and impacting natural communities. This study aimed to (1) characterise the environmental conditions at two different fish farms, (2) monitor the response of benthic foraminiferal species to this activity, and (3) assess the accuracy of existing foraminiferal biotic indices. In 2017, sea floor sediment was sampled in transects from two Corsican fish farms for living foraminiferal and sedimentary analyses. Four indices were calculated and compared: exp(H′_bc), Foram-AMBI, Foram Stress Index and TSI-Med. A significant increase in total organic carbon (TOC) has been shown, mainly below the fish cages. Communities were characterized by a shift from high density, opportunistic and tolerant species under the cages to lower densities and more sensitive species further away. According to their distribution patterns along the TOC gradient, we propose to update the ecological group classification of seven species to improve Foram-AMBI’s accuracy and sensitivity: Triloculina oblonga and Quinqueloculina lamarckiana to Ecological Group (EG) I; Rosalina bradyi to EGIII; and Bolivina dilatata, Bulimina aculeata and Quinqueloculina stalkeri to EGIV. We recommend prioritising the use of TSI-Med and Foram-AMBI with the updated list to assess ecological quality in coastal waters of the Mediterranean Sea.

1. Introduction

With an ever-increasing global demand for fish, aquaculture has been developed to counterbalance and to limit the impact of fisheries on marine biodiversity and biomass. Nowadays, aquaculture provides about 42.2% of worldwide fish production [1,2]. In the Mediterranean Sea, the sea bass (Dicentrarchus labrax) and the gilt-head sea bream (Sparus aurata) are among the most harvested fishes (Greece, Spain or France) [3]. In France, Corsica has had among the highest production for these two species for many years [4]. The impact of fish farms on the surrounding environment is now widely acknowledged, mainly due to the accumulation of organic matter (OM) under the cages; up to 30% of the organic matter in the sediment below the cages being non-consumed food [5]. In addition to food-related releases, feces of fishes also accumulate, corresponding to up to 12% of the OM in the sediment [5]. High nutrients and OM loads from fish farming may induce important changes in the environmental conditions of the habitat and thereby alter benthic communities [6,7,8,9]. Considering that the Mediterranean Sea is an oligotrophic sea, it artificially creates more eutrophic zones in areas adjacent to the fish cages, stimulating benthic macrofauna [10] and foraminiferal abundances [11]. Conversely, under the cages, decrease in the seagrass Posidonia oceanica meadows have been reported several times [12,13,14], up to 11 ha in 10 years [15]. Furthermore, OM enrichment decreases the number of macro-invertebrate species present under aquaculture cages [12]. Ecological assessment of benthic habitat under fish cages in the Mediterranean Sea concluded that benthic communities were heavily disturbed [15,16,17]. Specifically, there is an increase in macro-invertebrate opportunistic species near fish cages and sensitive species mostly occur away from the cages [15,18,19].

Since 1999, the European Water Framework Directive (WFD) and then the Marine Strategic Framework Directive (MSFD) emphasised the need to monitor and assess the health of coastal, estuarine and continental waters. According to these legislations, the evaluation of water masses’ quality involves determining the Ecological Quality Status (EcoQS), noticeably in human-modified habitats like fish farming areas. In this context, several biological quality elements have been chosen, for instance, seagrasses communities [20], macro-algae [21] and benthic macro-invertebrates [22].

Recently, a group of protists belonging to meiofauna, benthic foraminifera, are now used as indicators of pollution (e.g., [23,24,25,26,27,28,29,30,31]). Their natural features, i.e., fossilisable calcareous test [32,33,34] and great abundances in limited volumes of sediment, make them particularly interesting candidates for environmental monitoring [35]. They also have a short and fast life cycle that allows a clear observation of environmental modification [35]. Under fish farms, several studies have reported a change in foraminiferal communities [11,36,37,38,39,40]. Noticeably, tolerant and opportunistic species abundances increase close to the oyster and fish farms while sensitive ones tend to disappear [36,41,42,43]. To increase their use in environmental monitoring surveys, biotic indices based on benthic foraminifera indices have been developed (see review [44]). They are based either on their diversity i.e., the exp(H′_bc) [45] or on the species-specific response to total organic carbon content i.e., Foram-AMBI [41,46,47], TSI-Med [48] and Foram Stress Index [49]. These indices were successfully applied to assess the health of benthic habitats in the Mediterranean Sea [31,37,41,48,49,50,51,52,53,54]. However, benthic foraminifera are not yet acknowledged as a biological quality element for the implementation of the WFD and MSFD.

In this context, the main objective of this study was to understand the distribution pattern of foraminiferal assemblages under fish farming cages in Corsica, and to further monitor the health of the benthic habitats underneath. Specifically, this study aimed to (1) characterize the environmental conditions under fish farming cages, (2) monitor the response of benthic foraminiferal species to the presence of fish aquaculture, and (3) assess the accuracy of foraminiferal biotic indices to reflect the impact of fish farms.

2. Materials and Methods

2.1. Sampling Strategy

Two fish farms were sampled in August 2017 by scuba divers from the STARESO (Station de Recherches Sous-marines et Océanographiques based in Calvi). Both farms are located along the west coast of Corsica (Figure 1) at 28 m (from 27 to 29 m) and 20 m (from 18 to 22.5 m) water depth for site B and D, respectively. Both fish farms have been awarded the French Label Rouge quality mark, indicating that they follow environmentally sound practices. They produce sea bass (Dicentrarchus labrax), gilt-head sea bream (Sparus aurata) and meagre (Argyrosomus regius, only site B). Fish farm B is located at the entrance of a bay, which is protected by islands from western and north-western winds, storms and it further allows for good seawater renewal (Figure 1). Outward seawater currents are oriented towards the west. Fish farm D is located in the inner part of a bay. It offers natural protection from western and south-western winds. The main seawater currents follow the coast and are oriented west- and eastward.

Five stations were sampled along a transect at each farm: under the cages (0) and then at 50, 100, 200 and 300 m away from the farm following a presumed gradient of organic matter enrichment at the sea floor (Figure 1). These aquacultures are the most active fish farms in Corsica, producing 850 and 150 tons of fish per year, respectively, for B and D.

Sampling was performed by pushing a core liner into the sediment (diameter: 7.4 cm) in patches of soft sediment between Posidonia oceanica meadows (Figure S1). Samples for benthic foraminifera analyses (sediment layer: 0–1 cm, n = 3) were preserved in a mix of 96% ethanol with rose Bengal dye (2 g L⁻¹) to distinguish living from dead specimens [34]. At each station, three replicate cores (sediment layer: 0–1 cm, n = 3) were sampled for organic matter analyses (e.g., total organic carbon, carbon-to-nitrogen ratio) and one additional core (sediment layer: 0–1 cm, n = 1) was sampled for grain size analyses.

2.2. Grain Size Analyses

A sample dedicated to grain size analysis was taken from the surface sediment (0–1 cm). For each station, the bulk sediment was sieved over a 2000 μm mesh size. Then, the grain size distribution was measured in three aliquots with a Malvern Mastersizer 3000 (Malvern Inc., Malvern, UK). When the three aliquots had a similar grain size distribution, a single measurement was chosen randomly to describe the grain size distribution at each station.

2.3. Organic Matter Analysis

Sediment samples for organic matter analyses were preserved at −20 °C at the laboratory. Total organic carbon and nitrogen contents were obtained after sediment grinding by hand, by quantifying the total amount of carbon and nitrogen (measured in dried samples). The amount of inorganic carbon and nitrogen (measured in samples burned at 550 °C for 5 h) was subtracted. It was determined with an elemental analyzer (Thermofisher Flash 2000, Laboratory of Oceanology and Geosciences in Wimereux, France) and expressed as the % of total organic carbon (TOC) and total organic nitrogen per total weight of dry sediment. The carbon-to-nitrogen (C:N) ratio was calculated at each station to determine the terrestrial or marine origin of the organic matter.

2.4. Foraminiferal Analysis

Samples were washed with tap water through a succession of 500, 125 and 63 µm mesh sieves to remove clay, silt and any excess of dye. Each residual fraction was dried at 50 °C for 12 h. Quantitative analysis of benthic foraminifera was performed on the 125–500 μm fraction. Density separations were performed for all samples with sodium polytungstate (SPT) solution of 2.3 density with 161 g of SPT powder in 60 mL of Milli-Q^® water [55]. For all samples except for B0-3, B50-2, D0-2 and D0-3, all living specimens were dry picked. For these 4 samples showing high densities, 300 stained specimens per sample were dry picked [45]; the total number of individuals being then estimated using the precise weight of sediment picked compared to the total weight of sediment in the sample. Only specimens with dense and brightly red/pink stained protoplasm were considered as alive [35] and identified until species level, when possible, following taxonomical books on Mediterranean benthic foraminifera (e.g., [56,57]). Species names follow the World Register of Marine Species [58].

2.5. Faunal Analyses

We calculated specific richness, absolute densities (indiv.50 cm⁻²) and relative abundances per station (%) after pooling the three replicate samples together. Then, four biotic indices were calculated based on benthic foraminifera: the diversity index exp(H′_bc) [45], and three sensitivity-based indices, i.e., Foram-AMBI [41,46,47], FSI [49] and TSI-Med (Tolerant Species Index-Mediterranean, [48]).

As the Shannon index H′ is biased when there are unobserved species in the community, a common problem with under-sampling [59], we used the bias-corrected version of Shannon’s index (H′_bc), which has little bias [60]. Shannon’s index is an entropy rather than a diversity index. The entropy gives the average uncertainty of the identity of an individual picked from the community, not the number of species in the community [61,62]. It is calculated as follows (Equation (8) in [59]):

$${{\mathrm{H}}^{\prime }}_{\mathrm{bc}}=-{\displaystyle \sum }_{i=1}^S\frac{{\tilde{\pi }}_i\log \left({\tilde{\pi }}_i\right)}{1-{\left(1-{\tilde{\pi }}_i\right)}^n}I\left(A_i\right)$$

where ${\tilde{\pi }}_i=\left(\raisebox{1ex}{$X_i$}\!\left/ \!\raisebox{-1ex}{$n$}\right.\right)\widehat{C}$, $\widehat{C}=1-\left(\raisebox{1ex}{$f_1$}\!\left/ \!\raisebox{-1ex}{$n$}\right.\right)$, A_i denotes the event that the ith unit is included in the sample and I(A_i) is the usual indicator function (i.e., I(A_i) = 1 when the event A_i is true and I(A_i) = 0 otherwise), S is the number of species in a community (labeled from 1 to S), $X_i$ is the number of individuals of the ith species observed in the sample, n is the number of individuals, $f_1$ is the number of species with 1 individual in the sample [59].

It can be converted to true diversity, the effective number of species, with the exponential function N1 = exp(H′_bc) [63]. Exp(H′_bc) gives the number of species that would, if each were equally common, produce the same H′_bc as the sample. All data represent the pooled counts from three replicates per station (rather than the average).

For Foram-AMBI, species from open waters in the Mediterranean Sea were assigned to five ecological groups according to their response to TOC content [47]. Ecological Group I (EG-I) is composed of taxa sensitive to small increases in TOC; EG-II of taxa indifferent to TOC content; EG-III of third-order opportunists to excesses of TOC; EG-IV of second-order opportunistic species present at high TOC levels; and finally, EG-V of first-order opportunistic species able to resist strong disturbances and excesses in TOC. The following formula is used to compute Foram-AMBI:

$$\mathrm{FForam}-\mathrm{AMBI}=\frac{\left(0\times \%\mathrm{EG} \mathrm{I}\right)+\left(1.5\times \%\mathrm{EG} \mathrm{II}\right)+\left(3\times \%\mathrm{EG} \mathrm{III}\right)+\left(4.5\times \%\mathrm{EG} \mathrm{IV}\right)+\left(6\times \%\mathrm{EG} \mathrm{V}\right)}{100}$$

Following the procedure described in [22], species not assigned to any of the five ecological groups (see Table S1) were not considered in the calculation of the relative proportions of the five EGs. Consequently, the cumulative proportion of EGI to EGV is equal to 1. Foram-AMBI values range between 0 (bad ecological quality) and 6 (high ecological quality). In this study, we first used the species classification proposed in 2018 [47] for open waters in the Mediterranean Sea. In a second step, we proposed modifying the species assignment of six of the major species identified in our study area. Three of these species (Triloculina oblonga, Quinqueloculina lamarckiana, Q. stalkeri) were not assigned yet by Jorissen et al. [47], and exhibited a clear distribution pattern along the TOC gradient. The three other species (Rosalina bradyi, Bolivina dilatata and Bulimina aculeata), originally classified by Jorissen et al. [47], exhibited a response to the TOC gradient different from their original assignments in our study area. According to these study results and to their known ecology, we assigned or re-assigned these species to EGs. This is further explained in the discussion. We calculated therefore the “Foram-AMBI” according to the original list [47] and the “Foram-AMBI-updated” according to the proposed updated list based on the present study findings.

The FSI considers only two groups of indicative species: sensitive and tolerant species [49]. This classification is based on literature review from their study area in the Eastern Mediterranean and we decided to follow this original classification to apply FSI to our dataset. FSI is defined by the following equation:

$$\mathrm{FSI}=\left(10\times \%\mathrm{Sensitive}\right)+\left(\%\mathrm{Tolerant}\right)$$

As for the Foram-AMBI, some species were not assigned to any group (see Table S1) however they were not excluded from the dataset since no specifications were given for this aspect [49]. According to these authors, all taxa should be classified either as “Sensitive” or as “Tolerant”, but this was, in our opinion, not possible to decide for all species in our study area. FSI values range between 0 (high ecological quality) and 10 (bad ecological quality).

For TSI-Med, the index value is based on the cumulative proportion of a group of stress-tolerant taxa and is corrected for the expected proportion of this group under unpolluted reference conditions [48]. The latter value is estimated on the basis of a grain size parameter (% grains < 63 μm). We used the most recent species list provided in 2021 [31] which has been defined based on literature review, to identify tolerant species to increased OM and potential oxygen depletion. TSI-Med is calculated using the following equation:

$$\mathrm{TSI}-\mathrm{Med}=\frac{\%TSx-\%TSref}{100-\%TSref}\times 100$$

where

$$\%TSref=0.1547\times \%63 \mathsf{\mu }\mathrm{m}-0.3638$$

and %TSx is the percentage of tolerant species found in the sample. TSI-Med values range between 0 (high ecological quality) and 100 (bad ecological quality).

As the present study aims to evaluate the accuracy of the different biotic indices tested, we decided to not use the ecological quality status but rather the quantitative data given by these indices. Indeed, the transformation of these data into 5 EcoQS will strongly depend on the boundary limits proposed. This issue, i.e., how to set the boundary limits, is another scientific question that will need to be addressed separately.

2.6. Numerical Analysis

Kruskal-Wallis tests were performed on the TOC data as no normality and homogeneity has been found by Shapiro and Bartlett tests. Dunn’s post-hoc tests were computed when a significant difference was detected by the Kruskal-Wallis tests.

Correlation between indices and environmental parameters was investigated using Kendall’s coefficient of rank correlation (τ). Kendall’s coefficient of correlation was used in preference to Spearman’s coefficient of correlation (ρ), because Spearman’s ρ gives greater weight to pairs of ranks that are further apart, while Kendall’s τ weights each disagreement in rank equally [64].

For the species community data, initial detrended correspondence analyses indicated that the gradient is short (2.2SD) in the foraminiferal data. Consequently, principal correspondence analysis (PCA) was applied to the foraminiferal data. Replicate samples were pooled, and relative abundances data were log(x + 1) transformed prior to analysis. A heatmap was further used to visualize the distribution of species between stations (relative abundance > 5%, untransformed pooled data).

Statistical tests were performed with the statistical language R version 4.0.3 [65] and statistical packages (entropy v1.2.1, vegan v2.5-7, irr v0.84.1, gplots v3.1.1 and ggplot2 3.3.5).

3. Results

3.1. Environmental Parameters Analayses

Both sites are characterized by sandy sediment (sand fraction (>63 µm) = 92 ± 7%, all stations combined, Figure 2). Grain size analyses showed that fine fraction (<63 µm) percentages ranged from 0.9 (station B300) up to 23.9% (station D0, Figure 2). At site D, the fine fraction (<63 µm) decreased from the station below the cage (D0) towards more distant stations (Figure 2). Site B showed an increase in the fine fraction until 200 m from the fish farm with a lower percentage of this fraction at 300 m. A similar trend was observed for the fraction 500–2000 µm, with an increase along the transect (Figure 2).

Sites B (2.02 ± 1.04%) and D (1.71 ± 1.28%) show relatively similar total organic carbon percentages with intra-site differences (Figure 3). Indeed, at site B, TOC values at station 0 (i.e., under the fish cages) was significantly different from stations B100 (Dunn test, p < 0.05) and B200 (Dunn test, p < 0.001) with higher TOC content at station B0. For site D, station D0 had a significantly higher TOC content than stations D100, D200 and D300 (respectively Dunn test, p < 0.05, p < 0.05 and p < 0.001, Figure 3). There was also a statistically significant difference between station D50 and D100 (Dunn test, p < 0.05, Figure 3).

Conversely to TOC, the two sites exhibited different trends regarding organic C:N ratios (Figure 3). This ratio was rather stable at site B where no significant differences were found. On the other hand, C:N ratios increased slightly with a greater distance from the fish farm at site D with statistically higher value at site D300 (p < 0.05 at D300, and p < 0.001 at D50, Figure 3).

3.2. Foraminiferal Analyses

3.2.1. Abundance and Diversity

Foraminiferal analysis of the 125–500 µm fraction revealed a heterogeneity between samples in terms of density of foraminifera for 50 cm² (ind.50 cm⁻²). Densities varied from 126 ind.50 cm⁻² (station B300) to a maximum of 732 ind.50 cm⁻² (station D0; Figure 4, Table S1). Stations under the fish farms (i.e., B0, B50, and D0) exhibited higher densities than stations further away from the fish farms (Figure 4).

We identified 102 different species among which 93 were identified to the species level (Table S1). A variation of species richness was observed among all of the stations with a maximum richness of 75 at station B50 and a minimum at station D300 with 37 species (Figure 4). However, species richness was clearly influenced by the total abundances at the stations (see Figures S2 and S3 for Kendall’s correlation). The exp(H′_bc) exhibited a pattern that is similar between site B and D with higher biodiversity observed at intermediate stations of the transects (i.e., B100 and D100).

3.2.2. Faunal Composition

According to the PCA performed on relative densities of the total fauna, site B and D were well separated on axis 1 of the PCA representing 24.9% of the total variance of the fauna (Figure 5), meaning that both sites are characterized by different assemblages, probably because of their different geographical location. For site D, stations D100, D200 and D300, further away from the cages, have very similar faunal composition (Figure 5). For site B, the faunal composition seems to be quite different between all of the stations.

The heat map based on the 20 major species, present in more than 5% of the total fauna in at least one of the stations from sites B and D, showed different distribution patterns according to the species considered (Figure 6). For example, Bolivina dilatata and Bulimina aculeata exhibited a similar pattern at both sites with decreasing relative densities with increasing distance from the fish cages. On the other hand, species such as Adelosina cliarensis (site D), Astrononion stelligerum (site B), Lagenammina atlantica (site B, and D to a lesser extent), Pyrgo oblonga (site D), Triloculina oblonga (and T. tricarinata and T. trigonula to a lesser extent in site D) showed opposite trends with decreasing densities towards the cages. Eggerelloides scaber had a maximum of density at 50 m away from the cages at both sites. For Quinqueloculina spp., mainly present in site D, all species seemed to behave differently with a maximum density for Q. stalkeri under the fish cage whereas all other species (Q. bosciana, Q. lamarckiana and Q. seminula) showed increasing abundances further away. At site B, the epiphytic species Lobatula lobatula and Rosalina bradyi dominated the fauna under the cage.

In detail, Quinqueloculina lamarckiana and Triloculina oblonga (Figure 7a,b) exhibited maximum relative abundances when the TOC content was low with a distribution pattern similar to species from EGI (Figure 7g). Rosalina bradyi (Figure 7c) showed a maximum around 3% TOC; resembling to distribution pattern of species from EGIII (Figure 7g). Bulimina aculeata, Bolivina dilatata and Q. stalkeri (Figure 7d–f) had exponential relative abundances along the increasing gradient of TOC, a behaviour similar to species from EGIV (Figure 7g). These observations provide evidences for the assignment/re-assignment of these species. This is further addressed in the Discussion Section.

3.2.3. Proportions of Ecological Groups

At site D, the relative densities of tolerant to opportunistic species increased towards the fish cages whatever the index considered (Figure 8). The group of sensitive species identified in Foram-AMBI and FSI showed a slight increase with distance (Figure 8). For site B, there was no apparent trend in the different EGs from Foram-AMBI with distance from the cage. Actually, EGIII (3rd order opportunist species) dominated at the intermediate station B50. For TSI-Med and FSI, only station B0 below the cage exhibits higher proportions of tolerant species compared to all the other stations. Surprisingly, the proportion of sensitive species is also higher at B0 for FSI.

It is important to highlight that some species identified in this study were not assigned to any group for Foram-AMBI and FSI (“Not assigned species” (NA), Table S1). These proportions varied from 14.9% (D0) to 30.7% (D200) for Foram-AMBI, with 7 of the 10 stations studied with %NA > 20%. For FSI, it varied between 4% (B0) and 45% (B300), with 4 out of 10 stations with %NA > 20% (Table S1).

3.2.4. Biotic Indices and Environment Parameters

At site B, exp(H′_bc) and Foram-AMBI did not exhibit a clear trend with respect to TOC or distance from the fish cages. The highest values were observed at station B50 for Foram-AMBI and at B100 for exp(H′_bc). FSI and TSI-Med showed higher index values below the fish cages compared to more distant stations, suggesting opposite trends in terms of ecological quality evolution.

At site D, exp(H′_bc) had its maximum 100 m away from the fish farming cages, like at site B (Figure 9). Foram-AMBI values were quite stable along the transect, although the three stations closer to the aquaculture exhibited higher values. Regarding FSI and TSI-Med, the former had increasing values while the later had decreasing values (Figure 9).

The index exp(H′_bc) did not correlate with distance and TOC (Figures S2 and S3). Foram-AMBI did not correlate either with distance or TOC at both sites (Figures S2 and S3). FSI was significantly positively correlated with TOC at site B (p < 0.05, Figure S2) (i.e., better quality closer to the fish cages), although this is the opposite trend to what would be expected, and with distance at site D (p < 0.05, Figure S3). TSI-Med was significantly correlated with TOC at site B (p < 0.01) but not at site D. FSI was significantly positively correlated with Foram-AMBI at site D (p < 0.01), with TSI-Med at site B (p < 0.05) and negatively at site D (p < 0.05).

Foram-AMBI computed with this study-updated list (Foram-AMBI-updated) had higher values (Figure 10) than Foram-AMBI based on the original list [47]. Moreover, the percentage of NA species largely decreased in Foram-AMBI-updated, varying from 9 to 16.3% in the study area. None of the stations exhibited %NA > 20%. Specifically, values were higher at stations B0 and B50 for Foram-AMBI-updated, there was a shift then at station B100, B200 and B300 (Figure 10). At site D, the highest value was recorded under the cages, and there was a significant decrease toward site 300 (Figure 10). Foram-AMBI-updated was further significantly correlated with TOC at site D (p < 0.05, Figure S3). Foram-AMBI-updated was significantly positively correlated with TSI-Med and negatively with FSI at site D (p < 0.01 and p < 0.05, respectively, Figure S3).

4. Discussion

4.1. A Limited Effect of Fish Farming on Habitat Features

In Corsica, sites B and D correspond to large aquacultures with high annual productions and consequently higher fish densities and food inputs compared to other fish farms on the island. However, note that producers are engaged in sustainable fish-farming practices, further limiting the harmful effects on the benthic habitat. Except in a few stations, content in the fine fraction remained low compared to several studies made on aquaculture impacts (e.g., [36,37,53]) or on TOC content impacts [31]. In fact, the studied stations were mostly composed of fine or coarse sand, a composition typical of high energy hydrodynamic regimes. However, it has been shown that in some cases, accumulation of TOC can be seen in sandy bottoms without an increase in silty fractions because of the removal of fine particles with the currents [12,37,66].

Noticeably higher TOC values were reported under the cages. Since the different sites are located within Posidonia meadow habitats, part of the TOC content could correspond to Posidonia leaves and roots. However, C:N ratio values under the cages at sites B and D, respectively 13.39 (±2.32) and 11.59 (±1.94), are closed to characteristic values from aquaculture-derived organic matter [5] and not to the ones induced by Posidonia oceanica, typically with a C:N ratio of about 27.4 [67], corresponding to refractory organic matter. This suggests that the enrichment is due to OM produced by aquaculture activities (i.e., fish feeding and feces). The organic matter accumulation remains largely constrained under the cages. This confirmed previous studies’ findings reporting a moderate and spatially-limited impact of aquaculture in terms of direct OM loads [5,8,36,37,68].

4.2. Diversity-Based Indices May Not Be Adapted to Coastal Waters in the Mediterranean Sea

The foraminiferal faunas respond clearly to the TOC enrichment observed in the stations close to the fish cages with increasing absolute densities and specific richness. These results suggest that TOC under fish farming cages may most likely be composed of labile OM on which benthic foraminifera are feeding rather than refractory OM derived from Posidonia meadows. Furthermore, the positive response of foraminiferal assemblages to a moderate increase in TOC content is in line with previous studies in the Mediterranean Sea [11,31,37,52,69]. Specifically, in the early stages of enrichment, TOC can stimulate the size and specific richness of benthic communities, thanks to an increase in the biotic capacity of the environment [70]. Furthermore, we observed the same trend with the diversity index exp(H′_bc) which is less sensitive to sampling effort than the specific richness; noticeably, maximum values were reported at intermediate stations (50 and 100 m). In mesotrophic to eutrophic environments, a TOC enrichment would result in a decrease in diversity. However, as above-stressed, in oligotrophic environments in the first stage of the TOC enrichment, diversity would be stimulated. Then diversity would decrease when the TOC content become too high with species tolerant to degraded conditions taking over sensitive ones. The non-monotonic trend observed here confirms observations from previous studies in the Mediterranean [31,48]. This succession follows the model of [10] the often named “the hypertrophic zone theory” in other studies [11,38,39]. Oligotrophic characteristics of the Mediterranean Sea habitats further enhance this phenomenon [69] and induce a strong dichotomy between stations under the influence of fish farming and stations not influenced. The present findings further suggest that the diversity index exp(H′_bc) may not be accurate to evaluate environmental quality in oligotrophic coastal ecosystems such as coastal waters in the Mediterranean Sea; conversely, this index seems to be better adapted to more sheltered environments in the Mediterranean Sea [50,54].

4.3. Species Response to Fish Farming Cages: Improving Foram-AMBI Original Species List Assignment to Ecological Group, the Foram-AMBI Updated

Our study area shows a clear pattern of TOC content according to the distance from the fish cages and is therefore a good opportunity to study the response of major species to this gradient and propose an update to the EG assignments originally published in [47]. This work would potentially help reducing the percentage of not-assigned species and improve the reliability and performance of Foram-AMBI. The heatmap of major species densities at both sites already highlighted that some species exhibited clear distribution patterns along the distance from the fish cages (Figure 6). To follow the same approach as performed in [47], the relative densities of seven major species were plotted to determine their response along the TOC gradient (Figure 7 and Figure S4). Then, the distribution patterns were compared with the theoretical model of succession for species of the five EG along the TOC gradient from [47] (Figure 7g) and their original assignments. Among the 20 major species, three were not assigned to any EG yet: Triloculina oblonga, Quinqueloculina lamarckiana and Q. stalkeri. For ten of these species, the distribution patterns observed in our study area were in agreement with the assignments proposed by Jorissen et al. [47]: Adelosina cliarensis (EGI), Eggerelloides scaber (EGIII), Neoconorbina terquemi (EGI), Q. stelligera (EGIII), T. trigonula (EGI), and to a lesser extend Gavelinopsis praegeri (EGI), Pyrgo oblonga (EGII), Q. bosciana (EQII), Textularia agglutinans (EGI) and T. tricarinata (EGI). For the other seven species, their responses to TOC in our study area were somewhat different from their initial assignments. However, we only proposed a new classification for three of them, Rosalina bradyi, Bolivina dilatata and Bulimina aculeata, as explained thereafter. For Astrononion stelligerum, Lagenammina atlantica, Lobatula lobatula and Q. seminulum, we considered that the arguments were not strong enough to modify the original assignments proposed by Jorissen et al. [47] (e.g., low relative abundances, absence of literature evidences to support our observations). We also proposed an assignment for the three originally not-assigned species (T. oblonga, Q. lamarkiana and Q. stalkeri).

Adelosina cliarensis, Gavelinopsis praegeri, Neoconorbina terquemi, Quinqueloculina lamarckiana, Triloculina oblonga, T. tricarinata and T. trigonula exhibited increasing relative abundances away from the fish farming cages, responding negatively to TOC (Figure 6). These species were assigned as sensitive species (EGI, [47]), except for T. oblonga (NA) and Q. lamarckiana (NA). In detail, Q. lamarckiana and T. oblonga (Figure 7a,b) exhibited maximum relative abundances when the TOC content was low with a distribution pattern typical of EGI (Figure 7g). We propose assigning the non-assigned species Q. larmarckiana and T. oblonga [47] to EGI; which is consistent to their classification as sensitive species in FSI [49].

Model of distribution pattern of Rosalina bradyi (Figure 7c) had a bell-shaped curve with a maximum around 3% TOC. Furthermore, high relative abundances were observed for Rosalina bradyi under the fish farming cages, particularly at site B. The observed distribution pattern for R. bradyi is in contradiction with its assignment to EGI by [47] and to the sensitive group by [49]. It responded positively to the first stage of TOC enrichment. Rosalina bradyi is usually recognised as a sensitive species because of its epiphytic behaviour. More precisely, this species is predominantly motile, living in the sediment, but is able to leave temporarily attached to leaves of Posidonia oceanica, with a short life cycle [71]. One hypothesis to explain the high relative abundances under the cages could therefore be that specimens were transported by the current from Posidonia oceanica meadows and accumulated under the fish farms. However, high abundances for R. bradyi were also reported in fish farms in Tunisia [37] and in abalone farms in Korea [72]. Furthermore, Rosalina spp. are often abundant in organic-rich habitats [73], being influenced by the quality of the OM [74]. This further suggests that R. bradyi may benefit from the organic matter released by aquaculture.

The distribution pattern of Rosalina bradyi along the TOC gradient may be compared with that of other species classified as EGIII by [47] and present as major species in our study area. This is the case for Quinqueloculina stelligera, Eggerelloides scaber and Textularia agglutinans. They were dominant at stations B50 (8%, 26.4% and 7.2%, respectively) and D50 (3.5%, 7.5% and 6.6%, respectively), but their abundances dropped under the fish farming cages. It is interesting to notice that E. scaber and T. agglutinans are also considered to be tolerant species in FSI, whereas TSI-Med did not consider these species. The behaviour of these two species along the TOC gradient observed in this study gives a good argument to add these species to the list of tolerant species proposed by [48]. Similarly, the miliolid species Q. stelligera can tolerate the first stage of TOC enrichment but then its relative abundance dropped when TOC values reached about 4%, which is a typical behaviour of third-order opportunistic species [47].

According to literature observations, the density distribution of Rosalina bradyi along the present study’s TOC gradient and the similarity of this distribution with other EGIII species and theoretical models of succession of species from [47], R. bradyi should be assigned to EGIII, in our opinion.

Bulimina aculeata exhibited a clear trend of increasing relative densities towards the fish farming cage (Figure 6). In the Mediterranean, this species was reported to positively respond to labile OM [75] and to increasing TOC content [49,76]. It is an opportunistic species characterising highly stressed environments [75,77,78] and is considered stress-tolerant in the FSI and TSI-Med species lists [48,49]. It was originally assigned to EGIII [47], however the distribution pattern of B. aculeata is different from the typical EGIII distribution in our study area (i.e., Eggerelloides scaber, Quinqueloculina stelligera). Specifically, its abundance increased when species from EGIII started to decrease, with maximum densities at the highest TOC values in our study area, which does not correspond to the extreme conditions of the TOC range presented in the theoretical model of EG distribution (Figure 7g). According to the present findings, B. aculeata should be re-assigned to EGIV in Jorissen et al.’s species lists.

Bolivina dilatata showed the same pattern as Bulimina aculeata, with increasing relative abundance towards the cages (Figure 6). Like other bolivinids, this species is able to flourish under harmful conditions [27] and is often considered a stress-tolerant or opportunistic species [41,48,49]. However, this species was assigned to EGII in open marine environments from the Mediterranean Sea [47]. According to the theoretical Foram-AMBI model of species succession along the TOC gradient, its presence in relatively high densities just below the fish cage and its distribution pattern along the TOC gradient suggests that these species should be re-assigned to EGIV like B. aculeata in the species list [47].

Quinqueloculina stalkeri is not yet included in the species list that we used originally. In the FSI and TSI-Med lists, it is respectively classified as a sensitive species like all miliolids [49] and not included in the stress-tolerant group [48]. Actually, this species exhibited the same pattern as Bolivina dilatata and Bulimina aculeata with low relative abundance at low TOC contents and a sudden increase in relative abundance up to 15% at site D under the cages when TOC reached 3.8%. Like B. dilatata and B. aculeata, this species bloomed when EGIII species’ relative abundances dropped. This suggests that Q. stalkeri is an opportunistic species, which, according to the theoretical model, can be assigned to EGIV.

4.4. Comparison of Sensitivity-Based Indices’ Efficiency in Monitoring the Impact of Fish Farming

Biotic indices based on groups of indicative species better reflect the impact of fish farms than diversity indices. However, we obtained contradictory results between indices. Noticeably, TSI-Med and Foram-AMBI-updated performed well compared to FSI and Foram-AMBI, since they correlated statistically with TOC and distance from the fish cages.

Higher values of FSI reported at stations with the highest TOC content at site B would indicate a better quality at these stations [49], when in fact, it should be the opposite; this shows that FSI did not perform well in highlighting the impact of the fish farming cages. One bias in the species classification proposed by [48] is that they considered all miliolids as sensitive species whereas [47] and this study identified that some miliolid species have tolerant or opportunistic behaviour, i.e., Quinqueloculina stelligera and Q. stalkeri. As FSI is giving more weight to the sensitive species (10 times more than to tolerant species), the assignment issue of the above-mentioned species impacted FSI more than the other indices.

As for Foram-AMBI, using the original list [47], we did not observe any clear trend. The high number of not-assigned species may explain this. Indeed, it is highly recommended not to calculate AMBI, and, hence, Foram-AMBI, when the percentage of NA species is higher than 20% [22]. In our study area, seven of the ten stations at site B and D exhibited more than 20% NA species, hampering the use of this index at the study sampling sites. Conversely to TSI-Med where only the tolerant species are assigned, the main challenge with the FSI and Foram-AMBI biotic indices is that all identified species are considered in the calculation of the index. Therefore, a huge effort has to be made to assign as many species as possible despite the difficulty due to the ecological and biological plasticity of the species [79,80,81]. Foram-AMBI values obtained in the present study may not be robust enough to give accurate information on the impact of the fish farms. For instance, Foram-AMBI does not detect the lower quality at B0 and D0 as does TSI-Med. It shows the urgent need to assign more species to ease the implementation of this index. The assignment of three new species to EGs (Triloculina oblonga, Quinqueloculina lamarckiana, Q. stalkeri) and the re-assignement of three others (Rosalina bradyi, Bolivina dilatata and Bulimina aculeata) allowed us to improve the accuracy of Foram-AMBI. Specifically, Foram-AMBI computed with this study-updated list (Foram-AMBI-updated) had higher values (Figure 10) than Foram-AMBI based on the original list [47]. This was mainly due to the attribution of EGIII to R. bradyi and the re-assigment of B. dilatata and B. aculeata to EGIV. Foram-AMBI-updated allowed us to better distinguish between the stations under the cages and further away from the cages; hence this index better reflected the deleterious effect of the fish farming cages on the benthic environment. It further suggests that more work is still required to improve the species’ assignment to EG, and to further facilitate the implementation of Foram-AMBI.

Higher values for Foram-AMBI-updated [41] and lower values for TSI-Med [48] reflected the environmental degradation in the fish farming area. Noticeably, according to TSI-Med and Foram-AMBI-updated, the impact of fish farming on the ecosystem is constrained to just below the fish cages at site B (or within 50 m according to Foram-AMBI) and within 100 m from the fish farm at site D. This corroborated quite well with the higher TOC values under the cages. In general, the effects of aquaculture are mostly limited in the sediment under the fish farms [11,16,36,82]. The abundances of opportunistic, tolerant and sensitive species were clearly constrained by the presence of fish cages; opportunistic and tolerant species were more abundant in the assemblages impacted by the farm and, conversely, sensitive species were more abundant outside the fish farm, confirming previous studies in fish farm areas [11,36,37,43,83]. Moreover, according to TSI-Med and Foram-AMBI-updated, which seems to work particularly well in this case study, the impact of the fish farm on the ecosystem quality would be slightly higher at site D than at site B, the former being more sheltered in a bay than the later which is at the entrance of a bay and which may contribute to the stronger impact observed at site D. Although the aquaculture located at site B produces roughly 5 times more tons of fishes than the aquaculture at site D, the fish cages are located at the outer part of the bay where the currents are strong enough to sweep away the OM discharges produced by the aquaculture activity, hydrodynamics playing a key role in limiting or enhancing the impact of aquaculture [5]. Conversely, aquaculture site D is located at the end of a bay where the hydrodynamic conditions are low therefore allowing the OM input from the fish farms to settle down and remain in the sediment.

Biotic indices based on species sensitivity were able to detect the more or less important decrease in the ecosystem quality below the fish cages. The TSI-Med and Foram-AMBI (using the updated list) highlighted quite well the ecological degradation of the benthic habitat under the fish farming cages. It confirms that foraminifera can be relevant for aquaculture impacts monitoring [36,37,43,84,85], and should be considered as a Biological Quality Element for the implementation of marine legislations.

5. Conclusions

This study reveals how benthic foraminiferal communities are affected by off-coast Corsican aquacultures and how biotic indices based on species sensitivity are able to detect and evaluate the ecological quality of the ecosystem. Generally, the impacts of the studied farms were restricted to a small area around the farms (50m), the TOC enrichment being moderated. Foraminiferal assemblages are characterised by opportunistic and tolerant species under the cages, while more sensitive species occurred further away from the farm. This study confirmed that diversity indices such as exp(H′_bc) may not be adapted to detect organic matter pollution in oligotrophic systems like the open marine environments of the Mediterranean Sea. TSI-Med and Foram-AMBI-updated are the most relevant and sensitive indices amongst the sensitivity-based indices. Regarding Foram-AMBI, the here-proposed update of the original list improves the accuracy of this index to detect the decrease of the health of the ecosystem below the fish cages. The modifications to the original list are suggested and summarized in

Table 1. Summary of ecological group changes from the originally used list [47] to this study’s updated list.

Species	Originally Used List [47]	This Study
Bolivina dilatata	EGII	EGIV
Bulimina aculeata	EGIII	EGIV
Quinqueloculina lamarckiana	NA	EGI
Quinqueloculina stalkeri	NA	EGIV
Rosalina bradyi	EGI	EGIII
Triloculina oblonga	NA	EG I

.

These results emphasize the need to conduct more quantitative studies to improve the knowledge of benthic foraminiferal responses to TOC in order to improve the species’ assignment to EG and further facilitate the implementation of Foram-AMBI. They also highlight the advantage of TSI-Med which requires only the identification of tolerant species. At this stage, we recommend prioritising the use of TSI-Med and Foram-AMBI with the updated list to assess ecological quality in open marine environments of the Mediterranean Sea until further improvements are made to FSI.