Response of Benthic Macrofauna to Eutrophication in a Mesocosm Experiment: Ecosystem Resilience Prevents Hypoxic Conditions

A benthic-pelagic mesocosm experiment was performed for the study of the response of the benthic macrofaunal community under a eutrophication gradient. Nine mesocosms were deployed in the facilities of the Hellenic Center for Marine Research in Crete, in the eastern Mediterranean. The mesocosms were 4 m deep, contained 1.5 m3 coastal water, and included 85 liters of undisturbed sediment at the bottom. The experimental design included two eutrophication levels (Low and High) for the 58-day duration. Macrofaunal samples were collected at the end of the experiment from each mesocosm and compared to the ones collected at the beginning of the experiment from the sediment collection area. Results showed significant differences of the High eutrophication treatment in terms of macrofaunal species composition, diversity, ecological status and ecosystem processes. The increased availability of organic matter in the sediment caused differences in macrofaunal community structure by favouring deposit-feeding species with high bioturbation ability, which significantly increased their abundance. The increased bioturbation potential of the new community combined with the high organic matter consumption contributed to the oxygenation of the sediment within the mesocosm, preventing the creation of hypoxic conditions in the sediment and maintaining ecosystem health despite the highly eutrophic conditions and significant changes in sediment geochemical variables.


INTRODUCTION
The effects of organic enrichment on abundance, diversity, and biomass of marine macrobenthic communities were initially described in detail by Pearson and Rosenberg (1978), and have been proven in many publications since then. Based on that model, numerous methodologies have been proposed using benthic macrofauna as an indicator for assessing ecosystem health. In the context of the European Water Framework Directive (WFD), a series of benthic indicators were developed and have been widely used, including the AZTI Marine Biotic Index (AMBI) (Borja et al., 2000) and its multivariate version M-AMBI (Muxika et al., 2007), the Benthic Index (BENTIX) (Simboura and Zenetos, 2002), the Benthic Quality Index (BQI) (Rosenberg et al., 2004) and the BQI-Family (Dimitriou et al., 2012). As reported in numerous publications (Borja et al., 2003(Borja et al., , 2008Simboura and Reizopoulou, 2007;Gremare et al., 2009;Subida et al., 2012;Karakassis et al., 2013;Moraitis et al., 2013;Simboura et al., 2016), these indices have been found to be suitable tools for assessing the Ecological Status (ES) of coastal marine ecosystems.
Organic enrichment-induced changes in benthic species composition may also affect the functional diversity of benthic communities (Brown et al., 1987;Solan et al., 2004;Tillin et al., 2006). Functional diversity refers to the range of organismic characteristics affected by ecosystem processes (Hooper et al., 2002). Changes in the functional properties of benthic communities represent the adaptation of the organisms to the environment and their response to disturbance and pressure (de Juan et al., 2007). In healthy communities, the functional roles of individual species are often complementary rather than competitive (Rice et al., 2012). In the case of eutrophication, the response of benthic fauna to organic enrichment depends on the biological traits of the organisms. Different species have different activity patterns and the importance of faunal activities for ecosystem maintenance is closely related to individual functional traits (Papageorgiou et al., 2009). A moderate increase of organic matter in the sediment has been reported to increase macrofaunal biomass by favoring deposit-feeding organisms that take advantage of the available food (Grall and Chauvaud, 2002), thus reducing the chance of hypoxia occurring. Species and functional diversity have been reported to be strongly inter-correlated (Hewitt et al., 2008), and Biological Traits Analysis (BTA) (Bremner, 2008) has proven to be a powerful tool for assessing anthropogenic impacts on ecosystem functional diversity (Beauchard et al., 2017).
Macrofaunal communities influence biogeochemical processes, especially in soft bottom sediments through bioturbation, i.e., the process of modification of sediments through particle reworking and burrow ventilation carried out during foraging, feeding and burrow maintenance activities (Volkenborn et al., 2010;Queirós et al., 2013). Bioturbation also affects sediment oxygenation, pH, redox potential (Eh), and other variables (see Queirós et al., 2013) and therefore reduces the negative effects of organic enrichment on marine sediments (Martinez-Garcia et al., 2015). A measurement of the bioturbation potential of a benthic community can be provided by the index Community Bioturbation Potential (BPc) proposed by Solan et al. (2004). This index combines abundance and biomass data as well as the functional traits of the macrofaunal community in order to provide a quantification of the bioturbation activity in a particular study area.
All the above-mentioned methodologies have been widely applied in filed data. However, in order to study the effects of specific pressures on biological responses, a controlled experimental environment may be more suitable (Mostajir et al., 2013). An experimental design where an adequate volume of water and sediment is isolated and exposed to specific (or a combination of) pressures is ideal in order to study the interactions between organisms and their physical environment (Petersen et al., 2009). Such enclosed ecosystems, called mesocososms, have been used in a wide range of studies but mainly in those focusing on the water column (see Dimitriou et al., 2017b). Typically, a mesocosm should contain a large volume of water (typically >1.45 m 3 ) with sufficient depth in order to allow the study of biological responses (Petersen et al., 2009). Adding a sediment compartment with a fully functioning benthic community allows the study of the whole coastal ecosystem, enabling the creation of doseresponse curves between the water column and benthos. Largescale mesocosm experiments focusing on eutrophication in the water column have been performed in numerous cases (see Dimitriou et al., 2017a). However, only a few of those studies include a sediment compartment as they mainly use smaller-scale boxcosms (Olsgard et al., 1998;Suomela et al., 2005;Widdicombe et al., 2009;Valdemarsen et al., 2010). A common finding in those studies is that bioturbation plays a very important role in the resuspension of nutrients and sediment metabolism. However, none of these studies focus on changes in macrofaunal community structure, biodiversity, and ES, and directly relating them to the pressure of eutrophication.
The application of well-established methods and tools used in field studies in an experiment is a novel approach to describe in detail the response of the macrofaunal community to eutrophication. Studying only biodiversity is no longer enough for understanding ecosystem health and response to anthropogenic pressures (Tett et al., 2013). Also, ecosystem resilience is more related to the functional traits of the benthic community rather than to biodiversity itself. Therefore, combining "traditional" biodiversity methods and ES assessment with BTA and bioturbation potential in a controlled experimental manipulation should provide reliable results for assessing the effects of marine coastal eutrophication on the benthic ecosystem.
The aim of the present study was to assess the effects of eutrophication on the macrofaunal community during a benthicpelagic coupling mesocosm experiment in the oligotrophic eastern Mediterranean region. More specifically, we wanted to quantify the dose-response relationship between Organic Matter (OM) sedimentation caused by various levels of eutrophication and the change of the benthic community species composition, ES and ecosystem functions after a 2-month mesocosm experiment.

Experimental Setup
The mesocosm experimental setup used in the present study took place at the CRETACOSMOS mesocosm facilities of the Hellenic Center for Marine Research in Crete, in the eastern Mediterranean (www.cretacosmos.eu). It included nine mesocosms (0.35 m radius) containing a water column of 4 m depth and 1.5 m 3 total volume, and at the bottom, undisturbed sediment of 0.3 m radius, 0.3 m height, and 85 l total volume. The experimental design included two levels of nutrient addition, with one single addition of dissolved nutrients-based on the ambient concentration of nutrients-on the first day of the experiment (codenamed day 0): 100x for PO 3+ 4 (5 µM) and 300x for NO − 3 (30 µM) for Low eutrophication treatment (codenamed "Low") and 200x PO 3+ 4 (10 µM) and 600x NO − 3 (60 µM) for High eutrophication treatment (codenamed "High"), with the aim of creating eutrophic conditions representative of the eastern Mediterranean region. Additionally, the experiment included triplicate control mesocosms (codenamed "Control") with no nutrient addition.
The water used in the experiment was pumped from the Hellenic Center for Marine Research coastal area in Heraklion, Crete, Greece (southern Aegean Sea, eastern Mediterranean) from 2 m depth and was immediately transferred into the mesocosms. The sediment was collected from an area of reduced anthropogenic impact within the harbor area of the port of Heraklion, Crete, where the water conditions become periodically eutrophic, ranging from Poor to Moderate ES according to the scale based on Chl-a concentration proposed by Simboura et al. (2005) and in which no other anthropogenic pressure is important. Triplicate boxcorer samples (0.017 m 2 surface) were taken for the determination of ambient macrofauna community status. Detailed information about the experimental setup can be found in Dimitriou et al. (2017b), while data on water column eutrophication, sedimentation rate and effects on sediment geochemical variables are presented in Dimitriou et al. (2017a).
At the end of the experiment, one boxcorer sample (0.017 m 2 surface) was taken from each mesocosm for the determination of the final status of the benthic community. All macrofaunal samples were sieved using a 0.5 cm sieve mesh size, preserved in 5% Formalin and stained with Rose Bengal to facilitate sorting and identification. Animals were then identified to the species level or to the lowest taxonomic level possible and then weighed for the determination of the total wet biomass of each species.

Statistical Analysis
Macrofaunal species data were used for a number of multivariate analyses including non-metric Multidimensional Scaling (nMDS), ANOSIM, and SIMPER analysis using PRIMER-E ver. 7 (Clarke and Gorley, 2015). The PRIMER-E software was also used to calculate a number of established diversity indices such as total number of species (S), total number of individuals (N), Hurlbert's expected number of species [ES(10)], log 10 Shannon index (H ′ ) and Pielou's (J ′ ). In addition, a series of benthic indicators for ES was calculated using the formulas provided by the authors in the respective publications. For BENTIX (Simboura and Zenetos, 2002) and BQI-Family (Dimitriou et al., 2012), the Biological Indices Calculation Tool (BICT) online calculation tool hosted on the Lifewatch Greece website (http://www.lifewatchgreece.eu) was used with the species inventory provided at http://www.hcmr.gr/en/the-bentixindex/ and at (http://www.sciencedirect.com/science/article/pii/ S1470160X12000544#MMCvFirst), respectively. Additionally, the BPc index (Solan et al., 2004) was calculated using the methodology and species inventory data provided by Queirós et al. (2013). Two-way Repeated Measures Analysis of Variance (RM-ANOVA) was performed using SPSS v22 in order to detect significant differences between treatments (ending point) and through time of all the above-mentioned indices and indicators (assumptions of normality were assessed using the Shapiro-Wilk, homogeneity using the Levene test and sphericity using the Mauchly's sphericity test). In case of significant two-way treatment and time interaction, simple effects of treatment (comparison between values of each treatment at day 58) and time (comparison between values at days 0 and 58 for each treatment separately) were assessed.
Biomass per individual −1 for each species was compared between the ambient conditions and the end of the experiment by dividing the biomass of each species for each treatment with the corresponding ambient value. If the ratio was smaller than 100%, it was indicative of weight loss while a ratio larger than 100% indicated weight gain.
The selection of functional traits for BTA should be determined based on their relevance to research objectives (Beauchard et al., 2017). In the present study, the impact of eutrophication-induced sedimentation on benthos was the most important factor considered for trait selection, therefore selected traits included Life Strategy, Body Size, Feeding Type, Age at First Reproduction, and Ecosystem Engineering. Life Strategy provided an estimation of a species' tolerance to disturbance as characterized by the benthic indicators AMBI, BENTIX, and BQI in the corresponding publications. Body Size was based on the maximum body size of each species; Feeding Type indicated the feeding habits of each species; Age at First Reproduction indicated the time required for each species from hatching to first reproduction; and Ecosystem Engineering indicated the type of ecosystem engineering of each species. For the determination of functional traits and modalities of macrofaunal species, the database POLYTRAITS (http://polytraits.lifewatchgreece.eu/) (Faulwetter et al., 2014), which contains information about functional traits for more than 950 different species along with the relevant publications, was used.

RESULTS
The nMDS plot of macrofaunal data (Figure 1) indicates that the samples are clustered into three groups: On the top right side of the plot, the 3 Ambient samples (ambient conditions) are grouped, on the left side the three Control samples are grouped together with the Low treatment samples, while on the bottom right side of the plot, the High treatment samples are grouped together. ANOSIM results (Table 1) indicate that there are significant differences between samples (Global R: 0.77, p < 0.01) while pairwise tests showed significant differences between all pairs of samples except between Control and Low treatment.
The values of the diversity and ES indices (Figure 2) indicate that in all cases, the High treatment samples showed major differences compared to the other samples in respect of increases in the total abundance and the dominance of certain species while reducing diversity as well as the values of the ES indicators. RM-ANOVA results (Supplementary Material Table 1) indicate that in all cases [except for the total number of species (S) and families (F)], the change was statistically significant (p < 0.01) ( Table 2).
The results of SIMPER analysis ( Table 3) indicated that five polychaete species were responsible for the greater part of the dissimilarity between the groups of samples. The species Paradoneis lyra (Fam. Paraonidae) was responsible for the largest percentage of dissimilarity, followed by Leiochone leiopygos and Chaetozone setosa (Fam. Cirratulidae), Monticellina dorsobranchialis (Fam. Maldanidae), and Notomastus latericeus (Fam. Capitellidae). Out of the five species, P. lyra exhibited the largest increase in abundance, starting from 65 ind per 0.1 m −2 and reaching 116 in Control, 156 in Low treatment and 350 ind per 0.1 m −2 in High treatment. In contrast, the species N. latericeus was the only one that presented a reduced abundance at the end of the experiment in all samples compared to the initial conditions.
Those species responsible for the greater dissimilarities between samples also changed their biomass per individual −1 ( Table 4). In the High treatment, all five species increased their biomass per individual −1 compared to the ambient samples. In contrast, in the Control, all five species exhibited a reduction in their biomass per individual −1 , while in the Low treatment three species showed an increase and two showed a decrease in biomass per individual −1 . The BTA of those five macrofaunal species (Table 5) showed that they have similar ecosystem functions. They are transitional or opportunistic species with medium or small body sizes, deposit feeders and show high contributions to the bioturbation potential of the community. A difference in the age at first reproduction was noted, with P. lyra and N. latericeus having fast reproduction cycles of <2 months, unlike the rest of the species.
The results of the BP c index (Figure 3) showed that the bioturbation potential of the benthic community of the High treatment was significantly higher (p < 0.01) while the corresponding of the Control was significantly lower (p < 0.05) (Supplementary Material Table 1).  Table 1.

DISCUSSION
Mesocosm experiments have proven to be powerful tools for studying the response of aquatic ecosystems to specific pressures.
The results of the present study indicated that the increased pressure caused by eutrophication on the benthic macrofauna (in the High treatment) resulted in a significant alteration of the community in terms of biodiversity, ES and functional characteristics.
In Dimitriou et al. (2017a), it was reported that in both nutrient-enriched treatments, OM sedimentation caused by eutrophication exerted a strong pressure on the benthos, inducing changes in sedimentary environmental variables.
In the middle of the experiment (days 18-36), OM, Total Organic Carbon and Sulfide emissions increased while Eh-values decreased, forming a less suitable environment for benthic macrofauna to inhabit. However, in the Low treatment, conditions in the sediment had improved by the end of the experiment, returning close to the initial point. In the High treatment, in contrast, eutrophic conditions in the water column and high Particulate Organic Carbon sedimentation rate persisted throughout the experiment. While the benthic community of the Low treatment did not seem to be affected (no significant variation from the Control was observed), the community structure of the High treatment altered significantly. The new community was characterized by lower diversity, a greater dominance of specific species and a shift from "Good Environmental Status (GES)" to "Unacceptable" conditions, as indicated by the 3 benthic indices.
As OM settles on the sediment surface, benthic infauna contributes to its multidimensional distribution between the different layers of the sediment through feeding or tube construction and its recycling through the process of digestion and defecation (Heilskov et al., 2006). Solan et al. (2010) argue that different benthic species affect ecosystem functions differently by not having the same ecosystem engineering capabilities. Of the 23 species found in the samples, SIMPER analysis identified five that were responsible for the greater dissimilarity between the different levels of nutrient addition. BTA analysis, aimed to relate the functional profile of those species to the environmental condition of each treatment, showed that all five species exhibited common characteristics-conserving their life strategy and being tolerant  or opportunistic species-according to the scale proposed by Glémarec and Hily (1981), with a similar size and feeding type through deposit feeding, which, according to Weigelt (1991), is favored in eutrophic conditions, and all were "ecosystem engineers." However, only two species had life cycles fast enough to be able to reproduce within the mesocosms during the duration of the experiment. Furthermore, the species P. lyra had significantly higher starting abundances (65 vs. 10 individuals per m 2 ). As a result, P. lyra competed and could increase its abundance, dominate the community and eliminate N. latericeus from the Control and Low treatments, where food availability was limited. In contrast, in the High treatment where food was abundant, N. latericeus abundance was not affected.
In a benthic mesocosm experiment in Rhode Island, Oviatt et al. (1986) reported that adding 32x nutrients in the water column increased benthic respiration, causing the system to become temporarily anoxic at regular intervals. In the current experiment, highly eutrophic conditions were created, especially in the High treatment where the ES of the water column was "Bad" throughout the experiment (Dimitriou et al., 2017a). Despite the increased phytoplankton-generated OM sedimentation and the resulting increased OM content of the sediment, even in the highly eutrophic conditions of the High treatment, no hypoxia in the benthic boundary layer or in the sediment was observed. Dissolved oxygen concentration remained above the ambient values throughout the experiment (Dimitriou et al., 2017a). Even if the benthic boundary layer is well-oxygenated, due to OM decomposition in the sediment, a reduction in the Eh-values would be expected based on the results of other studies concerning organic enrichment, i.e., in aquaculture (Karakassis et al., 2000;Hyland et al., 2005) or the effects of eutrophication in the eastern Mediterranean (Dimitriou et al., 2015). However, in this experiment, the lowest recorded Eh-value was −8 mV in the High treatment near the end of the experiment and the values quickly improved (Dimitriou et al., 2017a). In this case, the macrofaunal community played a very important role in maintaining oxic conditions: All species indicated by SIMPER analysis, and most importantly P. lyra, which reached very high abundances in the High treatment, were ecosystem engineers (modifiers or upward conveyors). As indicated by the BP c index, the benthic community of the High treatment had a significantly higher bioturbation potential. High bioturbation potential, combined with an increased abundance of deposit-feeding species (again, all five species were deposit feeders), contributed to the OM consumption and oxygenation of the sediment. On the contrary, in the Low treatment, the macrofaunal community was not affected significantly as there were no substantial changes in diversity, ES and bioturbation, with only an increase in the biomass of some species, as also reported by Widbom and Frithsen (1995). The relationship between nutrient enrichment and bioturbation is also highlighted in different studies. Suomela et al. (2005) reported that in a eutrophication benthic mesocosm experiment in the Baltic Sea, the absence of tube-forming macrofaunal species reduced bioturbation while the high abundance of suspension-feeding bivalves contributed to increased nutrient flux, especially ammonium, from the sediment. Valdemarsen et al. (2010), in a smaller-scale benthocosm experiment, measured bio-irrigation and found increased rates in nutrient-enriched treatments.
As stated by Griffin et al. (2009) and Culhane et al. (2014), the inclusion of functional traits in ecological studies contributes to the better understanding of ecosystem processes. Benthic communities with a GES are normally characterized by a few abundant species and many rare ones, and the roles played by individual species are often complementary rather than competitive (Rice et al., 2012). The benthic community used in the experiment originated from an area of limited anthropogenic impact and the ES classification was "Good" in respect of all the benthic indicators used, with the water column becoming both locally and periodically eutrophic. Such communities show high resilience potential to pressure because biodiversity buffers ecosystem processes by altering species composition and functions (Rice et al., 2012). Under the pressure of increased sedimentation of phytoplankton-generated OM, sediment geochemical variables showed irregular succession patterns, especially in the High treatment (Dimitriou et al., 2017a). Furthermore, benthic fauna responded by increasing the abundance and the biomass of species capable of consuming OM from the sediment, as opposed to the Control where no such increase was recorded. Despite the fact that the new benthic community was characterized by lower diversity and a shift to "Moderate" ES, it had a higher bioturbation potential, reducing the detrimental effects of eutrophication while increasing sediment oxygenation. Tett et al. (2013) argued that "healthy ecosystems show resilience to external pressures and can be maintained without the need of human intervention." They contain adaptation mechanisms and the necessary vigor to maintain ecosystem balance. In the face of a strong but short-term disturbance event caused by eutrophication, as in this experiment, deposit feeding and bioturbation were used as feedback mechanisms that prevented the collapse of the benthic communities, thereby maintaining the ability of the benthic system to provide its basic ecosystem services despite the significant changes in sediment geochemistry.