First Field-Based Evidence That the Seagrass-Lucinid Mutualism Can Mitigate Sulfide Stress in Seagrasses

Seagrass meadows form vital ecological components of coastal zones worldwide, but are rapidly declining. Large-scale seagrass diebacks have been related to accumulation of toxic sulfide in the sediment, a phenomenon predicted to occur more frequently in the near future due to ongoing global warming and increasing organic loading of coastal systems worldwide. Recently, a facultative mutualism between seagrasses and lucinid bivalves with endosymbiotic sulfide-consuming gill bacteria was discovered that may prevent toxic sulfide accumulation in seagrass sediments. Yet, direct field-based evidence for the importance of this mutualism in alleviating sulfide stress in seagrasses is currently lacking, as well as how its role may change when sediment sulfide levels increase due to environmental change. Here, we investigated the sulfide detoxification function of this seagrass-lucinid mutualism and its resilience to organic-loading induced sulfide stress in a temperate lagoon system (Thau lagoon, France), using a correlative field survey and a full factorial field experiment. The field survey revealed a strong positive correlation between seagrass above-ground biomass and lucinid densities, and pore water sulfide concentrations close to zero at all sites. Furthermore, the field experiment revealed that addition of organic matter (starch mixed with sucrose) increased sedimentary sulfide intrusion into seagrass (Zostera noltei) leaves (a proxy for sulfide stress), while experimentally enhanced lucinid densities reduced sulfide intrusion, regardless of addition of organic matter. Moreover, addition of organic matter reduced seagrass rhizome biomass and increased pore water sulfide levels, lucinid tissue sulfur content, lucinid condition (expressed as flesh/shell dry weight ratio), and total lucinid biomass, while enhancement of lucinid densities reduced lucinid condition. These results provide the first field-based evidence that lucinid bivalves and their sulfide-oxidizing gill symbionts mitigate sulfide stress in seagrasses, and suggests that the dependence of seagrass on this seagrass-lucinid mutualism will increase under conditions of enhanced sediment sulfide production, as predicted for the near future. Therefore, we suggest that awareness of the ecological importance of the seagrass-lucinid mutualism may be instrumental for designing new measures for improving long-term restoration success and seagrass resilience to global change.


INTRODUCTION
Seagrass meadows rank amongst the ecosystems providing the greatest value in terms of ecological importance and ecosystem services including enhanced biodiversity, carbon storage, nutrient cycling, and coastal protection (Costanza et al., 1997;Hemminga and Duarte, 2000;Orth et al., 2006). Yet, despite their ecological and societal importance, seagrass meadows are rapidly declining worldwide (Waycott et al., 2009), with losses often characterized by sudden collapse. Ecological studies from a wide range of terrestrial, freshwater and marine ecosystems suggest that such rapid declines are often associated with the disruption of internal positive feedback mechanisms due to environmental changes (Scheffer et al., 2001;Kiers et al., 2010). Indeed, mounting evidence suggests that seagrass meadows are often mediated by multiple self-reinforcing feedback mechanisms, that can cause non-linear responses to environmental changes, and even result in alternative stable states if they are strong enough (bistability) (Van Der Heide et al., 2007;Carr et al., 2010;Maxwell et al., 2017). Yet, despite these recent insights, both protection and restoration of seagrass meadows remain extremely difficult (Orth et al., 2006;Van Katwijk et al., 2016), suggesting that certain fundamental processes or interactions may be overlooked.
Indeed, a mutualistic sulfide-detoxification feedback was recently discovered in laboratory experiments that appears to be present in many seagrass meadows worldwide (Van Der Heide et al., 2012) and that could potentially be a key interaction mediating both seagrass productivity and ecosystem stability (De Fouw et al., 2016a, 2018. Seagrass meadows facilitate their own growth by trapping suspended particles and stabilizing sediments, which increases water clarity (Van Der Heide et al., 2007;Folmer et al., 2012). This self-reinforcing feedback, however, also results in a negative feedback as organic matter accumulates in the sediment, which is primarily decomposed by sulfate-reducing bacteria with toxic sulfide as a metabolic end-product (Jørgensen, 1982;Marbà et al., 2006). Although seagrasses transport oxygen to their belowground tissues (radial oxygen loss, ROL) to maintain an oxic microsphere around roots and rhizomes (Pedersen et al., 2004;Lamers et al., 2013;Brodersen et al., 2015), sulfide production can outpace oxygen release in organic matter-rich sediments and under warmer conditions (Calleja et al., 2007;Koch et al., 2007). This can eventually result in sedimentary sulfide intrusion into plant tissues (i.e., sulfide stress), with negative effects on photosynthesis, growth and seagrass survival (Frederiksen et al., 2007;Mascaró et al., 2009;Garcìa et al., 2013;Borum et al., 2014;Kilminster et al., 2014). However, a recent laboratory experiment revealed that seagrasses can create a positive feedback by engaging in a mutualistic interaction with lucinid bivalves and their endosymbiotic sulfide-oxidizing gill bacteria to reduce sulfide stress (Van Der Heide et al., 2012). In return, the lucinid bivalves and their endosymbionts profit from the seagrasses not only due to enhanced sulfide production from the organic matter trapped by the plants, but also from oxygen released by seagrass roots (Van Der Heide et al., 2012). Hence, while the laboratory experiment revealed that this facultative mutualism can strongly enhance seagrass productivity (Van Der Heide et al., 2012), observational and theoretical investigations in a West African tropical intertidal seagrass system suggested that breakdown of this mutualistic feedback has amplified seagrass collapse after a severe drought event (De Fouw et al., 2016a, 2018. These recent studies suggest that this mutualistic sulfidedetoxification feedback could be a key interaction mediating both seagrass productivity and ecosystem stability. Yet, the ecological importance of this sulfide-detoxification mutualism, as well as how its role may change in the face of environmental change, remain to be experimentally tested under field conditions. Because eutrophication and global warming are predicted to cause higher sulfide production rates in coastal shallow sediments in the near future (Garcìa et al., 2013;Lamers et al., 2013), such knowledge is essential, not only to validate the potential of using this sulfide-detoxification mutualism for seagrass conservation and restoration purposes, but also to predict future scenarios based on eutrophication and climate change prospects.
Hence, in this study we aim to (1) assess whether the sulfide-detoxification mutualism can mitigate sulfide stress in seagrasses under field conditions, and (2) test whether mutualism-dependency of seagrass (and simultaneously the risk of mutualism breakdown) will increase under conditions of enhanced sediment sulfide production. We use the subtidal seagrass beds dominated by Zostera noltei in Thau lagoon (France), which are known to be inhabited by the lucinid bivalve Loripes orbiculatus (synonyms L. lucinalis and L. lacteus) (Rossi et al., 2013). Like all lucinid bivalves investigated so far, L. orbiculatus mainly lives off carbon products provided by sulfide-oxidizing chemoautotrophic bacteria [i.e., Candidatus Thiodiazotropha endoloripes (Petersen et al., 2016)] living inside its gills . In this symbiotic association, the bivalve host favors chemosynthesis by its gill-symbionts by facilitating the supply of sulfide, carbon dioxide and oxygen to its gills. In exchange, the bacterial gill-symbionts fix inorganic carbon, fueling their own energetic and biosynthetic needs, in addition to those of their host (Stewart et al., 2005).
As a first exploration of the ecological importance of the seagrass-lucinid mutualism in our study system, we carried out a field survey to examine in situ pore water sulfide levels and the correlation between seagrass biomass and L. orbiculatus densities. Next, we investigated in situ how natural vs. experimentally enhanced L. orbiculatus densities affected sedimentary sulfide intrusion into Z. noltei leaf tissues (a proxy for sulfide stress) and L. orbiculatus performance, in sediments with either natural or experimentally enhanced sulfide production rates (induced by addition of organic matter). We hypothesized that (i) under field conditions, the seagrass-lucinid mutualism is apparent and mitigates sulfide stress in seagrasses, which should be reflected by a positive correlation between seagrass biomass and lucinid densities (and vice versa) and overall low pore water sulfide concentrations, and (ii) that sedimentary sulfide intrusion would increase in plots with experimentally enhanced sulfide production rates (induced by addition of organic matter), which in turn (iii) would be mitigated in plots where we also experimentally increased lucinid densities.

Study System
Both the field survey and the field experiment were conducted in the subtidal Zostera noltei-dominated seagrass meadows of Thau lagoon, which are known to be inhabited by the lucinid bivalve Loripes orbiculatus (Rossi et al., 2013). This lagoon is located on the French Mediterranean coast (Figure 1) and covers an area of about 69 km 2 and has a mean depth of 4 m. Seawater inputs from the Mediterranean Sea enter by artificial inlets. Precipitation shows large variation between years, ranging from 200 to 1300 mm year −1 (Collos et al., 2009). Like for most lagoon systems, nutrient concentrations in Thau lagoon depend mainly on rainfall events and watershed runoff, with nitrate and soluble reactive phosphorus ranging between undetectable (<0.05 µM) and 22 µM and undetectable and 10 µM, respectively (Collos et al., 2009). Seawater temperature ranges between 4 • C in February and 29 • C in August, and salinity ranges between 29 in February and 40 in August (Collos et al., 2009).
Events of sudden seagrass die-offs in Thau lagoon have been reported up to the late 1990's (Plus et al., 2003b). These die-offs were attributed to toxic sulfide accumulation in the sediment and water column due to the anaerobic degradation of green algae and probably organic matter from aquaculture, accelerated by high temperatures (Souchu et al., 1998;Chapelle et al., 2001).

Field Survey
In late spring (between 28 April 2016 and 15 June 2016), we conducted a field survey at Thau lagoon (France) to test the strength of the relation between seagrass biomass and lucinid bivalve density and vice versa. Seven sampling sites (Figure 1) were chosen, which aimed to be representative for the variable conditions experienced by shallow (0.3-0.8 m depth range) Z. noltei meadows within the lagoon. At each sampling site, we haphazardly selected at least 2 sampling stations (n total = 21) to be able to predict within-site variation in response variables. At each sampling station a sediment core (internal diameter 10 cm) was taken to a depth of 20 cm and the content was sieved over a 1 mm mesh sieve. Lucinid bivalves and seagrass shoots were sorted from the material retained on the sieve. Lucinid bivalves, all of which were identified as Loripes orbiculatus, were counted and shell length (precision 0.1 mm) was determined for each clam. Next, seagrass shoots were scraped for epiphytes and divided into above and below-ground biomass and dried for 72 h at 60 • C to constant weight.
In addition, at each sampling station a sediment pore water sample was collected by suction, using 60-ml vacuumed syringes connected to ceramic rhizon samplers (Eijkelkamp Agrisearch Equipment, Giesbeek, Netherlands). The rhizon samplers were placed at a depth of 5 cm in the sediment (the observed average rooting depth of Z. noltei). For 5 of the 7 sampling sites, i.e., S1, S2, S5, S6, and S7 (see Figure 1B) we succeeded in preventing oxygen contamination during transport and were able to measure free sulfides in the pore water within 6 h after sampling. This was measured in a mixture of 50% sample and 50% Sulfide Anti-Oxidation Buffer (SAOB), using a calibrated ion-specific silver-sulfide electrode (detection limit 0.1 µM).

Field Experiment
A field experiment was developed to test the importance of the sulfide-detoxification mutualism under field conditions and how its role may change under experimentally enhanced sulfide production. Densities of L. orbiculatus and pore water sulfide production (stimulated by organic matter addition) were manipulated in situ in a 2 × 2 factorial design, and their effect on sedimentary sulfide intrusion into Zostera noltei leaf tissues (a proxy for sulfide stress) was monitored, using stable sulfur isotope analysis . The field experiment was carried out in summer between 27 June and 16 August 2016 (i.e., 50 days) when sulfide production rates are expected to be high due to relatively high seawater temperatures and when Z. noltei growth rates are high (Plus et al., 2003a). As the leaf plastochron interval (time period between production of consecutive new leaves from the same meristem) for Z. noltei is typically 10-12 days (Brun et al., 2002;Peralta et al., 2002), the 50day duration of our experiment allowed the treatments to induce effects on isotopic ratios in newly formed leaf tissues in this period. The field experiment was conducted at a homogeneously Z. noltei-vegetated area covering 200 m 2 (water depth ∼ 0.5 m) at sampling site S5 (see Figure 1). Here, we haphazardly selected 6 blocks of 1.5 m 2 (>5 m apart). At the onset of the experiment on 27 June 2016, we first took a sediment sample (diameter 2.5 cm) in the middle of each block to a depth of 7 cm for analyses of initial sediment characteristics in the rhizosphere. Within each block, we then allocated four experimental plots (diameter 30 cm) spaced 0.9 m apart, to which we randomly applied two treatments, addition of lucinids bivalves and addition of organic matter, according to a 2 × 2 factorial design. This resulted in four experimental groups: control (C), only lucinids added (L), only organic matter added (OM), and both lucinids and organic matter added (L + OM). Organic matter additions were realized by inserting custom-made OM sticks into the sediment allowing the delivery of labile organic matter. The OM sticks (20 g each) contained 2 parts starch and three parts sucrose and water and were baked in silicon trays at 100 • C for 1 h (modified from Govers, 2014). In each OM and L + OM-plot we inserted in total 9 OM sticks in the sediment (every 10 × 10 cm, to a depth of ∼7 cm). The organic matter was only reactive when it dissolved from the sticks; we thus created a slow-release carbon source (Govers, 2014). To compensate for the physical disturbance of placing the OM sticks in the sediment, we applied a similar disturbance to the C and L-plots, with the only difference that we used nine inert sticks that were removed immediately after insertion. Our reasoning to remove these inert sticks was that it would better mimic the addition of organic matter treatment, because after some time the inserted OM sticks in the OM and L + OM-plots would also be no longer present due to dissolution.
For the lucinid bivalve addition treatment, we haphazardly collected 1200 L. orbiculatus clams from a Z. noltei patch at ∼25 m distance from our experimental blocks. These clams were equally divided over 12 buckets, such that size-frequency distributions of the 100 clams per bucket were similar to those in the field, which resulted in a total clam wet weight of 32 g in each bucket. Next, at each experimental plot to which a L or a L + OM treatment was randomly assigned, we emptied a bucket with 100 L. orbiculatus clams (thus increasing natural L. orbiculatus densities within our plots with 1415 individuals per m 2 ).

Water, Pore Water and Sediment Handling and Analysis
After 49 days, on 15 August 2016, at each experimental plot (n = 24), we collected sediment pore water samples (60 ml) as described above between 11:15 AM and 1:15 PM. A subsample of the pore water sample was analyzed within 1 h after sampling for total sulfide concentration as described above. The remainder of the pore water sample (20 ml) was transferred to a vial for analysis of nutrient concentrations, to investigate whether our treatments had any unforeseen effects on pore water nutrient concentrations.
Next, on 16 August 2016 (i.e., after 50 days), at each experimental plot a sediment core (internal diameter 2.6 cm) was taken to a depth of 7 cm for stable sulfur isotope analyses after which a 15 ml subsample was transferred into a 50 ml tube with 15 ml of 1 M ZnAc (vol:vol). Subsequently, seagrasses and lucinid bivalves were collected by taking a benthic core (internal diameter 10 cm) to a depth of 20 cm in the center of each experimental plot, which was sieved over a 1-mm mesh sieve, after which all material retained on the sieve was transferred to a plastic bag. In addition, duplicate GF/F filtered-bottom water samples of 20 ml each were collected for stable sulfur isotope analysis. All samples were kept refrigerated during transportation to the laboratory where they were preserved frozen until analysis.
Pore water ammonium and ortho-phosphate concentrations were measured colorimetrically (Skalar and Seal autoanalyzer), using ammonium-molybdate and salicylate. Nitrate was determined by sulfanilamide, after reduction of nitrate to nitrite in a cadmium column (Wood et al., 1967). All nutrients were measured at the analytical lab of the Radboud University in the Netherlands.
The water samples collected for stable sulfur isotope analyses (δ 34 S water ) were prepared by boiling under acidic conditions followed by precipitation of sulfate with BaCl 2 as BaSO 4 . Next, δ 34 S was measured at the University of Southern Denmark, by weighing grinded BaSO 4 (and reference material) into tin capsules together with vanadium pentoxide, which were analyzed using an EA-IRMS (Delta V Advantage Isotope Ratio MS with Thermo Scientific EA). The stable isotope signatures were reported in standard delta notation (units per mill, ) relative to the international standard for δ 34 S of Canyon Diablo Troilite (CDT).
The rhizosphere sediment collected at each block at the onset of the experiment was freeze-dried, gently homogenized, and sieved over a 700-µm mesh. Next, sediment organic matter content was estimated as mass loss by ignition at 550 • C for 4 h, while grain size distribution was measured by laser diffraction on a Beckman Coulter Model LS13320 particle size analyzer (Geosciences Montpellier). From the grain size distribution, the mean grain size (MGS) was calculated. Total iron (total Fe) was obtained by boiling combusted sediment in 1 M HCl and analyzed according to Stookey (1970) and Sørensen (1982). The carbonate content was estimated from the weight loss of sediment after acid digestion. Based on these analyses, the initial sediment conditions within the experimental blocks can be characterized by sediment with mean grain size of 201.6 ± 5.9 µm (mean ± SE), organic matter content of 3.4 ± 0.1%, carbonate content of 38.7 ± 1.0%, and iron content of 30.0 ± 0.3 µmol g −1 .
The sediment collected for sulfur stable isotope analysis was distilled according to the two-step procedure of Fossing and Jørgensen (1989), but only the chromium reducible sulfur (CRS) fraction consisting of FeS 2 and S 0 was used, as too little material was retrieved in the acid volatile (AVS) fraction. Briefly, an aliquot of 5-10 g of homogenized ZnAC-treated sediment was transferred to a reaction flask and 10 ml 50% ethanol was added. After degassing with N 2 for 10 min, the sediment was distilled at room temperature with 8 ml 12 M HCl for 30 min to release AVS. CRS was obtained by reduction with 16 ml Cr 2+ in 0.5 M HCl followed by boiling and distillation for 30 min. The extracted sulfide was precipitated in traps as Ag 2 S and analyzed for stable sulfur isotopes (δ 34 S CRS ) as described above for δ 34 S water .

Seagrass and Lucinid Bivalve Handling and Analysis
Zostera noltei shoots and Loripes orbiculatus clams were sorted from each benthic sample. Shoot density was determined and epiphytes were removed from each shoot by gently scraping with a glass slide. Next, shoots were rinsed in deionized water and divided into leaf, rhizome and root biomass and dried for 72 h at 60 • C and weighed to the nearest 0.1 mg. Per benthic sample, the youngest leaves (i.e., two pairs of the most internal leaves) of two apical shoots were subsampled, freeze-dried and analyzed for stable sulfur isotopes (δ 34 S leaves ) as described above for δ 34 S water .
Only the youngest Z. noltei leaves were used for this purpose, as they typically represent the latest ∼2 weeks of growth (Brun et al., 2002;Peralta et al., 2002) so that their δ 34 S signals reflected the sulfide conditions experienced by Z. noltei during the 50-day experimental period (Holmer and Kendrick, 2013). Total sulfur concentration (% of dry weight) in leaf tissue (TS leaves ) was also obtained from this analysis.
The percentage of sulfur in the plant tissue originating from sediment sulfide (F sulfide ), i.e., sulfide intrusion, was calculated according to Frederiksen et al. (2006): where δ 34 S leaves is the value measured in the leaf, and δ 34 S water and δ 34 S CRS are the values measured in bottom water samples and sediment samples (CRS pools), respectively. The amount of sulfur incorporation (% of dry weight) into the leaves from sedimentary sulfides (SS leaves ) was then calculated as follows: Per benthic sample, L. orbiculatus clams were counted and shell length of each clam was measured to the nearest mm. Flesh and shell were separated and dried for 72 h at 60 • C and weighed to the nearest 0.1 mg. For some L. orbiculatus clams of which shells were broken (42 out of 510 clams) we could not determine shell dry weight. L. orbiculatus condition was expressed as flesh/shell dry weight ratio, which is a commonly used size-and-age independent measure of fitness in bivalves (Lucas and Beninger, 1985). Indeed, L. orbiculatus condition was independent of L. orbiculatus shell length (F 1,460.43 = 0.46, P = 0.50).
Per benthic sample, we haphazardly selected a minimum of two and if available three L. orbiculatus clams >7 mm in shell length, from which (dried) flesh was freeze-dried and grinded to measure sulfur stable isotope (δ 34 S Loripes ) and total sulfur content (TS Loripes ) as described above for seagrass.

Statistical Analysis
Due to the nested structure of the data collected during our field survey with multiple benthic samples taken per sampling site, repeated measures correlation (r rm values with 95% confidence intervals) were calculated for each sampling site using the "rmcorr" R package (Bakdash and Marusich, 2017) to determine within-site correlations between seagrass above-and belowground biomass, seagrass above-ground biomass and lucinid bivalve density, and seagrass below-ground biomass and lucinid bivalve density.
For data obtained during our field experiment linear mixed effects models with random intercepts and with Block and, where applicable, sampling core (nested in Block) as random effects, were used to investigate whether addition of L. orbiculatus clams, addition of organic matter, and the interaction between these two treatments explained variation in δ 34 S leaves , F sulfide , TS leaves , SS leaves , Z. noltei shoot density, Z. noltei leaf, rhizome and root biomass, δ 34 S Loripes , TS Loripes , L. orbiculatus condition, pore water sulfide, total nitrogen and phosphate concentration, L. orbiculatus biomass, L. orbiculatus density, and L. orbiculatus shell length.
Assumptions of normality and homogeneity of residuals were tested for all linear mixed effects models, and if assumptions were violated then non-parametric tests (Friedman and post hoc Conover's test) were applied. P-values for fixed terms in linear mixed effects models were computed through Satterthwaite approximation for denominator degrees of freedom using the "lmerTest" R package (Kuznetsova et al., 2017). This method has Type I error rates closer to 0.05 than other commonly used methods for deriving P-values for fixed terms in linear mixed effects models, especially for smaller sample sizes (Luke, 2017). Differences with P-values <0.05 were considered significant. The R-package "lsmeans" (Lenth, 2016) was used to calculate Least Squares Means (LSM, population means) of each response variable for the four experimental groups, based on the linear mixed effects models. These LSM-values were used to calculate the absolute and relative differences between experimental groups. When the statistical interaction between both treatments was not significant, the main effect of the lucinid addition treatment on each response variable was determined by comparing the mean of the LSM-values for C-and OM-plots with the mean of the LSM values for the L-and L + OM-plots. Likewise, the main effect of the organic enrichment treatment on each response variable was determined by comparing the mean of the LSM-values for Cand L-plots with the mean of the LSM-values for the OMand L + OM-plots.

Field Experiment
Zostera noltei leaf δ 34 S signatures (δ 34 S leaves ) were lower in plots to which organic matter was added (OM and L + OM) compared to plots where no organic matter was added (C and L) (10.5 vs. 12.0 , F 1,20 = 10.6, P = 0.004; Figure 3A). Addition of organic matter resulted in a relative increase in sediment sulfide r rm = 0.81 FIGURE 2 | Positive correlation between seagrass above-ground biomass and Loripes orbiculatus density. Benthic samples collected within the same sampling site are given the same color, with corresponding lines to show the repeated measures correlation fit for each sampling site. For the location of the different sampling sites within Thau lagoon see Figure 1B. DW = dry weight.
No significant statistical interactions were found for any of the response variables. For a complete overview of all Analyses of Variance tables see Supplementary Table S1. For each experimental response variable, the mean and standard error (SE) are presented per experimental group in Table 1. Shell length (mm) 7.2 ± 0.3 9.0 ± 0.2 5.9 ± 0.4 6.7 ± 0.3 C, control; L, addition of Loripes orbiculatus clams; OM, addition of organic matter; L + OM, addition of both L. orbiculatus clams and organic matter. DW, dry weight (n = 6 per treatment).*Only specimens selected that were >7 mm in shell length.

DISCUSSION
In this study we tested whether the facultative mutualism between seagrasses and lucinid bivalves with endosymbiotic sulfide-consuming gill symbionts that was revealed by laboratory experiments (Van Der Heide et al., 2012) can mitigate sulfide stress in seagrasses under field conditions. In addition, we tested whether the dependency of seagrass on this sulfide-detoxification mutualism increased under conditions of enhanced sediment sulfide production. The field survey revealed a strong positive correlation between seagrass above-ground biomass and lucinid (i.e., Loripes orbiculatus) densities (and vice versa) explaining 66% of their respective variation [repeated measures correlation coefficient (r rm ) = 0.81], and pore water sulfide concentrations close to zero at all sites. These results are to be expected when the mutualistic sulfide-detoxification feedback mediates seagrass productivity and lucinid bivalve densities by preventing accumulation of sulfide concentrations in the sediment. However, other mechanisms, such as lower predation risk for lucinids with increasing seagrass biomass and associated refuge provisioning (De Fouw et al., 2016b) may also have contributed to this positive correlation. The field experiment revealed that when organic matter was added to the seagrass sediments, pore water sulfide levels increased considerably and seagrass (i.e., Zostera noltei) suffered from both increased sedimentary sulfide intrusion into their leaves (a proxy for sulfide stress), and detrimental effects on rhizome biomass. This shows that under field conditions the mutualistic sulfide-detoxification feedback was unable to adequately mitigate the sudden dramatic increase in pore water sulfide production as a result of our experimentally enhanced organic loading. Yet, when lucinid (i.e., L. orbiculatus) densities were experimentally enhanced, sedimentary sulfide intrusion into seagrass leaves was significantly reduced, regardless of the addition or organic matter. This result provides the first fieldbased evidence that lucinid bivalves and their sulfide-oxidizing gill symbionts mitigate sulfide stress in seagrasses, and that their sulfide buffering capacity is resilient to conditions of enhanced sediment sulfide production.
Moreover, addition of organic matter had a positive effect on L. orbiculatus body condition, total L. orbiculatus biomass and sulfur content in L. orbiculatus flesh tissues, while experimental enhancement of L. orbiculatus densities had a negative effect on L. orbiculatus condition. This suggests that L. orbiculatus and their gill-symbionts benefited from the higher sulfide production rates caused by the organic enrichment, while some form of negative density dependence occurred when natural L. orbiculatus densities were increased, probably due to competition for sulfide. Negative density dependence was also reported for L. orbiculatus living in a tropical intertidal Zostera noltei system (Banc d'Arguin, Mauritania) (Van Der Geest et al., 2019). These results imply that in our temperate study system, the L. orbiculatus population was close to carrying capacity, and that under current environmental (i.e., sulfide) conditions, the experimentally enhanced L. orbiculatus densities cannot be sustained in the long term. As body condition of L. orbiculatus is known to be positively correlated with reproductive output  and likely also with growth and survival (Paine, 1976), the population size of L. orbiculatus is expected to increase when organic inputs and associated pore water sulfide production gradually increase over time. This could eventually lead to a new equilibrium where the higher sulfide productions rates are again buffered by the L. orbiculatus population that has now increased in size, thus resulting in overall low sediment L. orbiculatus condition expressed as the flesh/shell dry weight (DW) ratio per treatment (n = 6 per treatment) after 50 days. L, addition of L. orbiculatus clams; OM, addition of organic matter. Midline in box; median; box: 25th and 75th percentiles; whiskers: 1.5 × interquartile range; circles: outliers. Note that for each response variable we plotted the mean value per benthic core.
sulfide levels, and low sulfide stress for seagrasses. However, this new equilibrium can only be reached when this organic enrichment does not trigger a breakdown of the mutualistic feedback between seagrass and lucinid bivalves and their sulfideoxidizing gill-bacteria (De Fouw et al., 2016a). Experimental enhancement of lucinid densities resulted in a relative reduction of sedimentary sulfide intrusion in Z. noltei leaves (F sulfide ) of 9%, while addition of organic matter resulted in a relative increase in F sulfide of 21%, but there was no significant interactive effect between the two treatments. This suggest that the sulfide detoxifying capacity of the seagrass-lucinid mutualism was resilient to conditions of suddenly increased sediment sulfide production, but that this capacity did not increase as a function of sediment sulfide production. However, in absolute terms, the mean effect of experimental enhancement of lucinid bivalves on sulfide intrusion was 2.5-fold higher in plots enriched with organic matter [F sulfide from 25.3% (OM plots) to 22.3% (L + OM plots)] compared to plots where no organic matter was added [F sulfide from 20.3 (C plots) to 19.1% (L plots)]. Overall, these results provide the first field-based evidence that lucinid bivalves and their sulfide-oxidizing gill symbionts mitigate sulfide stress in seagrasses, and suggest that the dependency of seagrass on the seagrass-lucinid mutualism will increase under conditions of enhanced sediment sulfide production, as predicted for the near future.

Sulfide Intrusion Under Natural Conditions and During Organic Loading
In our field experiment, sulfide levels were not manipulated by addition of sulfide, but rather by addition of organic matter (starch mixed with sucrose). Our organic matter treatment successfully stimulated sulfide production in the sediment by sulfate-reducing bacteria, resulting in mean pore water sulfide concentrations of approximately 1200 and 700 µmol L −1 at the end of the experimental period in OM and OM + L plots respectively, compared to pore water sulfide concentrations close to zero in control plots (C and L) ( Figure 6B and Table 1). This suggests that the observed increase in F sulfide and TS leaves in plots to which organic matter was added (OM and L + OM) compared to their respective controls (C and L) (Figures 3B,C) was caused by the sulfidic environment created by the addition of organic matter. Thus, our hypothesis that sedimentary sulfide intrusion increases in plots with experimentally enhanced sulfide production rates (induced by organic enrichment) was confirmed. This implies that in plots enriched with organic matter, Z. noltei was not able to prevent sedimentary sulfide from entering its vulnerable meristematic leaf tissues, a phenomenon that is considered to be an indicator of sulfide stress, with negative impact on photosynthesis, growth and seagrass survival (reviewed by Holmer and Hasler-Sheetal, 2014;Kilminster et al., 2014). In our field experiment, sulfide stress was reflected by a reduction in Z. noltei rhizome biomass in plots to which organic matter was added (Figure 4 and Supplementary Table S1). However, none of the other seagrass fitness measures (i.e., Z. noltei shoot density, leaf biomass, and root biomass) were affected by the increased sulfide stress that was generated by our organic matter addition treatment. These results are somewhat surprising, since the mean pore water sulfide levels of almost 1200 µmol L −1 that were measured in plots to which only organic matter was added (i.e., OM-plots) were well above the previously reported threshold values at which sulfide toxicity effects were observed in Z. noltei, i.e., (i) decreased above-ground production at sulfide levels >200 µmol L −1 (Van Der Heide et al., 2012), (ii) decreased leaf elongation rate at sulfide levels >500 µmol L −1 (Lamers et al., 2013), and (iii) reduced patch expansion at sulfide levels >1000 µmol L −1 . This may suggest that the duration of the experiment (50 days) was too short for our organic enrichment-treatment to have a significant effect on these other seagrass fitness measures through enhanced sulfide intrusion. Alternatively, under natural conditions, sulfide stress may have been relatively low in our study system, which may have resulted in Z. noltei being more resilient to the sulfide stress induced by our organic matter treatment.
Whereas our benthic survey revealed that only seagrass aboveground biomass was positively correlated with lucinid densities, our field experiment revealed that only seagrass below-ground biomass (i.e., rhizomes) was negatively affected by experimentally enhanced sulfide stress (induced by addition of organic matter). These results suggest that the observed lack of a positive correlation between seagrass below-ground biomass and lucinid densities under natural conditions, could potentially be attributed to ambient sulfide levels being too low to be toxic to seagrass below-ground tissues.
The fact that addition of organic matter only affected rhizome biomass could be explained by the rhizomes acting as a primary buffer in detoxifying toxic sedimentary sulfide. Part of the sedimentary sulfide entering the below-ground tissues through passive diffusion will then be assimilated into organic sulfur compounds, and a part will be oxidized to elemental sulfur (S 0 ) and sulfate in the rhizomes (Hasler-Sheetal and Holmer, 2015;Holmer et al., 2017). However, when pore water sulfide levels become too high, the rhizomes will not be able to neutralize all of the toxic sulfide entering the below-ground tissues, resulting in rhizome tissues being exposed to sulfide, which may negatively affect rhizome biomass. If true, this indicates the presence of a sulfide detoxification mechanism in Z. noltei as was previously described for Z. marina (Hasler-Sheetal and Holmer, 2015), where underground tissues possessed the highest detoxification capacity (86%), with the rhizome (61%) acting as the primary buffer in detoxifying toxic sediment sulfide and protecting the vulnerable meristematic tissue in the leaves. Yet, more studies are needed to confirm the notion that the rhizomes of Z. noltei are acting as a primary buffer in detoxifying toxic sediment sulfide.

Organic Matter-Induced Effects on Stable Sulfur Isotope Signals
Addition of organic matter resulted in lower δ 34 S signatures of Z. noltei leaf tissues (δ 34 S leaves ) ( Figure 3A). The two potential sources of sulfur available for uptake by seagrass (i.e., sulfate and sulfide), have very distinct δ 34 S signatures. Average δ 34 S values for seawater sulfate are + 20.7 ± 0.6% (mean ± SE), while those for sulfide are −18.9 ± 8.7 for the AVS fraction and −22.4 ± 8.4 for the CRS fraction in the sediment (reviewed by Holmer and Hasler-Sheetal, 2014). Therefore, the lower δ 34 S leaves values in plots to which organic matter was added can be attributed to increased sedimentary sulfide intrusion into the Z. noltei leaf tissues ( Figure 3B).
In contrast, addition of organic matter resulted in higher δ 34 S signatures of L. orbiculatus flesh tissues (δ 34 S Loripes ) ( Figure 5A). This may be explained by enhanced sulfate reduction rates due to organic enrichment, which may have caused decreased discrimination of 34 S by the sulfate-reducing bacteria, leading to higher δ 34 S values for pore water sulfide (Hasler-Sheetal and Holmer, 2015). Subsequently, this 34 S-enriched pore water sulfide may have been used by the sulfide-oxidizing gill-symbionts of L. orbiculatus and incorporated into the lucinid tissues, which may have resulted in the observed enrichment of δ 34 S Loripes values in plots to which organic matter was added.

Sulfide-Driven Redistribution of Lucinid Bivalves
Addition of organic matter had an unforeseen positive effect on L. orbiculatus total biomass ( Figure 6A and Supplementary Table S1), which was mainly driven by higher L. orbiculatus densities of clams with smaller shell length (see Supplementary  Table S1 and Supplementary Figures S1A,B) in plots to which only organic matter was added. This indicates that during our experiment, especially small-sized (∼5 mm in shell length) L. orbiculatus clams immigrated into the OM-plot area, presumably to profit from the higher sulfide concentrations in these plots compared to the surroundings. This notion is supported by the observed positive effect of our OM-treatment on L. orbiculatus condition and L. orbiculatus tissue sulfur content, and the observed negative effect of our L-treatment on L. orbiculatus condition (Figures 5B,C and Supplementary Table  S1). This result also points toward a sulfide-limited L. orbiculatus population that is likely to increase its sulfide buffering capacity as a function of sulfide production.
As a result of this presumed sulfide-driven redistribution of L. orbiculatus during the experimental period, the observed effects of our OM-treatment on sulfide intrusion (and subsequent seagrass performance) are most likely conservative. Moreover, the notion that (small) L. orbiculatus clams may be able to redistribute themselves depending on sulfide availability suggests a locally adaptive sulfide-detoxification mutualism that can swiftly neutralize local pulses of sulfide caused by deposition of organic matter (e.g., algal mats) and eutrophication. However, more research is needed to support this claim.

Factors Influencing the Strength of the Seagrass-Lucinid Mutualism
Our benthic field survey revealed a positive correlation between seagrass above-ground biomass and L. orbiculatus densities explaining 66% of their respective variation. Nevertheless, it appears surprising that a similar relationship was not found with seagrass below-ground biomass. This result may suggest that the strength of the tripartite mutualism critically depends on the direct relationship between seagrass above-ground biomass and photosynthesis-driven oxygen release by the roots, and the indirect relationship between seagrass above-ground biomass and organic matter trapping, which results in higher sediment sulfide production rates. Yet, more empirical studies are needed in various seagrass-lucinid systems to support the generality of this correlation.
As anaerobic decomposition of organic matter and related sulfide production are strongly temperature dependent, it has been suggested that the vigor of the seagrass-lucinid mutualism will also depend on temperature. This notion is supported by a meta-analysis showing a lower proportion of studies that found seagrass-lucinid associations in temperate (56%) and subtropical seagrass beds (90%) compared to tropical meadows (97%) (Van Der Heide et al., 2012). However, in our subtidal temperate Z. noltei-dominated study system 66% of the variation in L. orbiculatus densities was explained by seagrass (above-ground) biomass (and vice versa), while this was only 42% for L. orbiculatus densities living in the tropical intertidal Z. noltei-dominated meadows of Banc d'Arguin, Mauritania (Van Der Heide et al., 2012). Moreover, in our subtidal temperate Z. noltei-dominated study system the mean density of L. orbiculatus was 1771 individual clams per m 2 , which is much higher than mean densities reported for the tropical intertidal Z. noltei beds of Banc d'Arguin, being 320-840 individual clams per m 2 (Honkoop et al., 2008;Ahmedou Salem et al., 2014). Note that mean shell length of L. orbiculatus collected during our benthic survey [6.9 ± 2.7 (mean ± SD) mm; range 2.4-16.5; n = 292] was very similar to those collected by Ahmedou Salem et al. (2014) for tropical Banc d'Arguin (7.4 ± 2.2 mm; range 1.7-12.5; n = 2169) (unpublished data).
This illustrates that other factors, such as light availability and sediment conditions (e.g., organic matter content and lability, bioturbation, grain size, iron availability) that affect production, oxidation and precipitation of sulfide in the rhizosphere, and species-specific sulfide tolerance, may also have a strong effect on the mutualism dependency. Thus, to get a better mechanistic understanding of the underlying factors that determine the strength of the sulfide-detoxifying mutualism, latitudinal studies are needed covering coastal systems with different organic loading regimes and varying sediment and light conditions, for a variety of seagrass species.

Ecological Considerations
Our results provide the first field-based evidence that the recently discovered tripartite mutualism between seagrasses, lucinid bivalves and their sulfide-consuming gill symbionts can mitigate sulfide stress in seagrasses. Moreover, our results suggest that the dependency of seagrass on this sulfide-detoxification mutualism will increase when sediment sulfide production rates increase, which is expected due to ongoing global warming and organic loading of coastal systems worldwide. These findings show that the seagrass-lucinid mutualism, when (historically) present, should be taken into account for seagrass conservation and restoration and that specific management should be developed to protect, enhance and use this mutualistic network for improving long-term restoration success and seagrass resilience to global change. Yet, previous studies also indicate that environmental disruption of this positive feedback can trigger mutualistic breakdown, potentially leading to seagrass ecosystem regime shifts (De Fouw et al., 2016a, 2018. Coastal zone management should be aware that it may be very difficult and costly to reverse such a regime shift and re-establish the seagrass community (Scheffer et al., 2001;Le Fur et al., 2019). We therefore recommend that future research focus on exploring the threshold environmental conditions and rates of environmental change under which the mutualism can still function, knowledge that is key to preserve the sulfide-detoxification mutualism and may contribute to the development of a safe operating space for seagrass ecosystems.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available from the Dryad Digital Repository doi: 10.5061/dryad.jm63xsj6n.

AUTHOR CONTRIBUTIONS
All authors designed the study and interpreted the data. MG and RW collected the data. MG analyzed the data and drafted the manuscript. TH, MH, and RW contributed to critical revision of the submitted version.

ACKNOWLEDGMENTS
Sonia Cherkaoui, Célia Graziani, Clémence Laborde, Grégoire Thiry, and Mariam Maki Sy assisted in collecting and processing benthic samples. Katrine Clement Kirkegaard is thanked for doing isotope analyses. Jaap van der Meer gave constructive comments on the statistics, and Tamar Lok and two reviewers are thanked for their constructive comments on a draft of the manuscript. We also thank Nathalie Barré at the Pôle Relais Lagunes Méditerranéennes for discussing and promoting the possible implications of our findings for the management of Mediterranean coastal lagoons.