2.1. Study sites and sampling

The study sites were located in macrotidal mudflats in Brittany (north-western France) under eutrophic conditions where green algae mats are observed (https://bretagne-environnement.fr/donnees-algues-vertes-bretagne). The mudflats that we selected for our study are the mudflats that are regularly monitored as part of the green tide monitoring programs set up by the Loire-Bretagne Water Agency (CEVA], final report, 2015]). Over the sampling period (between mid-April to mid-June 2019), the tidal range fluctuated from 2 to 9 m. Overall, 200 sediment samples were collected from 45 sites (Figure 1, Table S1) at mid–low tide. Sediment cores were sampled with a PVC core (diameter = 6 cm,h = 20 cm) in the upper 10 cm sediment layer for the benthic flux measurements, and another core (diameter = 9 cm, h = 5 cm) was sampled in order to analyze the surface sediment characteristics (Figure S1). For the benthic flux measurements, the cores were incubated in the dark directly on site in a mobile laboratory under controlled temperature (19 ± 2 °C) within one hour following sampling (see below). Based on the sample collected using the PVC core with a height of 5 cm, an aliquot of the wet sediment, with a known volume and weight, was maintained at 4 °C and used to determine the porosity. The remaining sediment was frozen at −20 °C in order to analyze the pigment content (chla and phaeopigment), phosphorus speciation, elemental composition and grain-size.

2.2. Benthic fluxes

The cores were incubated in the dark during 4 hours for the ${\text{NH}}_4^{+}$ and PO₄ efflux measurements. The overlying water was replaced by 150 mL of nutrient-free artificial seawater ([NaCl] = 33 g⋅L⁻¹, [NaHCO₃] = 0.2 g⋅L⁻¹, pH ≈ 8) and gently aerated by bubbling in order to preserve the redox conditions of the sediment. Two 0.2 μm-filtered water samples were collected in the overlying water after 2 h and 4 h of incubation and stored at refrigerator temperature (4 °C) for less than three days prior to the nutrient analysis. The ${\text{NH}}_4^{+}$, PO₄, ${\text{NO}}_2^{-}$ and ${\text{NO}}_3^{-}$ concentrations in the overlying water were measured with the colorimetric method using an automated photometric analyzer Gallery™, with a detection limit of 0.9, 0.2, 0.07 and 3.6 μM, respectively. All of the ${\text{NO}}_2^{-}$ and ${\text{NO}}_3^{-}$ concentration measurements were below the detection limit.

The ${\text{NH}}_4^{+}$ and PO₄ fluxes across the sediment–water interface were accessed by using the change in the molar concentration of the solute in the known volume of overlying water as a function of incubation time and the surface area of the sediment core [Aller et al. 1985]. If the rate of nutrient release from the sediment did not follow a linear trend over the incubation period, only the first sample was considered in the flux estimation. This was the case for PO₄ when the exchanges between the sediment and water reached an equilibrium state after two hours of incubation, in general, due to the re-adsorption onto particles [Sundby et al. 1992]. Note that these fluxes were conducted under standardized conditions (temperature, aeration) and therefore considered as potential fluxes.

All of the values are reported in Supplementary Table S1.

2.3. Analysis of the sedimentary characteristics

The porosity was calculated by drying a previously weighed wet aliquot sediment at 60 °C. The water loss, determined by mass difference, and the sediment density set at 2.55 were used to calculate the porosity [Krom and Berner 1980]. The particle-size distribution (<2 mm) was measured using a laser diffraction instrument (Malvern Mastersizer). The particles were classified as either sand (63 −−2000 μm), silt (3 −−63 μm) or clay (<3 μm) fractions [e.g. Keil and Hedges 1993; Pye and Blott 2004]. “Mud” is defined as the sum of the clay and silt particles. The percentage of mud was used as a sedimentary parameter thereafter.

The pigment content (chla and phaeopigment) was measured in freeze-dried sediment (0.5–1 g) extracted in 90% acetone (10 mL) in the dark (18–20 h) at 4 °C. Each sample was previously gently ground using an agate pestle and mortar. After centrifugation, the chla and phaeopigment contents were measured in the supernatant using the spectrophotometric method of Lorenzen [1967] at 665 and 750 nm (Uvikon spectrophotometer), and expressed in μg⋅g⁻¹ of dry sediment. The detection limits of the chla and phaeopigment contents were 2.7 and 12.9 μg⋅g⁻¹, respectively.

The total organic carbon and nitrogen (TOC and TN) contents were determined using an element analyzer (FLASH™ 2000 OEA). An aliquot of freeze-dried and crushed sediment was acid-treated with 2N HCl to remove the carbonate and was subsequently rinsed with deionized water. After centrifugation, the carbonate-free sample was dried at 60 °C, and ground before being placed into a tin capsule for the TOC analysis. A second aliquot without an acidification treatment was used to determine the TN analysis.

The iron oxide-bound P (Fe-P) content was determined using a Dithionite-Bicarbonate solution [Ruttenberg 1992], as described in Andrieux-Loyer et al. [2008]. The total P (P_tot) content was determined using a 1 mol⋅L⁻¹ HCl treatment overnight after sediment ignition at 550 °C (4 h), while the inorganic P content refers to the sum of the P-forms (Fe-bound P, Ca-bound P and detrital P), that was extracted with 1 mol⋅L⁻¹ HCl before sediment ignition [Aspila et al. 1976]. The organic P (Org-P) content was then quantified by calculating the difference between the total P and inorganic P contents [Andrieux-Loyer et al. 2008]. The phosphorus-form extracts were subsequently analyzed using segmented flow analysis (SFA) [Aminot and Kérouel 2007]. The Org-P and Fe-P represented the pool of potentially bioavailable P.

The TOC, TN, P_tot, Org-P and Fe-P contents were expressed as the mass of the carbon, nitrogen and phosphorus in the total dry mass of the sediment.

The C:N and TN:Org-P ratios (mol:mol) in SOM were calculated from the TOC, TN and Org-P contents.

All of the values are reported in Supplementary Table S1.

2.4. Statistical analysis

Pearson’s correlation matrix was calculated to establish the pairwise correlations between the sedimentary characteristics of the sediment. The relationships between the benthic nutrient fluxes and the sedimentary characteristics were assessed through multiple linear regressions. To better predict ${\text{NH}}_4^{+}$ and PO₄ fluxes, two models were built by multiple linear regression (MLR) using a stepwise selection procedure based on the Akaike information criterion (AIC). To ensure that the multi-collinearity did not skew our results, the variance inflation factors (VIF) were measured. None of the VIF values were higher than 3 (a threshold value was set at 5). In addition, when the selected variables had a Pearson’s correlation higher than 0.7 between them, we checked that there was no interaction effect in the model. To identify the main drivers of the benthic nutrient fluxes by MLR, the data were normalized beforehand, and a p-value < 0.05 for the F-test was used to consider a significant effect of one predictor variable on the model. In addition to the calculation of AIC, the K-fold cross validation was run to evaluate the performance of each model built by linear regressions. The data set was split into 10 folds. In order to compare two means, a Student’s t-test was performed when the sample size (n) was higher than 30. Each statistical analysis was performed using the R-studio software. The “aod”, “car” and “caret” packages were used to calculate AIC and VIF and to perform the k-fold cross validation, respectively. The “lm”, “step” and “scale” functions were used to perform the linear regressions, the stepwise selection procedure by MLR and the database normalization, respectively.

3.1. Benthic nutrient fluxes

Taking the entire data set into account, the average ${\text{NH}}_4^{+}$ and PO₄ fluxes were 101 ± 117 μmol⋅m⁻²⋅h⁻¹ and 17 ± 20 μmol⋅m⁻²⋅h⁻¹, respectively. Over 50% of the ${\text{NH}}_4^{+}$ fluxes ranged from 19 and 151 μmol⋅m⁻²⋅h⁻¹, with values up to 526 μmol⋅m⁻²⋅h⁻¹ (Figure 2A). With regards to the PO₄ fluxes, 50% of the data ranged from 3 and 24 μmol⋅m⁻²⋅h⁻¹, with a maximum flux of 172 μmol⋅m⁻²⋅h⁻¹ (Figure 2B). In general, the PO₄ fluxes were lower than those of ${\text{NH}}_4^{+}$. Over 50% of the N:P flux ratios ranged from 2 to 12, with values up to 765. The large variability in the ratios of the N:P fluxes reflects a weak correlation between the ${\text{NH}}_4^{+}$ and PO₄ fluxes (Pearson’s coefficient = 0.31, p-value < 0.05) (Figure 2C). The lowest N:P flux ratios (N:P ⩽ 2) were significantly explained by high PO₄ fluxes (mean = 25 μmol⋅m⁻²⋅h⁻¹) and low ${\text{NH}}_4^{+}$ fluxes (31 μmol⋅m⁻²⋅h⁻¹) (Student’s t-test for one sample, p-value < 0.05). The highest ratios for the N:P fluxes (N:P ⩾ 12) were significantly explained by low PO₄ fluxes (mean = 7 μmol⋅m⁻²⋅h⁻¹) and high ${\text{NH}}_4^{+}$ fluxes (192 μmol⋅m⁻²⋅h⁻¹) (Student’s t-test for one sample, p-value < 0.05) (Figure S2).

The variability of both ${\text{NH}}_4^{+}$ and PO₄ fluxes at the regional scale allowed to distinguish the locations where the sediment presented the highest average effluxes of nutrient. For the ${\text{NH}}_4^{+}$ effluxes, they were Gulf of Morbihan (mean: 185 μmol⋅m⁻²⋅h⁻¹), Vannes Estuary (mean: 173 μmol⋅m⁻²⋅h⁻¹) and Auray River (mean: 162 μmol⋅m⁻²⋅h⁻¹) (Figure 3). For the PO₄ effluxes, they were Lorient Bay (mean: 34 μmol⋅m⁻²⋅h⁻¹), Pont L’Abbé (mean: 34 μmol⋅m⁻²⋅h⁻¹) and Rance Estuary (mean: 30 μmol⋅m⁻²⋅h⁻¹) (Figure 4). By contrast, the three locations which presented the lowest average effluxes of ${\text{NH}}_4^{+}$ were Aber Wrac’h (mean: 19 μmol⋅m⁻²⋅h⁻¹), Morlaix Bay (mean: 48 μmol⋅m⁻²⋅h⁻¹) and Port La Forêt (mean: 79 μmol⋅m⁻²⋅h⁻¹) (Figure 3), and those with the lowest average effluxes of ${\text{PO}}_4^{3-}$ were Ria Etel (mean: 3 μmol⋅m⁻²⋅h⁻¹), Gulf of Morbihan (mean: 4 μmol⋅m⁻²⋅h⁻¹) and Port la Forêt (mean: 7 μmol⋅m⁻²⋅h⁻¹) (Figure 4). A high spatial variability of these fluxes was also observed between the sampling sites of the same location (Figures 3 and 4, Table S1). For example, in the Gulf of Morbihan, the lowest effluxes of ${\text{NH}}_4^{+}$ were measured at the site #45 (7.3 ± 1.4 μmol⋅m⁻²⋅h⁻¹), whereas the site #44 presented large effluxes of ${\text{NH}}_4^{+}$ (388.6 ± 189.5 μmol⋅m⁻²⋅h⁻¹) (Figure 3).

3.2. Surface sediment characteristics

The distribution and pairwise correlations with Pearson’s coefficient for the surface sediment characteristics measured in the present study are presented in Tables 1 and 2.

The particle-size distribution of the samples indicates that the median grain-size (D50) ranged from 16 to 519 μm. With an average mud content of 68.3 ± 17.2% for all of the samples, the clay and silt particles represented 4.7 ± 2.1 and 63.6 ± 15.4% of the particle-size distribution, respectively. According to the ternary diagram based on the sand/mud ratios proposed by Flemming [2000], a large proportion of the surface sediments collected in the present study were classified as sandy mud (38%) and slightly sandy mud (49%). As expected, the mud content was positively correlated with the porosity of the sediment (Pearson’s coefficient = 0.82, p-value < 0.05). The average porosity was 69.8 ± 10.3% with values that ranged from 46 to 87%.

The mud content was also positively correlated with the TOC (Pearson’s coefficient = 0.64, p-value < 0.05), TN (Pearson’s coefficient = 0.73, p-value < 0.05) and Org-P contents (Pearson’s coefficient = 0.72, p-value < 0.05) in the surface sediment. The TOC, TN and Org-P contents averaged 2.2 ± 1.2%, 0.24 ± 0.12% and 0.026 ± 0.014%, respectively. With an average P_tot content of 0.055 ± 0.022%, the amount of Org-P was 45 ± 12% of the P_tot content compared to 14 ± 6% for the Fe-P (average content of 0.008 ± 0.005%).

The C:N and TN:Org-P ratios in the SOM averaged 10.0 ± 2.4 and 23.1 ± 6.9, respectively. The C:N ratios reached up to 21.9 with half of the values ranging from 8.7 to 10.9. The highest values corresponded to the sediment samples collected in the Goulven Bay (mean = 17.8 ± 3.9).

The accumulation of OM from the algal biomass (which included microphytobenthos, pelagic phytoplankton, macroalgae) in the surface sediment was described by the chla and phaeopigment contents, with a mean value of 5.8 ± 4.8 and 28 ± 18.7μg⋅g⁻¹ of dry sediment, respectively. The chla and phaeopigment contents were positively correlated with the TOC, TN and Org-P contents (Pearson’s coefficient > 0.6, p-value < 0.05).

3.3. Relationships between the sedimentary parameters and benthic nutrient fluxes

The stepwise multiple linear regression (MLR) allowed to select the best combination of sedimentary characteristics to predict the ${\text{NH}}_4^{+}$ and PO₄ fluxes. With regards to the PO₄ flux, the parameters selected by the final MLR model were the phaeopigment, Fe-P and Org-P contents, and the porosity (Figure 5). Whereas the PO₄ flux was positively and significantly correlated with the phaeopigment and Fe-P contents (p-value < 0.05), it was negatively and significantly correlated with the porosity (p-value < 0.05). In addition, a positive relationship was observed between the Org-P content and the PO₄ flux, despite a less significant effect (p-value = 0.079). To predict the ${\text{NH}}_4^{+}$ flux, the porosity, chla content and C:N and TN:Org-P ratios in the SOM were selected by the final MLR model (Figure 6). The porosity, as well as the TN:Org-P ratio and chla content, were positively correlated with the ${\text{NH}}_4^{+}$ flux. In contrast, a negative relationship was observed between the C:N ratio and the ${\text{NH}}_4^{+}$ flux. The parameter with the highest significant effect was porosity (p-value < 0.001), followed by the TN:Org-P ratio (p-value = 0.01). The chla content and C:N ratio presented a less significant relationship with the ${\text{NH}}_4^{+}$ flux (p-value > 0.05).

The accuracy of the prediction remained low for both models, with only 18 and 24% of the ${\text{NH}}_4^{+}$ and PO₄ flux variations explained, respectively.

4.1. What are the main parameters of surface sediment driving the PO₄ fluxes?

Four parameters were selected via the final MLR model in order to predict the PO₄ fluxes: porosity, phaeopigment content, Fe-P content and Org-P content. Among the four selected parameters, the correlation between the Org-P content and the PO₄ flux was not significant (F-test p-value > 0.05). Therefore, the porosity, phaeopigment and Fe-P contents were considered as the key variables for predicting PO₄ fluxes.

The parameter presenting the most significant positive correlation for the PO₄ fluxes was the phaeopigment content (F-test p-value = 0.002). The phaeopigment content is an indicator of algal biomass detritus in the SOM composition [Dell’Anno et al. 2002; Therkildsen and Lomstein 1993]. In the present study, the significant and positive correlations between the phaeopigment content and that of the TOC, TN and Org-P contents indicated that the organic matter was enriched in the surface sediment due to the accumulation of algal detritus. In the Auray River (Brittany, France), Andrieux-Loyer et al. [2014] also observed that the increase in the phaeopigment, TOC, TN and Org-P contents in the surface sediment was fueled by the occurrence of green macroalgae. An accumulation of algal OM most likely results in intense SOM mineralization, leading to both oxygen depletion in porewater and a shift in anaerobic mineralization from iron to sulfate reduction thereby enabling iron-sulfide formation [Lehtoranta et al. 2009]. Under these anoxic conditions, the adsorption of PO₄ onto iron oxides would be less efficient, increasing the upward diffusion of PO₄ from the sediment [Ekholm and Lehtoranta 2012 and references therein; Lehtoranta et al. 2009; Rozan et al. 2002]. Thus, surface sediment that is rich in algal detritus would be a more effective source of PO₄ due to the creation of anoxic conditions and the subsequent desorption of P (Figure 7A). The solubilization of Fe-P has been shown to be a main process in PO₄ release in overlying water from sulfidic sediments during macroalgal blooms [Rozan et al. 2002]. Consequently, sediment that is rich in Fe-P would represent a potentially bioavailable source of P under anoxic conditions. This is in good agreement with our results which identify the Fe-P content as the second main driver, and which was positively correlated with the PO₄ flux.

In contrast to the previous variables, porosity, the third main driver, was negatively correlated with the PO₄ flux. An increase in porosity favors the diffusion of solutes through the sediment [Boudreau 1996], and as a result, a positive effect on the PO₄ fluxes could be expected due to an upward diffusion from the sediment. On the other hand, high porosity may also improve the oxygen diffusion across the sediment–water interface [House 2003]. The ferrous ion (Fe²⁺) oxidation followed by the re-precipitation of iron oxides at the oxygenated interface, would increase the adsorption sites for PO₄ [e.g. Krom and Berner 1980; Gunnars et al. 2002; Mayer and Jarrell 2000; Zhang et al. 2010]. This may limit the release of PO₄ into the overlying water, thereby explaining the negative correlation between the porosity and the PO₄ flux observed in the present study.

4.2. What are the main parameters of surface sediment driving the ${\text{NH}}_4^{+}$ fluxes?

The parameters selected by the final MLR model for predicting the ${\text{NH}}_4^{+}$ fluxes were the porosity, the TN:Org-P ratio, the chla content and the C:N ratio. Only the porosity and TN:Org-P ratio had a significant effect (F-test p-value < 0.05), and were thus considered as the key variables for predicting ${\text{NH}}_4^{+}$ fluxes.

Contrary to the PO₄ flux, the ${\text{NH}}_4^{+}$ flux was positively correlated with the porosity. When sediment has a high porosity, this likely enhances the diffusion of solutes across the sediment–water interface [Boudreau 1996] including that of ${\text{NH}}_4^{+}$. As discussed above, an increase in porosity could also enhance the oxygen diffusion at the sediment–water interface. This may favor the oxidation of ${\text{NH}}_4^{+}$ by nitrifying microorganisms, and could therefore increase the effluxes of oxidized nitrogen species (NO_x) at the expense of ${\text{NH}}_4^{+}$ [Capone et al. 2008 and references therein]. In the present study, the impact of nitrification on the ${\text{NH}}_4^{+}$ fluxes seemed to be negligible as there was neither a negative effect of the porosity on the ${\text{NH}}_4^{+}$ fluxes nor a release of NO_x in the overlying water ([NO_x] < detection limit) (Figure 7B). Consequently, the enhanced oxygen diffusion in the upper layer of the sediment due to the increase in porosity would poorly impact the efflux of ${\text{NH}}_4^{+}$ in contrast to that of PO₄. This could be explained by the fact that abiotic Fe²⁺ oxidation, fueling the precipitation of Fe-P, would be faster and less oxygen-limited compared to ${\text{NH}}_4^{+}$ oxidation [Canavan et al. 2006].

The positive correlation observed between the porosity and the ${\text{NH}}_4^{+}$ flux might also be explained by the fact that the sediments with higher porosity would be also associated with higher OM input flux, through the positive and significant correlation observed between the porosity and e.g. the TOC content (Table 2). Indeed, an increase in the deposition of OM in the sediment, supported by the river export of fine mud particles (with higher porosity, Table 2), can improve the mineralization rates of SOM, as observed in the Elorn and Aulne estuaries in Brittany [Khalil et al. 2018], and thus lead to a larger ${\text{NH}}_4^{+}$ regeneration in porewaters of sediment. Thereby, higher ${\text{NH}}_4^{+}$ recycling rates combined with higher sediment–water exchanges could enhance the diffusion of ${\text{NH}}_4^{+}$ to the overlying water.

In addition to the porosity, the ${\text{NH}}_4^{+}$ flux was significantly and positively correlated with the TN:Org-P ratio in the SOM. This suggests that the ${\text{NH}}_4^{+}$ fluxes were controlled by the elemental composition of the SOM. SOM that is rich in both N and easily biodegradable compounds enhances the production of ${\text{NH}}_4^{+}$ in the sediment via SOM mineralization [Capone et al. 2008 and references therein], releasing ${\text{NH}}_4^{+}$ into the overlying water (Figure 7B). In coastal systems, the SOM composition depends on a large variability of marine and terrigenous sources of OM [e.g. Carlier et al. 2008; Cook et al. 2004; Dubois et al. 2012]. As a result, the variability of the TN:Org-P ratio in the SOM observed in the present study may be controlled by different OM sources. An increase in the TN:Org-P ratio in the SOM ratio could indicate an accumulation of macroalgae in coastal systems. Such a relationship has been observed in the Venetian lagoon during the spring and summer [Sfriso et al. 1988]. An increase in this ratio could also reflect an anthropogenic N input from riverine particulate matter. This was based on imbalances in the N:P ratio observed in watersheds in previous studies, which depend on industrial and agricultural activities, and the management policy [Arbuckle and Downing 2001; Bouwman et al. 2017; Carpenter et al. 1998; Guenther et al. 2015; Sardans et al. 2012].

4.3. Sedimentary proxies: a good tool to predict ${\text{NH}}_4^{+}$ and PO₄ fluxes?

We hypothesized that the chla content and the C:N ratio, i.e. tracers of labile OM input fluxes explain the spatial variability in benthic nutrient fluxes and therefore allowing them to be used as proxies. Such correlations between these two parameters and the ${\text{NH}}_4^{+}$ and PO₄ fluxes were shown by Cowan and Boynton [1996] and Clavero et al. [2000]. In the present study, the chla content and the C:N ratio were selected by the final MLR model of the ${\text{NH}}_4^{+}$ flux, but without significant correlations (F-test p-value > 0.05). With regard to the PO₄ fluxes, the modeling approach did not identify the chla content and the C:N ratio as main drivers. However, other parameters related to the SOM composition presented significant correlations with the benthic nutrient fluxes: i.e. the TN:Org-P ratio for the ${\text{NH}}_4^{+}$ flux, and the phaeopigment content, a chla breakdown product, for the PO₄ flux. The characterization of SOM composition, largely dependent on the spatial variability of OM sources, may require more specific tools (e.g. isotopic or biomolecular markers) to better relate the benthic nutrient fluxes to the sedimentary composition.

In the present study, the stepwise MLR allowed to select the best combinations of parameters to predict ${\text{NH}}_4^{+}$ and PO₄ fluxes. Nevertheless, only 18 and 24% of ${\text{NH}}_4^{+}$ and PO₄ flux variations respectively were explained by our modeling approach. This highlights the difficulty to determine the controls on the benthic nutrient fluxes by sedimentary proxies in the intertidal mudflats at the regional scale where the rate deposition of particulate matter is influenced by river discharges, tidal fluctuations and anthropogenic pressures [Bauer et al. 2013].

In addition, the contribution of biological parameters was not considered in the present work, which could explain the weak accuracy of prediction for both models of ${\text{NH}}_4^{+}$ and PO₄ fluxes. The ${\text{NH}}_4^{+}$ and PO₄ fluxes may be impacted by bioturbation, through particle reworking and burrow ventilation by the benthic macrofauna [Kristensen et al. 2012]. Bioturbation enhances the solute exchanges between the porewater and the overlying water, affects the remobilization of burial OM, and modifies the redox conditions [Graf and Rosenberg 1997; Welsh 2003; Kristensen et al. 2012]. The stimulation of ${\text{NH}}_4^{+}$ and PO₄ fluxes varies according to the density, the burrowing depth, the ventilation of the benthic macrofauna, and the specific-site conditions (e.g. sediment components, nutrient concentrations in porewater). The study carried out in four European estuaries by Nizzoli et al. [2007] has shown that the early stage of sediment colonization by the polychaete Nereis spp. stimulated the ${\text{NH}}_4^{+}$ effluxes in all geographical areas (from 1.5 to 4-fold higher) and promoted the sediment to act either as a sink or a source of PO₄ according to the site-specific conditions. This thus leads to shifts in the magnitude of effluxes and ratios of N:P fluxes which depend on the spatial variability of environmental conditions.