Different paths, same destination: similar benthic nitrogen cycling sustained by distinct biogeochemical mechanisms in artificialized coastal systems

Human activities are increasingly affecting coastal ecosystems, with shoreline stabilization structures becoming a prevalent response to sea-level rise and extreme climatic events. While these structures aim to protect coastal communities and infrastructures, their effects on ecosystem functioning, and particularly nitrogen cycling and its removal processes such as denitrification, remain poorly understood. To assess the ecological impacts of breakwater construction, we employed a space-for-time substitution approach to examine changes in sediment biogeochemistry, macrobenthic community structure, and nitrogen cycling in an intertidal salt marsh in the St. Lawrence Estuary, Canada. We measured benthic fluxes, denitrification rates, and macrofaunal assemblages at five locations: two landward of the breakwater (impacted sites), one situated immediately seaward of the breakwater (intermediate site) and two reference sites (one vegetated and one unvegetated outside the breakwater system). Landward sites exhibited finer sediment with higher organic carbon content compared to other sites and supported distinct macrobenthic communities dominated by opportunistic Oligochaetes, while reference sites were dominated by Molluscs and Crustaceans. The intermediate seaward site consistently resembled the unvegetated reference site across all variables measured at both sites, including sediment characteristics, oxygen dynamics and benthic fluxes, thereby showing reference-like conditions. Surprisingly, despite these substantial physical and biological changes, dark total benthic oxygen uptake, nitrate and ammonium benthic fluxes showed no consistent separation among impacted, intermediate, and reference sites, while benthic denitrification rates were comparable among the impacted and reference sites where they were measured. This functional similarity reflects an apparent functional convergence arising from distinct combinations of abiotic and biotic factors, including root-mediated rhizosphere oxygenation, macrofauna-mediated sediment irrigation, and contrasting reliance on internally produced versus water-column nitrate. However, this convergence might not apply to all ecosystem services offered by natural salt marshes and may depend on specific biological communities, emphasizing the importance of careful consideration in coastal management decisions.

1 Introduction

Coastal areas provide essential ecological, economic, and social services. They host various ecosystems and support fisheries, aquaculture, recreational and cultural activities. Around 10 % of the global population lives in coastal areas less than 10 m above sea-level, with coastal development rising in recent decades, a trend expected to continue (Neumann et al., 2015). However, these areas face significant threats from sea-level rise and shoreline erosion (Mentaschi et al., 2018). For instance, a global temperature rise of 2 °C above pre-industrial temperatures could lead to a global sea-level rise exceeding 0.2 m by 2040 (Jevrejeva et al., 2016). The most common approach to these threats is shoreline hardening (Gittman et al., 2015) which involves the construction of rigid, engineered structures to stabilize shorelines and protect infrastructures. The main types of shorelines hardening structure include bulkheads, seawalls and breakwaters. Bulkheads are vertical walls typically deployed in environments with low to moderate hydrodynamic energy, built directly on the shoreline to retain soil and prevent land loss due to wave action. Seawalls are robust, vertical or sloped structures designed to protect the shoreline in high-energy environment by reflecting and/or dissipating wave energy and preventing coastal flooding and erosion. Breakwaters are offshore, barrier-like structures built nearly parallel to the coast to dissipate wave energy and reduce current velocity before the waves reach the shore. While effective at protecting human settlements, these structures artificialize coastal ecosystems by fundamentally modifying natural habitats, leading to significant socio-ecological changes such as loss of biodiversity and reduction of the intertidal habitat (Bernatchez and Fraser, 2012; Gittman et al., 2016).

Salt marshes are important coastal ecosystems developing in low-energy shorelines, in temperate and cooler regions. They provide some of the most valuable ecosystem services (Barbier et al., 2011). Indeed, these intertidal habitats play a crucial role in coastal protection by stabilizing sediments, attenuating wave energy, and enhancing sediment deposition, thereby helping to mitigate coastal erosion (Duarte et al., 2013; Gedan et al., 2009). Due to their high productivity and the capacity to enhance sediment deposition, salt marshes contribute to long-term carbon sequestration (Duarte et al., 2013), provide habitat for numerous organisms (Barbier et al., 2011), and play a significant role in nutrient cycling, particularly by acting as nitrogen sinks (Sousa et al., 2008; Valiela and Cole, 2002). This function has gained increasing importance in recent decades with the rise of coastal eutrophication (Pascal et al., 2022; Rabalais et al., 2009). Eutrophication, the excessive supply of water bodies with organic matter (Nixon, 1998), largely driven by nitrogen enrichment, can lead to significant ecosystem degradation (Rabalais et al., 2009).

In this context, understanding nitrogen removal processes in salt marshes is crucial when assessing the impacts of shoreline hardening on coastal ecosystem functioning (Sousa et al., 2008; Valiela and Cole, 2002). Shoreline hardening can alter sediment characteristics, hydrodynamic conditions, and benthic communities (Gittman et al., 2016), all of which influencing nitrogen cycling processes (Albert et al., 2021; Stief, 2013). Two key microbial processes drive nitrogen removal from soft-sediment systems. (i) Denitrification which involves the oxidation of organic carbon (OC) using nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$) as an electron acceptor under anaerobic conditions, resulting mainly in the production of dinitrogen gas (N₂), but can also release as nitrous oxide (N₂O), a potent greenhouse gas, together with carbon dioxide (CO₂). (ii) Anaerobic ammonium oxidation (anammox), produces N₂ through the oxidation of ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) coupled to the reduction of nitrite (${\mathrm{NO}}_{\mathrm{2}}^{-}$) although in many organic-rich coastal sediments its contribution can be smaller than that of heterotrophic denitrification (Dalsgaard et al., 2005; Devol, 2015). In addition to denitrification and anammox, dissimilatory nitrate reduction to ammonium (DNRA) can also occur, competing for nitrate and retaining nitrogen within the system.These processes are further influenced by benthic macrofauna, which play key roles in organic matter (OM) mineralization and nitrogen cycling through bioturbation (sensu Kristensen et al., 2012), enhancing solute exchanges and nutrient regeneration at the sediment-water interface (Pascal et al., 2016; Snelgrove et al., 2018). Therefore, alterations to coastal environments through shoreline hardening could significantly impact these nitrogen removal pathways, with potential consequences for coastal water quality and eutrophication mitigation. Given the vulnerability of salt marshes to such modification, understanding how nitrogen cycling is impacted is critical. Salt marshes face increasing threats worldwide (Mcowen et al., 2017) including shoreline hardening (Gedan et al., 2009) that can exceed 50 % of the total shoreline length in some regions (Bulleri and Chapman, 2010; Gittman et al., 2015). This extensive shoreline hardening contributes to a phenomenon known as coastal squeezing (Doody, 2013) where natural habitats are trapped between rising seas and hardened shoreline, resulting in habitat loss and degradation (Bernatchez and Quintin, 2016; Bozek and Burdick, 2005). Previous research on bulkhead implementation found no direct effects on nitrogen removal rates per unit area, but showed that these structures can reduce the overall nitrogen removal capacity of the ecosystem by decreasing the total area of salt marshes through coastal squeezing, rather than impacting nitrogen removal process rates (O'Meara et al., 2015). However, breakwaters differ from bulkheads as they are built offshore. When emergent and built with sufficient vertical distance between the crest of the breakwater and the high-tide level, breakwater create low-energy environments on their landward side (Wang et al., 2021). These conditions promote the sedimentation of fine, organic-rich particles (Gittman et al., 2016), which can substantially alter macrobenthic communities, and biogeochemical processes. Such changes may enhance favourable conditions for heterotrophic denitrification by increasing OC availability and the extent of anoxic sediment, potentially increasing nitrogen removal rates (Albert et al., 2021; Stief, 2013). Yet, the extent and underlying mechanisms of these impacts remain poorly understood, limiting our ability to predict how different shoreline hardening strategies will affect nitrogen removal processes and associated ecosystem services.

Figure 1 Study site. (a) Location of the study area (open square) in the St. Lawrence system (©Stamen Design (CC BY 4.0)/©Stadia Maps; OpenStreetMap), (b) location of the study site (black rectangle) in the bay of Rivière-du-Loup (©Gouvernement du Québec (CC BY 4.0; MNRF, 2025)) and (c) aerial view of the sampling location around the breakwater (from Zouaghi et al., 2024) with saltmarsh areas highlighted by yellow contours. Reference sites: Vegetated Sediment (VS), and Unvegetated Sediment (US); Impacted sites: Southern area landward of the breakwater (LBWS), Northern area landward of the breakwater (LBWN); and Intermediate site: Seaward of the breakwater (SBW).

This study aims to compare sediment biogeochemistry and nitrogen cycling across contrasting areas within an intertidal salt marsh seascape, including areas dominated by salt marsh vegetation, unvegetated sediments, and areas located directly seaward and landward of a breakwater. We hypothesize that the presence of the breakwater causes organic-rich sediment accumulation on its landward side, which alters macrobenthic community structure, increases the sediment oxygen demand and the rates of nitrogen removal processes compared to areas farther from the breakwater.

2 Methods

2.1 Study site

The St. Lawrence Estuary area (Québec, Canada) has already seen an accelerated sea-level rise since the early 1990s (Gibeault et al., 2016), with anticipated relative sea-level rise to 0.67 m by 2100 (Savoie-Ferron et al., 2021), threatening infrastructure and human populations with erosion and flooding. As in many regions worldwide, coastal hardening has been implemented in response to these threats; recent estimates indicate that 10 % of the Estuary and Gulf of St. Lawrence coastline is now artificialized by rigid structure (Sauvé et al., 2020). Among these, the Rivière-du-Loup embayment (47°50^′29.30^′′ N, 69°33^′11.84^′′ W; WSG 4), located on the south shore of the Upper St. Lawrence Estuary (Fig. 1), provides a representative example of coastal hardening.

The Rivière-du-Loup embayment receives freshwater inputs from a river located in its northern part, with a mean annual discharge of 14.7 m³ s⁻¹ (from 3.4 to 141.1 m³ s⁻¹) and mean annual concentrations of 2.7, 13.7 and 0.5 µmol L⁻¹ for ${\mathrm{NH}}_{\mathrm{4}}^{+}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$, respectively (MELCCFP, 2026). This bay is covered by sea ice during winter and is characterized by a semi-diurnal tide with an amplitude between 0.5 and 3.2 m. A salt marsh partially covers its intertidal area with cordgrass (Sporobolus alterniflorus, formerly known as Spartina alterniflora) in the high foreshore. The shoreline has experienced erosion rates ranging from 2.3 to 6.4 m yr⁻¹ since 1984 (Dionne, 2007), endangering a major highway in addition to reducing the salt marsh area by approximately 62 750 m² (Zouaghi et al., 2024). In 2009, the Quebec Government implemented a 300 m-wide breakwater made of boulders to mitigate erosion. It was estimated that the construction of the breakwater per se induced a loss of 6318 m² of salt marsh cover between 2009 and 2011 (Zouaghi et al., 2024). A distinctive feature of this breakwater is that the access path for its construction (perpendicular and towards the center of the breakwater; Fig. 1c) was never removed and now acts as a jetty that is only submerged during high tide events (Zouaghi et al., 2024). In 2010, following the implementation of the breakwater, a storm event coupled with high tides transported and deposited a substantial amount of sediment, which became trapped landward of the breakwater (Zouaghi et al., 2024).

2.2 Experimental design

To assess the impacts of the breakwater structure on ecosystem functioning and nitrogen cycling, we adopted a space-for-time approach comparing sites with and without breakwater influence within the same coastal system. This approach assumes that spatial differences in biogeochemical processes and ecosystem functioning reflect the long-term impacts of the breakwater. Five sites were selected to assess how the breakwater affects sediment characteristics and biogeochemical processes within the intertidal zone (Fig. 1). Two unvegetated sites located landward of the breakwater, the southern (LBWS) and northern (LBWN) areas, were classified as “impacted”. Two other sites outside the breakwater, a vegetated (VS) and an unvegetated (US) sediment areas served as “reference” conditions. A site located seaward of the breakwater (SBW) was designated as “intermediate” to help to distinguish the breakwater's effects from natural environmental gradients between reference and impacted sites. We acknowledge however that the SBW site could also capture potential offshore effects of the breakwater.

2.3 Seawater, sediment and macrobenthic community sampling procedure

The main sampling campaign was conducted in June 2021, at the beginning of the salt-marsh growing season following winter ice-cover breakup. An additional campaign was conducted in June 2022 to address the specific data gaps described below. This period was chosen because salt marsh plants are entering their active growth phase, during which root oxygen release and associated rhizosphere processes are expected to influence sediment nitrogen cycling (Koop-Jakobsen and Giblin, 2010; Zheng et al., 2016). Additionally, the relatively short leaf canopy at this stage allowed intact vegetated sediment cores to be incubated without clipping the aboveground shoots, thereby avoiding disturbance associated with shoot removal. Sea water temperature and salinity during sampling were similar between years: 9 °C and 19.6 in 2021 and 10 °C and 18.7 in 2022. A 20 L volume of seawater was collected at high tide on the sampling day in the bay. This seawater was filtered on-site using a GF/F filter and aerated in an opaque container for later use. At each of the five sites, eight large sediment cores (10 cm diameter, ∼ 15 cm of sediment and ∼ 15 cm of air headspace) were manually collected at low tide using transparent acrylic liners (30 cm long), with compaction < 0.5 cm. The sediment surface was immediately capped with Styrofoam^® to prevent disturbance during transport. Additionally, five small sediment cores (2.5 cm diameter, 1 cm long) were collected at each site for analysis of surficial sediment characteristics (see Sect. 2.4), and stored at − 20 °C until analysis. Six of the eight large sediment cores were kept intact and immediately transported to the laboratory to assess total benthic fluxes. Of the two remaining cores, one was used to obtain O₂ microprofiles at the sediment-water interface and subsequently sliced for particle phase analysis, while the other was dedicated for porewater extraction. To assess macrobenthic community structure, sediment from three box corers (40 cm × 40 cm × 15 cm depth) was collected at each site. This sediment was sieved through a 0.5 mm mesh and preserved in a 4 % buffered formaldehyde solution for further analysis. Due to sampling failure at SBW, total benthic fluxes were assessed on only four sediment cores, and no isotope pairing incubation, porewater extraction, or macrofauna sampling were performed at this site. To address these data gaps, an additional campaign was conducted in June 2022 to obtain O₂ microprofiles at SBW and to provide a second measurement of selected sediment characteristics at all sites. All benthic flux and isotope-pairing incubations, porewater extractions, and macrofauna sampling were conducted in 2021. The sediment characteristics measured in both years showed consistent spatial patterns, supporting the temporal representativeness of the site comparisons (Figs. S1 and S2 in the Supplement).

2.4 Characterization of the sediment

Sediment characteristics were analyzed using samples from the top layer (0–0.5 cm) of the five small sediment cores collected at each location. One core was dedicated to measuring sediment grain size (D₅₀), porosity (φ), and organic carbon (C_org) and total nitrogen (TN) contents, as well as C_org natural isotope abundance (δ¹³C_org). From this core, D₅₀ was assessed by laser diffraction (Malvern Instruments, 2 µm detection limit) after OM digestion with hydrogen peroxide. Porosity was determined by measuring water loss after freeze-drying, corrected for sea salt content, using a dry sediment density of 2.65 g cm⁻³. Prior to the C_org and δ¹³C_org analyses, inorganic carbon was removed by fuming of sediments in a closed container in the presence of an open beaker of concentrated HCl. C_org and TNof the sediment were determined using an elemental analyzer (Carlo Erba NC2500), after which δ¹³C_org was measured using a Micromass Isoprime 100 isotope ratio mass spectrometer coupled to an Elementar Vario MicroCube elemental analyzer. Carbon isotope ratios are reported using the delta notation, relative to the Vienna Pee Dee Belemnite standard. The remaining sediment from this first core, along with the other four cores, was used to determine the OM content via loss on ignition (LOI, n=5), as the weight loss after combustion at 550 °C during 4 h. Due to storage failure, LOI was not assessed at US.

2.5 Oxygen and nutrient porewater concentrations

A 500 mL volume of aerated seawater was gently added at the top of the two large cores dedicated to O₂ microprofiles and porewater extraction and the overlying water was then aerated. The vertical oxygen concentration within the sediment was measured within an hour of sampling using Clark-type microelectrodes (OX100, Unisense with a repeatability better than 2 % as specified by the manufacturer) in one core. The core was maintained at in situ temperature (i.e., 9 °C) using ice packs and the overlying water was stirred and kept at 100 % of air saturation by gently bubbling it during all the following operations. For each location, we conducted five depth profiles of O₂ concentration, randomly distributed within a single large sediment core. Microelectrodes were inserted into the sediment at 100 µm intervals using a motor-driven micromanipulator (MU1, Pyroscience). The sediment surface was identified as the depth corresponding to the abrupt change in O₂ concentration and was assigned a depth of 0. A two-point linear calibration of the electrodes was achieved from measurements of the O₂ concentration of air bubbled overlying water (100 % saturation) and the anoxic zone of the sediment (0 %). Oxygen diffusive fluxes (DOU in $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) were computed from the O₂ gradient between 0 and 500 µm depth below the sediment-water interface using Fick's first law of diffusion, as denoted in the following equation:

$$\begin{array}{}\text{(1)}& \text{DOU}=\mathit{\varphi }\times {\displaystyle \frac{D_{\mathrm{0}}}{(\mathrm{1}-\ln (\mathit{\varphi }){)}^{\mathrm{2}}}}\times {\displaystyle \frac{\mathit{\delta }C}{\mathit{\delta }z}}\end{array}$$

where ϕ is the sediment porosity, D₀ the diffusion coefficient (m² h⁻¹) and $\frac{\mathit{\delta }C}{\mathit{\delta }z}$ is the concentration gradient ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{m}}^{-\mathrm{1}}$).

Oxygen penetration depth (OPD) was determined from the micro-profiles as the depth at which O₂ concentration was < 1 % of the concentration at the sediment-water interface. Porewater extraction was performed on another large core (< 2 h after the core recovery) using CSS rhizons (Rhizosphere Research Products). Rhizons allow the extraction of porewater under a closed vacuum, minimizing the contact with atmospheric oxygen. After removing the overlying water to minimize potential contamination of surface porewater samples, Rhizons were horizontally inserted in the sediment core at 0.5, 1, 1.5, 2, 3, 4, 5, 7, and 10 cm depth. Porewater extraction was limited to 8 mL per depth to reduce overlap between adjacent sampling volumes. Nevertheless, some overlap between the volumes sampled by neighbouring Rhizons cannot be excluded and may have smoothed concentration gradients in the uppermost centimetre. After extraction, porewater samples were stored at −20 °C until nutrient analyses.

2.6 Reactive Fe oxyhydroxide concentrations

Following O₂ microprofiling, the same large core was sectioned, and homogenized freeze-dried sediment from each depth interval was used for reactive Fe oxyhydroxide extractionusing a buffered ascorbate solution. The term “reactive” refers to the easily reducible fraction of Fe oxides, typically amorphous Fe(III) oxides (Hyacinthe et al., 2001; Kostka and Luther, 1994). The ascorbate reagent consisted of 10 g of sodium citrate and 10 g of sodium bicarbonate dissolved in 200 mL of Nanopure^® water, with 4 g of ascorbic acid added to reach pH 8 (Kostka and Luther, 1994). Approximately 250 mg of ground, freeze-dried sediment was digested in 15 mL of the ascorbate reagent. After 24 h of shaking at room temperature, the sediment slurries were centrifuged, and a fraction of the supernatant was diluted with a 1 M HNO₃ solution before Fe was analyzed using Microwave Plasma Atomic Emission Spectroscopy (4210 MP-AES, Agilent). Fe_asc concentrations were corrected for sediment porosity and adjusted for salinity, with a measurement reproducibility better than 2 %.

2.7 Intact sediment core incubations and isotope pairing experiments

Six sediment cores were fitted with PVC lids equipped with a magnetic stirrer. Prior to the initial incubation, the overlying water in the sediment cores was carefully filled with previously collected and bubbled seawater, and the cores were left open and bubbled overnight for acclimation. The following day, the overlying water was completely replaced with fresh, previously collected and bubbled seawater, and left open while being gently bubbled with air for two hours before starting the experiment to restore initial oxygen and nutrient conditions. The dissolved O₂ concentration was measured in real time using an optode (O₂ sensor SP; PreSens with a repeatability better than 0.5 % as specified by the manufacturer) connected to a multichannel O₂ meter (OXY-10 SMA; PreSens) through optical fibre cables. Incubations were performed in the dark at the in situ water temperature (i.e., 9 °C). The incubation duration was adjusted so that O₂ saturation in the overlying water never dropped below 80 % of its initial value (typically between 1.5 and 4 h). During the incubations, discrete water samples were collected each hour, filtered (0.2 µm) and stored at −20 °C for subsequent nutrient analyses. Water removed by sampling was automatically replaced with previously collected and bubbled seawater. Overall, cumulative volume of replacement waters corresponded to < 5 % of the overlying water and its dilution effect was thus considered as negligible. After the initial incubation period, two thirds of the overlying water volume in each sediment core was replaced with freshly aerated site water to restore initial oxygen and nutrient conditions and minimize carry-over from the preceding flux incubation. The sediment cores were then left open, and the overlying water was gently bubbled with air for 2 h. Then, an initial water sample was extracted from the overlying water of all six sediment cores to determine the ambient ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration before the addition of ${}^{\mathrm{15}}\mathrm{N}\text{-}{\mathrm{NO}}_{\mathrm{3}}^{-}$ stable isotope. A calculated volume of a 20 mM ${}^{\mathrm{15}}\mathrm{N}\text{-}{\mathrm{NO}}_{\mathrm{3}}^{-}$ solution (99 % atom ¹⁵N, Sigma-Aldrich) was then introduced into the overlying water of each core. For this tracer-addition experiment, the added volume was designed to achieve a nominal concentration equivalent to either 50 % or 100 % of the ambient concentration, with three cores per level. A 3 h pre-incubation period was implemented to ensure sufficient diffusion of ${}^{\mathrm{15}}\mathrm{N}\text{-}{\mathrm{NO}}_{\mathrm{3}}^{-}$ to the sediment nitrate reduction zone (Nielsen, 1992), including shallow rhizosphere and bioturbated layers. During this period, an air bubbling system was used to keep the overlying bottom-water homogeneously saturated with air, maintaining the in situ conditions. Following the pre-incubation, water samples were collected to determine the exact concentration of ${}^{\mathrm{15}}\mathrm{N}\text{-}{\mathrm{NO}}_{\mathrm{3}}^{-}$ and the initial ²⁹N₂ and ³⁰N₂ concentrations. Then, the incubations were performed as described above. Discrete water samples were collected using 12 mL Exetainers, poisoned using ZnCl₂ and analysed within a week for ²⁹N₂ and ³⁰N₂ measurements. At the end of the incubations, the remaining water was mixed with the top 5 cm of the sediment column. From the slurry, 12 mL of sample was transferred to an Exetainer to which 200 µL of 7 M ZnCl₂ were added for further ${}^{\mathrm{15}}\mathrm{N}\text{-}{\mathrm{NH}}_{\mathrm{4}}^{+}$ analysis to assess DNRA rates.

For porewater sample, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (quantification limit 0.6 µM), ${\mathrm{NO}}_{\mathrm{2}}^{-}$ (quantification limit 0.05 µM), ${\mathrm{NH}}_{\mathrm{4}}^{+}$ (quantification limit 0.1 µM) and ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ (quantification limit 0.1 µM) concentrations were analyzed using an autoanalyzer AA500 (Seal analytical) following Aminot et al. (2009). For water incubation samples, where sample volume was limited to avoid dilution due to replacement water, ${\mathrm{NH}}_{\mathrm{4}}^{+}$ was determined manually according to Hansen and Koroleff (1999), while NO_x (${\mathrm{NO}}_{\mathrm{2}}^{-}+{\mathrm{NO}}_{\mathrm{3}}^{-}$) were analyzed by continuous flow analysis following Anderson (1979) both with measurements reproducibility better than 5 %. Analytical accuracy was assessed using certified reference materials (MOOS-4, NRC Canada).

The isotope ²⁸N₂, ²⁹N₂ and ³⁰N₂ measurements were completed within one week of collection using a membrane inlet mass spectrometer (MIMS, Bay Instruments) (Kana et al., 1994)following Kana et al. (2006). Calibration factors for standard sample were computed based on solubility equations (Hamme and Emerson, 2004). Analytical precision for ${\mathrm{N}}_{\mathrm{2}}/\mathrm{Ar}$ ratios (coefficient of variation of five replicates standards) was consistently < 0.5 %. To quantify the production of ${}^{\mathrm{15}}\mathrm{N}\text{-}{\mathrm{NH}}_{\mathrm{4}}^{+}$, the water samples were first gently bubbled with He for 5 min to remove all ¹⁵N-N₂ produced during the incubation. The water samples were then treated with 200 µL of 2 M KCl to desorb any particle-bound ${}^{\mathrm{15}}\mathrm{N}\text{-}{\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ followed by the addition of alkaline hypobromite iodine to oxidize ${}^{\mathrm{15}}\mathrm{N}\text{-}{\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ to N₂ prior to analysis using MIMS as described in Yin et al. (2014).

The slope (in $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{h}}^{-\mathrm{1}}$) from the linear regression of the changes in O₂, nutrients and ^29,30N₂ concentrations over incubation duration was verified and used to compute benthic fluxes as follows:

$$\begin{array}{}\text{(2)}& J_{\text{benthic}}=\text{ slope }\times {\displaystyle \frac{V}{A}}\end{array}$$

where J_benthic is the flux across the sediment-water interface ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$), V is the volume of overlying water (m³), and A is the surface of the sediment-water interface (m²).

The validity of the original isotope pairing technique (Nielsen, 1992) was tested by checking the absence of increase of D₁₄ with the amount of added ${}^{\mathrm{15}}\mathrm{N}\text{-}{\mathrm{NO}}_{\mathrm{3}}^{-}$ as well as the linear increase of the D₁₅ with the amount of added ${}^{\mathrm{15}}\mathrm{N}\text{-}{\mathrm{NO}}_{\mathrm{3}}^{-}$ (Risgaard-Petersen et al., 2003). The tracer addition experiment provides no evidence of a substantial anammox contribution. DNRA was assessed independently from ${}^{\mathrm{15}}\mathrm{N}\text{-}{\mathrm{NH}}_{\mathrm{4}}^{+}$ production in the terminal sediment slurry. Since no ${}^{\mathrm{15}}\mathrm{N}\text{-}{\mathrm{NH}}_{\mathrm{4}}^{+}$ production was detected DNRA was considered negligible. Both processes are therefore not considered further in the manuscript and denitrification rates, including the contributions supported by ${\mathrm{NO}}_{\mathrm{3}}^{-}$ from the overlying water (D_w) and by ${\mathrm{NO}}_{\mathrm{3}}^{-}$ produced through sedimentary nitrification (D_n), were computed on the six sediment cores following Nielsen (1992).

2.8 Macrofauna identification

Macrofauna (epi- and infauna) was sorted and identified to the lowest possible taxonomic level (typically species level). Biomass (blotted wet mass; ± 0.0001 g) was determined as well. Using this macrofauna dataset, total density, total biomass and species richness were calculated. In this study, we identified six biological traits characterized by multiple categories to describe the organisms' functional behaviors (Table 1). Taxa were assigned trait modalities using a standard fuzzy coding approach which allowed each taxon to represent multiple trait modalities and thus account for intraspecific variation. Scoring ranges from 0 to 3 were adopted with 0 corresponding to no affinity for a given trait modality, 1 and 2 reflecting a low and high affinity, respectively, and 3 an exclusive affinity. Information on trait expression for all taxa were extracted from Queiros et al. (2013) and Degen and Faulwetter (2019). When information was not available for a given taxon, information was obtained from that of closely related species at the nearest taxonomic level. The community trait expression in each sample was computed from the multiplication of the species density matrix by trait categories matrix.

Table 1 Biological traits. List of traits, trait categories and associated abbreviations.

2.9 Data analyses

Data analyses were performed using R software (R Core Team, 2024). Data values are given as means ± SD. First, linear models (lm function from stats package; R Core Team 2024) were fitted to test the effects of Site (fixed factor; five levels) on the total macrofauna density, biomass and species number, DOU, OPD, as well as O₂, nutrients and N₂ total benthic fluxes. Test assumptions were assessed for each model (plotQQunif function from DHARMa package; Hartig, 2026). Effect sizes were calculated as generalized eta squared (ges) to quantify the magnitude of Site effects, with ges values of 0.01, 0.06, and 0.14 representing small, medium, and large effects, respectively. Both p-values and effect sizes are reported (Table S2 in the Supplement) given the limited replication (n=3 per site), which may affect statistical power. If a significant effect of the Site factor was observed (p < 0.05), pairwise tests (Tukey contrast) were performed to characterize the differences between its levels. Second, macrofauna density data were transformed using Hellinger transformation. Sites were grouped using a hierarchical clustering analysis with a ward agglomeration method (hclust function from stat package; the number of clusters was determined using the HCPC function from FactoMineR package; Lê et al. 2008). Patterns of macrobenthic community assemblage were visualized using principal component analysis (PCA, RDA function from vegan package; Oksanen et al. 2022). In addition, the difference in macrobenthic community between clusters was tested using a multivariate analysis of variance (anova.cca function from vegan package). Finally, functional diversity was assessed using the dbFD function from the FD package (Laliberté and Legendre, 2010), computing functional richness (FRic) and functional evenness (FEve). Community-weighted mean trait values (CWM) were calculated and visualized using PCA as described above. Site effects on functional indices were tested using linear models as described above.

2.10 Methodological considerations

Sediment characteristics were measured with variable replication levels (n=1 to 5) due to analytical constraints. LOI was replicated five times, some parameters were duplicated, while others were measured once per station. The representativeness of these measurements is supported by two observations: (1) low variability among replicated measurements where available, and (2) strong consistency between highly replicated LOI values and corresponding Corg and TN measurements (Figs. S1 and S2).

Oxygen microprofiles and porewater nutrients were each measured on a separate single sediment core per location. This approach was chosen because biogeochemical heterogeneity of the sediment is often associated with microenvironments (Stockdale et al., 2009). Such approach is considered qualitatively robust (Bull and Taillefert, 2001). Nevertheless, since this sampling design may involve pseudo-replication, statistical results on these measurements should be interpreted cautiously. Measurement representativeness is supported by: (1) consistency among sediment characteristics, oxygen microprofiling, and porewater nutrient profiles measured in three independent cores per location, and (2) absence of large infauna in cores, limiting potential bias from bioturbation.

Our benthic flux incubations were conducted in the dark to standardize experimental conditions and ensure robust comparisons among sampling sites. This approach inherently excludes potential in situ effects of microphytobenthos or S. alterniflorus photosynthesis, both of which can significantly influence sediment oxygenation and nutrient uptake during daylight. However, in intertidal zones near breakwaters, high turbidity and suspended sediment loads, often limit light penetration and thereby reduce photosynthetic activity (Pratt et al., 2014). Moreover, belowground oxygen release by marsh roots correlates with photosynthesis but can persist in darkness (Koop-Jakobsen et al., 2018). As such, our dark incubations represent standardized dark benthic fluxes, under low-light conditions typical of these habitatsbut do not capture the full spectrum of benthic fluxes. Future studies combining dark and natural-light incubations would better characterize daytime flux dynamics, particularly where benthic primary producers receive substantial light during immersion. Denitrification rates were derived from the linear accumulation of ²⁹N₂ and ³⁰N₂ in the overlying water, indicating consistent transfer of labelled N₂ to the overlying water during the incubation. Because the terminal inventory of labelled N₂ dissolved in sediment porewaters was not quantified, we cannot exclude that some labelled N₂ produced during the incubations had not reached the overlying water by the end of the experiment. Nevertheless, all incubations were conducted and analysed using the same standardized protocol, providing consistent measurements of labelled N₂ accumulation across sites. Future studies could quantify the end point porewater inventory of labelled N₂ to assess its contribution to total production.

Table 2 Surficial Sediment characteristics. Median sediment grain size (D₅₀, in µm), porosity (ϕ, volume ratio), loss on ignition (LOI in %), organic carbon content (C_org, in %), total nitrogen content (TN, in %), δ¹³C_org signature (in ‰ vs. VPDB). Values are means ± SD with replication and sampling year (SY) indicated by superscripts (^a n=2, SY: 2021–2022; ^b n=5, SY: 2021; none: n=1, SY: 2021).

3 Results

3.1 Sediment characteristics and biogeochemical properties

Among all sampling locations, surface sediments consisted of silty clay (Table 2). Overall, impacted sites (LBWN and LBWS) showed finer median grain size compared to reference sites (VS and UV). Correspondingly, sediment OC and TN contents tended to be higher in impacted sites, while VS exhibited the lowest OC and TN contents among sites. The δ¹³C_org signatures remained relatively uniform among all sites, suggesting a common mixture of terrestrial and marine carbon sources through the system. The intermediate sites (SBW) showed similar sediment characteristics to the two reference sites.

Figure 2 Sediment biogeochemical propoerties. (a) Diffusive O₂ uptake (DOU) by the sediment and (b) O₂ penetration depth (OPD). (c) O₂ concentration microprofile. (d) porewater ${\mathrm{NO}}_{\mathrm{2}}^{-}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ concentrations at each site (n=1 with 5 pseudo replicates for DOU and OPD). (e) Solid-phase iron oxyhydroxide concentration at each site. Different letters indicate significant (p < 0.05, Tuckey contrast) differences among the different levels of Site (a, b).

Vertical porewater O₂ concentrations showed contrasting dynamics among sites (Table S1 in the Supplement). DOU rates tended to be higher at impacted sites compared to intermediate and reference sites (Fig. 2a). Correspondingly, OPD were systematically shallower (Fig. 2b). Similarly, porewater nutrient distributions varied among sites (Fig. 2d). While ${\mathrm{NO}}_{\mathrm{2}}^{-}$ concentration was near zero throughout the sediment column at all sites, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration consistently reached its peak at 0.5 cm depth, and then sharply decreased below the detection limit around 4 cm deep in impacted sites and around 6 cm deep in reference sites (Fig. 2d). The peak concentration was particularly high at VS. ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ concentration increased with depth with a larger increase at impacted sites compared to reference ones (Fig. 2d). The increase in ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ concentration started from the sediment surface at impacted sites. In comparison, it started deeper (5 cm deep) at US and ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ concentration was near zero along the sampled sediment column at VS (Fig. 2d). Fe_asc showed a decreasing trend from the sediment surface to the bottom for LBWN and LBWS, no trend for US and an increase in the first 3.5 cm then a decrease for VS (Fig. 2e).

Figure 3 Univariate descriptors of macrofauna communities. Total biomass, density and species richness at the sampled sites (n=3 per site). Different letters indicate significant (p < 0.05, Tuckey contrast) differences among the different levels of Site.

3.2 Macrobenthic community structure

A total of 2657 individuals from 12 samples (4 sites × 3 replicates) were sorted, and 14 taxa were identified (Table S2). Total density, biomass and species richness significantly differed among sites (Table S1). However, these univariate metrics do not allow a clear separation between impacted and reference sites; comparisons must instead be made at the site level. Total densities tended to be higher at the sites LBWS than at LBWN, VS and US (Fig. 3). Total biomass was significantly lower at LBWN, LBWS and VS than at US. Finally, species richness was lower at LBWN than at LBWS, VS and US.

Figure 4 Macrofauna community assemblages. (a) Hierarchical clustering and (b, c) principal component analyses (PCA) of Hellinger-transformed species densities. Panels show (b) scaling 1, and (c) scaling 2. Only variables representing at least 20 % of the total variation are shown. Abbreviations: MYA = Mya sp.; MAC = Macoma balthica; GAM = Gammarus spp.; OLI = Oligochaeta; GAS = Gastropoda.

Based on hierarchical clustering on species density, the sites could be separated into three groups (Fig. 4a). Cluster 1 included the impacted sites, cluster 2 corresponded to VS, and cluster 3 corresponded to the US. The statistical significance of these clusters was confirmed by MANOVA (Table S3 in the Supplement). The three clusters exhibited distinct macrofaunal assemblages (Fig. 4b and c). The impacted sites (cluster 1) were characterized by a high density of Oligochaetes and Polychaete Hediste diversicolor, while VS (cluster 2) were characterized by relatively high density of small body grazers (Gammarus spp. and Gastropods). US (cluster 3) were distinguished by a relatively high density of the bivalve Mya sp. Both clusters 2 and 3 showed notable densities of the bivalve Macoma balthica as compared to cluster 1. It is worth noting that while gathered in the same cluster, LBWS showed greater densities of H. diversicolor and Oligochaete than LBWN.

Figure 5 Macrofauna community trait expression. (a) Functional richness (FRic) and (b) functional evenness (FEve); (c, d) PCA of community-weighted mean (CWM). Panels show (b) scaling 1, and (c) scaling 2. Only variables representing at least 20 % of the total variation are shown. See Table 1 for trait abbreviations.

Functional richness (FRic) differed significantly among sites (Table S1) and was higher at VS and US, intermediate at LBWS and lower at LBWN (Fig. 5a). Functional evenness (FEve) did not differ significantly among sites, although a large effect size (Table S1) suggests the functional evenness differences may exist but require greater sampling effort to confirm. Overall, FEve tended to be higher at LBWN, VS and US than at LBWS (Fig. 5b).

Beyond diversity metrics, the community-weighted mean trait composition revealed contrasting spatial patterns (Fig. 5c). Although LBWN and US differed markedly in species composition, they exhibited partial convergence in trait expression (Fig. 5d). US was characterized by larger, more mobile organisms dominated by filter-feeding, predatory and scavenger trait expressions, while LBWN showed a tendency toward this trait profile. In contrast, LBWS and VS maintained clear segregation in trait expression. VS was characterized by epifauna with minimal or surface-restricted bioturbation, while LBWS trait expression was dominated by small-bodied infaunal organisms characterized by subsurface deposit feeding and biodiffusive bioturbation.

Figure 6 Benthic fluxes. (a) O₂, ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ fluxes (in $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) and (b) N₂ flux (D_tot in $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) across the sediment-water interface. The dashed bar plots depict the denitrification of nitrate produced within the sediment through nitrification (D_n) and the percentage represents the proportion of D_tot explained by D_n (b). Different letters indicate significant (p < 0.05, Tuckey contrast) differences among the different levels of Site. Note that positive values correspond to an uptake by the sediment, and negative values correspond to a release by the sediment. n=6 per site except at SBW where n=4.

3.3 Benthic fluxes

Total O₂, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ fluxes differed significantly among sites (Table S1). Total O₂ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ fluxes were always negative, indicating an uptake by the sediment while positive total ${\mathrm{NH}}_{\mathrm{4}}^{+}$ flux indicated a release by the sediment (Fig. 6a). Similarly to the macrobenthic community structure results, comparison should be made at the site level. Overall, all benthic fluxes followed the same pattern with higher fluxes at LBWS and VS than at LBWN, SBW and US (Fig. 5a). Total denitrification (D_tot) also differed significantly among sites (Table S1). Similarly to benthic fluxes, D_tot was highest at sites VS and LBWS and lowest at sites LBWN and US (Fig. 6b). The contribution of sedimentary nitrification to total denitrification (coupled nitrification–denitrification; D_n) was higher than 50 % at sites LBWN, LBWS and VS and lower than 10 % at site US.

4 Discussion

4.1 Impact of breakwater construction on benthic habitats

Our study highlights marked alterations in surficial sediment characteristics, with impacted sites (i.e., LBWS and LBWN) accumulating finer and more organic-rich sediment than reference and intermediate sites. These observations are consistent with expectations for low-energy depositional environments landward of high-crested breakwaters (Da Silva et al., 2022; van Rijn, 2011; Wang et al., 2021). In contrast, sediments at the intermediate site (SBW) resembled those at the reference sites (i.e., US and VS), suggesting minimal breakwater influence in this area. The consistent sedimentary δ¹³C_org signature among sites suggested a shared origin of organic carbon, likely reflecting a mixture of terrestrial and marine inputs typical of intertidal zones. This homogeneity suggests that low-energy hydrodynamic conditions landward of the breakwater enhance deposition of fine, organic-rich particles without selective retention of organic matter from specific sources. Therefore, differences among sampled areas are driven by organic carbon quantity, not origin.

Geochemical measurements further support this pattern. Sediment cores from impacted sites showed up to 2.3 times higher diffusive O₂ uptake (DOU) and shallower O₂ penetration depth (OPD) than reference sites (Fig. 2a and b), consistent with enhanced organic carbon mineralization driven by higher organic carbon content (Glud, 2008). These results align well with our porewater nutrient profiles extracted from another sediment core. The decrease in porewater ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration with depth in all sites indicates a consumption by by microbial processes within the sediment (Fig. 2d), consistent with our measured ${\mathrm{NO}}_{\mathrm{3}}^{-}$ benthic fluxes (Fig. 6). This consumption results from the denitrification process, which has been detected at all sites (Fig. 6). VS exhibited a subsurface ${\mathrm{NO}}_{\mathrm{3}}^{-}$ peak, suggesting active nitrification consistent with rhizosphere oxygenation (Koop-Jakobsen et al., 2018), this pattern is consistent with the relatively high proportion of coupled nitrification–denitrification detected in incubation experiment. This interpretation is further corroborated by the solid-phase ascorbate extractable iron (Fe_asc) profiles (Fig. 2e). At VS, Fe_asc increase progressively from the surface to a peak at 3–4 cm depth, well below the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ maximum (Fig. 2d), suggesting subsurface oxidation consistent with rhizosphere oxygenation. Although nitrification likely occurs within this oxygenated rhizosphere zone, the absence of a corresponding accumulation of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ at the depth of maximum Fe_asc suggests that newly produced ${\mathrm{NO}}_{\mathrm{3}}^{-}$ is rapidly consumed by denitrification in adjacent anoxic sediment, preventing its accumulation in porewater. (Koop-Jakobsen et al., 2018; Zheng et al., 2016). The quicker ${\mathrm{NO}}_{\mathrm{3}}^{-}$ consumption in the porewater of impacted sites in the upper four centimeters compared to six centimeters in reference sites suggests a more reduced condition in impacted sites. Similarly, deeper sediment layers in impacted sites exhibit greater increases in ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ compared to reference sites, indicating larger ammonification driven by their higher organic matter content, along with ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ release, typically associated with the reductive dissolution of iron-oxyhydroxides in anoxic sediment layers (Slomp, 2011). Consistently, Fe_ascconcentrations at impacted sites are highest at the sediment surface and decrease rapidly with depth (Fig. 2e). This indicates that aerobic processes, including nitrification, are restricted to the very surface of the sediment at these sites. Overall, these O₂ and nutrient profiles suggest more reduced conditions and greater nutrient accumulations in the sediment column in impacted sites compared to reference sites. Despite these changes in sedimentary conditions, denitrification rates were not substantially altered, suggesting that distinct sediment biogeochemical conditions can sustain similar nitrogen removal rates.

4.2 Impact of breakwater construction on macrobenthic community structure

Consistent with previous studies on the ecological impact of shoreline stabilization structure (Bulleri and Chapman, 2010; Dugan et al., 2018; Gittman et al., 2016), impacted sites exhibited distinct assemblages compared to reference sites. Impacted sites were notably dominated by high densities of small-bodied H. diversicolor and Oligochaetes, the former being tolerant to organic enrichment, and the latter considered opportunistic under such conditions and often associated with fine-grained, organic-rich sediments (Borja et al., 2000). Both H. diversicolor and Oligochaetes are significant bioturbators (Queirós et al., 2013) that influence sediment structure, biogeochemistry and microbial activity at the sediment-water interface (Kristensen et al., 2012). In contrast, reference sites supported less dense but more diverse communities in terms of species and function. Notably, US appeared to be colonized by larger organisms than VS, as suggested by species traits analysis (Fig. 5d) and the relatively low total density and high total biomass in US (Fig. 3). This difference may reflect the dense root and rhizome network of S. alterniflorus in VS, which could hamper the movement and burrowing of larger infauna (Pascal et al., 2020; van Wesenbeeck et al., 2007). The higher density of biodiffusive bioturbation at impacted sites, and particularly LBWS (Fig. 5c), suggests that some ecosystem processes, such as sediment oxygenation and nutrient regeneration, may persist despite habitat alteration. However, the lower biodiversity and functional richness and evenness, at these sites may limit functional redundancy thereby reducing the resilience of ecosystem functions to additional stressors.

4.3 Potential convergence in ecosystem functioning

Despite marked differences in sediment characteristics, macrobenthic community structure and trait expression between impacted and reference sites (i.e., LBWN vs. US and particularly LBWS vs. VS), total benthic fluxes and denitrification rates remained similar. A similar convergence in process rates despite contrasting biological drivers has been reported in other benthic systems (Woodin et al., 2016). Here, however, it results from distinct mechanisms, namely root-driven oxygenation and bioturbation. These findings, derived from dark incubations, suggest that distinct benthic environments can yield similar functional outcomes regarding oxygen uptake, nutrient cycling and notably denitrification rates. Our measured denitrification rates ranging from 6 to 75 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ fall within the range of reported denitrification rates for temperate intertidal mudflats and salt marshes (Hinshaw et al., 2017; Koch et al., 1992; Poulin et al., 2007; Sousa et al., 2012). Coupled nitrification–denitrification (D_n) was the dominant N-removal pathway at impacted sites and VS, indicating that denitrification relies on locally produced ${\mathrm{NO}}_{\mathrm{3}}^{-}$. Conversely, at US, denitrification was primarily driven by water-column nitrate. The high density of subsurface bioturbators at LBWS (Fig. 5c) likely enhances oxygen and nutrient exchange at the sediment-water interface. The sediment irrigation enhances aerobic mineralisation of the OM, stimulates nitrification fuelling coupled nitrification–denitrification, a pattern supported by our isotope pairing results (Fig. 6). In VS, although macrofauna is present, bioturbation is absent or constrained to the sediment surface (Fig. 5c). However, as evidenced by the deep subsurface Fe_asc peak (Fig. 2e; Sect. 4.1), radial O₂ loss from S. alterniflora roots actively creates deep oxidized zones (Koop-Jakobsen et al., 2018). This rhizosphere oxygenation sustains subsurface nitrification, which is tightly coupled to denitrification, thereby achieving high denitrification rates similar to those driven by bioturbation at LBWS. By contrast, LBWN showed a trait expression tendency toward US that could yield similar benthic fluxes. The low denitrification rates at US likely reflect the relatively low organic carbon availability compared to impacted sites, which limits ammonium regeneration through organic matter mineralisation and consequently limits nitrate production through nitrification, thereby reducing coupled nitrification–denitrification. Moreover, macrofauna at US is dominated by a few large burrowers (e.g., Mya sp.) whose local influence may not be captured at the core scale, leading to low benthic fluxes. Interestingly, despite being both impacted sites and having similar sediment characteristics, LBWN and LBWS exhibited different ecosystem functions based on total benthic fluxes and denitrification rates. This most likely results from the higher density of macrobenthic organisms at LBWS than at LBWN, highlighting the critical role of macrobenthic community activity in shaping biogeochemical processes (Snelgrove et al., 2018). In the absence of dense bioturbators, the high organic carbon content at LBWN promotes reducing conditions at LBWN (Fig. 2d and e) that limit oxygen availability and thereby nitrification, ultimately restricting the nitrate supply available for denitrification. However, the factors driving the variation in macrobenthic community structure between these two areas remain elusive. The absence of a detectable breakwater effect on denitrification rates therefore does not reflect an absence of biogeochemical change, but rather a functional compensation mediated by distinct biological drivers at each site.

Placed in a broader context, although the impacts of coastal protection structures on benthic nitrogen cycling remain poorly documented, our results suggest that these effects are structure specific. Indeed, vertical, hard-engineered structures like seawalls induce coastal squeeze and habitat loss. This reduces the nitrogen-removal capacity at the scale of the seascape, even if unit-area denitrification rates remain unaffected (O'Meara et al., 2015). Conversely, in the present study, our offshore breakwater accumulates organic-rich sediment. This accumulation, associated with dense populations of bioturbators (e.g., H. diversicolor), sustains a bioirrigation pathway that functionally substitutes for the lost sediment oxygenation by roots. While this functional equivalence aligns with observations from nature-based protections like living shorelines (Onorevole et al., 2018; Schoell et al., 2023), our study shows that it can also occur in a completely different, fully artificial environment. Taken together, these results show that even in artificialized settings, nitrogen-removal processes can remain functionally equivalent, as long as biogeochemical mechanisms capable of coupling nitrification and denitrification are maintained. Whether such compensatory mechanisms occur in other artificialized systems will depend on local environmental conditions and biological assemblages.

From a coastal-management perspective, one of the main issues is whether areal benthic fluxes and denitrification rates measured in LBWS can offset, at the seascape level, the functional loss caused by salt marsh vegetation removal. LBWS area (∼ 18 800 m²) is about three times larger than the vegetated area directly lost due to the breakwater implementation (∼ 6300 m²). As a first-order estimate, and assuming that the measured fluxes capture the whole area heterogeneity, LBWS could compensate for these specific functional capacity losses of that part of the salt marsh. However, it is worth noting that this area is still more than three times smaller than the total vegetated area lost due to combined erosion and breakwater implementation (∼ 62 700 m²) so any functional compensation remains only partial at the scale of the whole seascape.

5 Conclusion

Our study highlights the complex relationships between sediment characteristics, macrobenthic communities, and ecosystem functioning in coastal environments impacted by shoreline hardening structures. We observed that some ecosystem functions (i.e., benthic respiration, nutrient cycling, and denitrification rates) can be maintained through several combinations of abiotic and biotic factors. This finding reflects a convergence in ecosystem functioning that can arise from distinct underlying mechanisms. However, it is important to interpret these results carefully. While we observed functional similarities in specific ecosystem functions, natural coastal habitats provide an extensive collection of ecosystem functions and services (Friess et al., 2020; Gedan et al., 2009) that are far from being fully captured by our measurements. The observed functional convergence may contribute to the resilience of some ecosystem functions at the seascape level. However, these localized similarities do not mean that these functions are unaltered, as the underlying mechanisms are completely modified and potentially more vulnerable to additional environmental pressures. Furthermore, it does not imply that artificialized environments can sustain all ecosystem functions and services provided by natural habitats. These findings highlight the challenges in predicting and managing ecosystem responses to shoreline hardening, coastal squeeze and probably other coastal disturbances. It also emphasizes that physical, geochemical, and biological factors should be considered when assessing the functional impact of coastal development projects. Given the increasing pressures on coastal ecosystems (Halpern et al., 2019), future research should assess the ability of artificialized ecosystems to maintain their functions under multiple stressors and investigate their long-term stability under various climate change scenarios. Finally, broader ecological implications of apparent functional convergence in coastal ecosystems including potential cascading effects on higher trophic levels should be explored. This understanding will be crucial for developing more accurate predictive models and effective management strategies in the face of ongoing coastal development and climate change.