Structural complexity and benthic metabolism: resolving the links between carbon cycling and biodiversity in restored seagrass meadows

Due to large losses of seagrass meadows worldwide, restoration is proposed as a key strategy for increasing coastal resilience and recovery. The emergence of a seagrass meadow is expected to substantially amplify biodiversity and enhance benthic metabolism by increasing primary productivity and respiration. Nevertheless, open questions remain regarding the metabolic balance of aging seagrass meadows and the roles benthic communities within the seagrass ecosystem play in overall metabolism.

To address these questions, we investigated a chronosequence of bare sediments and adjacent Zostera marina meadows of 3 and 7 years since restoration alongside a natural meadow located within a high-temperate marine embayment in Gåsö, Sweden. We combined continuous measurements of O₂ fluxes using underwater eddy covariance with dissolved inorganic carbon (DIC) and O₂ fluxes from benthic chambers during the productive season (July). Based on the ratio between O₂ and DIC, we derived site-specific photosynthetic and respiratory quotients, enabling the conversion of eddy covariance fluxes to DIC. We assessed benthic diversity parameters as potential drivers of metabolic flux variability.

We observed high rates of gross primary productivity (GPP) spanning −18 to −82 $\mathrm{mmol}\mathrm{DIC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, which increased progressively with meadow age. Community respiration (CR) mirrored the GPP trend, and all meadows were net heterotrophic (GPP < CR), with net community productivity (NCP) ranging from 16 to 28 $\mathrm{mmol}\mathrm{DIC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. While autotrophic biomass did not increase with meadow age, macrophyte diversity did, elucidating potential effects of niche complementarity among macrophytes on community metabolism. These findings provide valuable insights into how community composition and meadow development relate to ecosystem functioning, highlighting potential tradeoffs between carbon uptake and biodiversity.

1 Introduction

Climate change and concurrent biodiversity loss have motivated restoration of natural ecosystems that can contribute to climate change mitigation and adaptation and at the same time strengthen local biodiversity. One such ecosystem is seagrass meadows, which have suffered substantial losses worldwide during the last century (Waycott et al., 2009; Mckenzie et al., 2020). Due to their foundational role in structuring benthic communities, high productivity and their ability to sequester large amounts of carbon, the restoration of previously lost meadows has been proposed as a low-regret option to address both the climate and biodiversity crises concomitantly (Duarte et al., 2013; Gattuso et al., 2018; Orth et al., 2020; Unsworth et al., 2022). Nevertheless, few studies have assessed whether both these goals are mutually attainable within the same restoration projects or if there are tradeoffs between biodiversity conservation and carbon sequestration.

The mechanisms through which a seagrass meadow modifies carbon flows are manifold, influencing import, export and burial of autochthonous (i.e., seagrass biomass) and allochthonous (i.e., organic matter from other sources) carbon (Duarte and Krause-Jensen, 2017). While sedimentation of allochthonous carbon is largely a passive process ultimately governed by local hydrodynamics, autochthonous carbon sequestration is coupled to the productivity of the seagrass meadow and is thus a function of its metabolic fluxes on timescales ranging from minutes to years (Smith and Key, 1975; Smith and Hollibaugh, 1993; Duarte and Cebrian, 1996). Seagrass community metabolism comprises gross primary productivity (GPP) and community respiration (CR) constituted by autotrophic and heterotrophic respiration. The balance between GPP and CR on a daily basis reflects the net metabolism, hereafter termed net community productivity (NCP = GPP − |CR|). The magnitude and direction of GPP, CR and NCP determine all subsequent carbon flows whereby a positive NCP (net autotrophy) equals the net carbon fixed available for remineralization, burial or export (Duarte and Krause-Jensen, 2017). Contrarily, if NCP is negative (net heterotrophy), the meadow is respiring more organic carbon than is fixed and relies on external or historic inputs of organic matter to sustain metabolism. Empirically assessing community metabolism is thus imperative in order to constrain a carbon budget and infer the potential net effect of a seagrass meadow on carbon sequestration.

The vast majority of metabolism studies in seagrass ecosystems to date are based on oxygen fluxes (Ward et al., 2022). Converting these fluxes into carbon currency often relies on assuming a constant stoichiometric 1:1 ratio between oxygen and dissolved inorganic carbon (O₂:DIC) fluxes which may significantly under- or overestimate actual metabolism (e.g., Barron et al., 2006; Duarte et al., 2010; Turk et al., 2015). For marine sediments, this ratio has been estimated to range between 0.8–1.2 on annual timescales (Glud, 2008, and references therein), but the variability is poorly constrained and likely higher in seagrass systems and on shorter timescales (Turk et al., 2015; Trentman et al., 2023). The discrepancy from a 1:1 ratio between benthic oxygen and DIC fluxes can stem from a wide range of processes, including anaerobic sediment processes, nitrate assimilation, photorespiration, and differences in solubility and air–sea gas exchange rates between O₂ and CO₂ (Weiss, 1970; Trentman et al., 2023). In seagrasses, storage in tissues and transport of oxygen to roots and subsequent radial oxygen loss (ROL) can also contribute to deviations from the theoretical 1:1 relation (Borum et al., 2007; Ribaudo et al., 2011; Berg et al., 2019). Assessing carbon cycling in seagrass meadows without characterizing the marine carbonate chemistry system can thus lead to erroneous conclusions regarding their role in carbon cycling and ultimately their climate change mitigation potential.

Despite the growing number of seagrass restoration projects worldwide, assessments of the effect on benthic metabolism are lacking. To our knowledge, the only research effort that has specifically addressed benthic metabolism in restored seagrass was carried out in the Virginia Coast Reserve, USA (Rheuban et al., 2014a), where a large-scale Zostera marina restoration project commenced in 2001 (McGlathery et al., 2012). Rheuban et al. (2014a) employed a chronosequence approach comprising a bare site and two stages of development since restoration (5 and 11 years) and measured benthic metabolism on diel and seasonal timescales. The authors found that GPP and |CR| increased up to 25- and 10-fold, respectively, with meadow age, and this was consistent through seasons. However, NCP was similar, and slightly negative, between the bare site and the oldest restored meadow on an annual basis, despite the vast differences in autotrophic biomass between the two sites (Rheuban et al., 2014a). Notably, summer metabolism revealed a net autotrophic state in the 5-year-old meadow (NCP), whereas the older, mature meadow (11 years) had much higher metabolic fluxes and net heterotrophy on the order of about −50 $\mathrm{mmol}{\mathrm{O}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (Hume et al., 2011; Rheuban et al., 2014a).

Although GPP often substantially increases during summer in temperate seagrass meadows, so does CR to a similar extent (Ward et al., 2022). Consequently, despite large seasonal variability in photosynthesis and respiration, the metabolic state (NCP) is often relatively stable on an annual basis, granted there are no major ecosystem shifts. Interannual variability in NCP has been related to seagrass die-off and recovery episodes (Berger et al., 2020), and seagrass phenology typically dictates fluxes and metabolic state on intra-annual timescales (e.g., Champenois and Borges, 2012; Rheuban et al., 2014a). However, a seagrass meadow comprises several components that contribute to community metabolic fluxes. Aside from the seagrass itself, these include primary producers such as macro- and microalgae and heterotrophic organisms ranging from macrofauna to bacteria. Together, these make up the fluxes of O₂ and DIC measured in the overlying water column by methods such as aquatic eddy covariance, benthic chambers or open-water mass balance. Isolating fluxes deriving from a single meadow component is difficult in situ, although promising efforts have been made to estimate the role of benthic fauna in meadow metabolism (Rodil et al., 2019, 2020, 2021, 2022). When planting seagrass with the stated goals of obtaining both a carbon sink and a biodiversity hotspot, it is essential to understand the relationship between these two and over what timescales it may change as a meadow develops. It is therefore necessary to employ a holistic approach and to assess biogeochemical and biodiversity parameters in tandem across multiple stages of seagrass growth. Importantly, both autotrophic and heterotrophic components of biodiversity are relevant as they are expected to have contrasting effects on metabolism.

Figure 1 Aerial view and seagrass development stages after restoration. (a) Map showing study location (inset) and drone image of Gåsö bay outlining the approximate locations of the sites. (b) A schematic illustration of seagrass meadow development at the four sites, Bare, 3-year, 7-year and Nat, which represent different stages of meadow development as indicated by the arrow and text in italics. Photo: Eduardo Infantes.

The overarching goal of this study was thus to evaluate the alterations in metabolic fluxes and biodiversity across the transition from bare sediment to a mature seagrass meadow subsequent to active seagrass restoration. We hypothesized that, in an early stage, autotrophic biomass is dominating but total biomass is relatively low, resulting in small diel variability in metabolic fluxes and an overall net autotrophic state. As the meadow grows, fauna colonization occurs alongside organic matter accumulation, shifting the system towards a more balanced metabolic state as |CR| increases relative to GPP. Finally, when the meadow has reached maturity, CR and GPP are tightly coupled in a system characterized by high turnover and a balanced NCP.

To test these hypotheses, we utilized a chronosequence of four stages of seagrass development since restoration, all situated within the same sheltered bay. We employed non-invasive, high-resolution aquatic eddy covariance (EC) alongside benthic chamber (BC) incubations. This allowed us to simultaneously monitor fluxes of O₂ and carbonate chemistry parameters from which we could evaluate daily metabolic fluxes of both oxygen and carbon. Additionally, we investigated multiple aspects of taxonomic and functional diversity among both macrophytes and benthic fauna and assessed surficial sediment carbon stocks to infer short-term impacts of seagrass restoration on both carbon cycling and biodiversity.

2 Methods

2.1 Site description

The study took place between 4 and 20 July 2022 on the island of Gåsö (lat 58.2325° N, long 11.3984° E) located at the mouth of the Gullmar fjord on the NW coast of Sweden (Fig. 1). The area has microtidal characteristics, with an amplitude of 20–30 cm, while the Gullmar fjord is stratified, featuring three pycnoclines occurring between 10 and 50 m (Sundbäck et al., 2004). Surface water salinity naturally fluctuates between 15 and 30 in this region due to the alternating currents of brackish Baltic Sea water and saline North Sea water (Lindahl et al., 1998). The bay of Gåsö is a semi-enclosed bay spanning ∼ 0.3 km² with two narrow inlets and outlets and lacks major surface freshwater sources. Its sheltered position results in a small fetch and a mild wave climate (Fig. 1).

The benthos consists of a natural, subtidal continuous eelgrass (Z. marina) meadow and large patches interspersed with bare sediment occurring between 1 and 4 m depth (Fig. 1; Huber et al., 2022). In 2015 and 2019, as part of the seagrass restoration program ZORRO (https://www.gu.se/en/research/zorro, last access: 1 September 2023), two plots of Z. marina were planted at the same depth (∼ 2 m), using the same planting methodology (single-shoot) and shoot density (16 shoots m⁻²) (Gagnon et al., 2023). These plots thus provided a chronosequence of seagrass meadow ages spanning 3 and 7 years since planting, while the bare-sediment area and the natural meadow corresponded to a “before” state and a mature meadow, respectively. The part of the natural meadow we sampled was estimated to have been naturally colonized by meadow expansion 13–15 years ago (Eduardo Infantes, personal observation, 2021). Altogether, this yielded four sites within 100 m distance from each other representing four different stages in the development of a seagrass meadow (Fig. 1; Table S1 in the Supplement). Importantly, the validity of applying a chronosequence methodology to investigate age-related differences in seagrass metabolism relies on assumptions that the sites compared experienced similar abiotic conditions after planting and during sampling (Fig. S2 and Table S2 in the Supplement). Utilizing adjacent sites within a semi-enclosed bay addresses most of those matters, but, to further control assumptions, we monitored in situ flow velocity, photosynthetic active radiation (PAR), temperature, turbidity, salinity and wind conditions during all deployments and assessed the variation in oxygen fluxes explained by each variable using linear mixed-effects models (see below).

2.2 Benthic fluxes

2.2.1 Aquatic eddy covariance (EC)

The EC system (Berg et al., 2003) consisted of a stainless-steel tripod frame with an acoustic Doppler velocimeter (ADV Vector, Nortek) and a fast-responding oxygen microsensor (430-UHS, PyroScience GmbH) programmed to log data continuously at 16 Hz from co-located measurements of velocity and oxygen concentration. In addition, two PAR sensors (LI-192, RBR) were mounted to the frame where one was facing upwards to record incident light and one was directed downwards to record reflected light. This made it possible to calculate the fraction of absorbed light (fAPAR) during deployments. Dissolved oxygen optodes (miniDOT, PME, and HOBO U26, Onset) were mounted on the leg of the frame and recorded ambient dissolved oxygen (DO) concentration within the canopy at 1 min intervals. In addition, a salinity sensor (HOBO U24, Onset), a turbidity sensor (RBRsolo, RBR) and two light intensity loggers (HOBO Pendant MX, Onset) were located on the frame recording at 1 min intervals.

We deployed the EC at the center of the transplanted plots, resulting in a distance of 20 m between the 3-year and the 7-year site. In order to maintain the same water depth, the bare site was located 46 m from the 7-year and 29 m from the 3-year site, with the natural site 57 m away from the 3-year site, 47 m from the 7-year site and 58 m from the bare site (Table S1 in the Supplement). EC deployments lasted between 44 and 49 h (Table S1 in the Supplement). In between each deployment, data were offloaded, the sensors and frame were cleaned, and the batteries were exchanged as necessary.

2.2.2 Benthic chambers (BCs)

Simultaneously with EC, we deployed benthic chambers (n = 6 per site) within 5 m of the EC and randomly emplaced within 1 m of each other (Fig. S1 in the Supplement). We ensured that chambers were positioned downstream of the EC with respect to the dominant current direction so as not to influence the footprint of the EC. Benthic incubation chambers consisted of acrylic cylinders (inner diameter 12.45 cm, length 65 cm) with a custom-made motor running a propeller to mix the water within the chamber and avoid build-up of vertical concentration gradients. We employed pilot tests with dye injection in the laboratory and field to ensure sufficient mixing, and, during incubations, chambers were inserted approximately 20 cm into the sediment. We used transparent (n=3) and opaque (n=3) chambers to simulate day (photosynthesis and respiration) and night (respiration only), respectively. Upon deployment, chambers were left with lids off for about 30 min to allow for suspended sediment to settle.

We drew discrete samples of seawater at the onset and termination of each incubation using two 50 mL syringes attached to 30 cm Tygon^® tubing, inserted through a closable sampling port in the chamber lid. We immediately analyzed seawater in the syringes for pH and dissolved oxygen (DO). pH was measured using an InLab Micro pH electrode with a FiveGo handheld pH meter (Mettler Toledo). The electrode was calibrated using a two-point calibration with standard buffers (pH 7 and 10, Mettler Toledo) at both the onset and termination of the field campaign and was calibrated to a certified Tris buffer in synthetic seawater (Dr. A. Dickson, SIO) at the beginning and end of each sampling day. This was done to account for the effect of salinity and to yield values on the total hydrogen ion scale (pH_T). We measured salinity using a conductivity probe connected to a pH/Cond 340i multimeter (WTW).

We measured DO using a fiber optic oxygen sensor coupled to a FireSting^®-GO2 oxygen meter (PyroScience). A temperature probe was also connected to the FireSting^® to record temperature during each measurement. Seawater from the syringes was then filtered through 0.45 µm Minisart^® sterile syringe filters (Sartorius) and stored in 50 mL Falcon tubes on ice. Upon return to the laboratory, we placed samples for total alkalinity (TA) in a dark container at 4 °C, whereas samples for inorganic nutrients and DOC were frozen (−20 °C) immediately until subsequent laboratory analyses.

We determined TA by open-cell potentiometric titration using an 888 Titrando system with an Ecotrode Plus pH electrode (Metrohm). Samples (40–50 g) were titrated with prepared 0.05 M HCl in ∼ 0.6 mol kg⁻¹ NaCl, corresponding to the ionic strength obtained from the mean salinity of the samples. Accuracy and precision (−1.65 ± 3.76 µmol kg⁻¹) were determined using certified reference material (CRM, batch 200, n=8) provided by Dr. Andrew Dickson at Scripps Institution of Oceanography, San Diego.

We analyzed dissolved inorganic nitrogen (NH₄-N and NO₃-N) using flow injection analysis on a FIAstar 5000 analyzer (FOSS) and phosphate (PO₄–P) using ion chromatography on an 861 Advanced Compact IC (Metrohm). We analyzed dissolved organic carbon (DOC) and total nitrogen (TN) using a TOC-VCPH total organic carbon analyzer (Shimadzu).

We calculated DIC using the package “seacarb” in R (Gattuso et al., 2022) with measured values of pH_T and TA in conjunction with in situ temperature, salinity, pressure and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ as input parameters. We used dissociation constants $K_{\mathrm{1}}^{\ast }$ and $K_{\mathrm{2}}^{\ast }$ from Lueker et al. (2000). We also calculated the saturation state of the CaCO₃ mineral form aragonite (${\mathrm{\Omega }}_{Ar}=\left[{\mathrm{Ca}}^{\mathrm{2}+}\right]\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]/K_{\text{sp}}^{\ast }$) from each sample using seacarb. All solute concentrations were calculated to µmol kg⁻¹ using in situ pressure, temperature and salinity data.

Using incubations with discrete measurements to assess flux rates assumes that concentrations change linearly with time. We verified this assumption ex situ by bringing an intact chamber core from the natural meadow into the laboratory. The chamber was placed in a large water bath with running seawater, and, prior to each incubation, the chamber was saturated with oxygen by bubbling compressed air. We ran multiple dark incubations with continuous logging of dissolved oxygen and temperature (FireSting^®-GO2) combined with multiple discrete measurements of pH (n=4) and TA (n=2).

We used the abovementioned measurements at the onset of incubations to infer mean ambient seawater chemistry at each site (Table S3 in the Supplement).

2.3 Community components

2.3.1 Macrophytes and microphytobenthos

We evaluated seagrass shoot density by placing a 0.25 m × 0.25 m frame randomly in 10 areas of each site and counting seagrass shoots in subareas of 0.016 m² (n=10 per site). In addition, we collected seagrass shoots using a mesh net bag attached to a closable aluminum frame (opening area 0.1156 m², n=3 per site). From these samples, we measured aboveground biomass, shoot density, number of reproductive shoots, leaf length and number of leaves per shoot. We also assessed the taxa and biomass of macrophytes other than seagrass (e.g., red and brown macroalgae). Seagrass belowground biomass, including live and dead roots and rhizomes, was collected using sediment cores (see below). We dried biomass samples at 60 °C for 72 h, and values are reported as dry mass (g m⁻²).

We collected sediment samples to estimate microphytobenthos abundance, using sediment surface chlorophyll a as a proxy. From each sediment core, we used a cutoff 5 mL syringe (Ø = 12 mm) to collect 2 mL sediment from the surface layer. This was repeated 3 times for each core, and we pooled the three samples into one 6 mL sample per core and put them in a 50 mL centrifuge tube covered in aluminum foil to avoid light penetration. The samples were immediately frozen (−20 °C) until subsequent extraction and analysis. After the samples were thawed at 4 °C overnight, we drew a subsample of 2 mL sediment from each sample and weighed and dried it at 60 °C for 72 h to obtain wet weight (g), dry weight (g), dry bulk density (DBD, g cm⁻³) and water content (%). We extracted the chlorophyll using ethanol (99.5 %) and, after diluting and incubating overnight, measured fluorescence using a Turner TD-700 fluorometer (Turner Designs). We calculated chlorophyll a content (g m⁻²) using a modified equation from Hannides et al. (2014).

2.3.2 Benthic fauna

We targeted infauna and epifauna separately, where we collected seagrass epifauna from the mesh net bag samples described above (mesh size ∼ 0.2 mm, n = 3 per site). This approach allowed us to capture the entire community, with cores capturing infauna and slow-moving epifauna and the mesh net approach capturing fast-moving and larger epifauna present in the seagrass canopy. For infauna, we collected sediment cores using polycarbonate cylinders (inner diameter 7.4 cm, length 33 cm, depth 20 cm, n = 6 per site) for determination of infauna and seagrass belowground biomass. Upon return to the laboratory, samples were sieved (0.5 mm) and fixed in 95 % ethanol for subsequent counting and species identification. Fauna was identified to the lowest taxonomic level possible.

2.3.3 Sediment properties

In addition to the sediment cores used for infauna, we collected three additional sediment cores from each site to determine sediment properties. These cores were stored upright and immediately brought back to the lab and sliced into sections at 2, 4, 6, 8, 12 and 16 cm depth. We used the top 0–2 cm section for determination of water content, DBD and porosity (refer to Sect. 2.3.1). After removing all visible rhizomes, roots and shells, we dried all sections at 60 °C for 72 h, homogenized with a pestle and mortar, and analyzed subsamples (5 mL) for organic matter content using loss on ignition (4 h at 520 °C). Subsamples from the top 0–2 cm sediment layer (n = 12) were also analyzed for particulate organic carbon (POC), particulate inorganic carbon (PIC) and total nitrogen (TN) using a vario MAX TN elemental analyzer (Elementar). We pre-treated samples with HCl to remove carbonates, and PIC was obtained by subtracting POC from total carbon (TC). We obtained a linear relationship between OM and POC (POC = 0.47 × OM − 0.88; R² = 0.84, p < 0.001) which we used as a conversion factor to convert remaining OM values to POC and thereby obtain POC values for all core slices. This conversion is based on the assumption that the relationship persists with sediment depth, and this introduces uncertainty in the POC values at depth. We calculated carbon density for each slice between 0 and 12 cm by multiplying POC by surface DBD and integrated across 0–12 cm to obtain the organic carbon stock (POC_stock, g m⁻²) in the upper 12 cm of sediment. Using only DBD values for the top 0–2 cm introduces uncertainty in our depth-integrated POC_stock estimates, but a previous study by Dahl et al. (2023) from the same area showed similar DBD values from 0–11 cm (mean 0.43 ± 0.15 g cm⁻³) that were consistent with sediment depth.

2.4 Data analyses

2.4.1 Flux calculations

We calculated oxygen fluxes in the benthic chambers (BC) as the difference in solute concentration between the onset and termination of each incubation as

$$\begin{array}{}\text{(1)}& F_{{\mathrm{O}}_{\mathrm{2}}}\left(\mathrm{mmol}{\mathrm{O}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}\right)={\displaystyle \frac{\mathrm{\Delta }{\mathrm{O}}_{\mathrm{2}}}{\mathrm{\Delta }t}}\mathit{\rho }h,\end{array}$$

where ΔO₂ is the change in O₂ concentration in mmol kg⁻¹ between start and end of incubation, dt is the duration of the incubation in hours, ρ is the density of the seawater in kg m⁻³, and h is the height of the chambers from the top to the sediment surface in meters. We calculated the flux of salinity-normalized TA (nTA = $\mathrm{TA}/S_{\text{in\hspace{0.17em}\hspace{0.17em}situ}}\times S_{\text{mean}}$; Table S2 in the Supplement) in the same way:

$$\begin{array}{}\text{(2)}& F_{\mathrm{TA}}\left(\mathrm{mmol}\mathrm{TA}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}\right)={\displaystyle \frac{\mathrm{\Delta }n\mathrm{TA}}{\mathrm{\Delta }t}}\mathit{\rho }h.\end{array}$$

Similarly, we used DIC measurements to obtain fluxes as

$$\begin{array}{}\text{(3)}& F_{\mathrm{DIC}}\left(\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}\right)={\displaystyle \frac{\mathrm{\Delta }\mathrm{nDIC}}{\mathrm{\Delta }t}}\mathit{\rho }h-\mathrm{0.5}{F}_{\mathrm{TA}},\end{array}$$

where ΔnDIC is the change in salinity-normalized DIC in mmol kg⁻¹. The subtraction of 0.5 F_TA is to account for the effect of inorganic processes (i.e., calcification/CaCO₃ dissolution) on DIC according to the assumptions that net community calcification affects TA and DIC in a ratio of 2:1 and NCP only modifies DIC (Smith and Key, 1975). F_DIC thus represents the DIC flux stemming from primary production and respiration only.

We calculated the photosynthetic quotient (PQ) and respiratory quotient (RQ) from the average absolute fluxes (i.e., the magnitude of the flux, excluding the direction) in transparent and dark chambers, respectively, as

$$\begin{array}{}\text{(4)}& \mathrm{PQ}={\displaystyle \frac{\mathrm{|}{F}_{{\mathrm{O}}_{\mathrm{2}},\text{light}}-F_{{\mathrm{O}}_{\mathrm{2}},\text{dark}}\mathrm{|}}{\mathrm{|}{F}_{\mathrm{DIC},\text{light}}-F_{\mathrm{DIC},\text{dark}}\mathrm{|}}}\end{array}$$

and

$$\begin{array}{}\text{(5)}& \mathrm{RQ}={\displaystyle \frac{\mathrm{|}{F}_{\mathrm{DIC},\text{dark}}\mathrm{|}}{\mathrm{|}{F}_{{\mathrm{O}}_{\mathrm{2}},\text{dark}}\mathrm{|}}}.\end{array}$$

Due to issues with the dark incubations in the natural meadow, the RQ from this site was instead calculated as the average of the three other sites.

We computed EC fluxes from the high-frequency time series following a multiple-step protocol described in Attard et al. (2019). In short, we bin-averaged the time series to 8 Hz, extracted fluxes for consecutive 15 min time windows using linear detrending (Mcginnis et al., 2014) and corrected fluxes for oxygen storage within the canopy (Rheuban et al., 2014b). Subsequently, we bin-averaged 15 min fluxes to 1 h for interpretation. We defined F_light and F_dark based on when incident PAR was above or below 1 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, respectively. All sites experienced 19 light hours and 5 dark hours on average. We calculated daily metabolic parameters of gross primary productivity (GPP) as

$$\begin{array}{}\text{(6)}& \mathrm{GPP}\left(\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}\right)=(F_{\text{light}}+\mathrm{|}{F}_{\text{dark}}\mathrm{|})\times t_{\text{day}},\end{array}$$

where t_dayis the number of light hours. We calculated community respiration (CR) as

$$\begin{array}{}\text{(7)}& \mathrm{CR}\left(\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}\right)=F_{\text{dark}}\times \mathrm{24}\end{array}$$

and net community productivity (NCP) as

$$\begin{array}{}\text{(8)}& \mathrm{NCP}\left(\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}\right)=\mathrm{GPP}-\mathrm{|}\mathrm{CR}\mathrm{|}.\end{array}$$

We converted oxygen-based daily metabolic fluxes to DIC fluxes by multiplying F_light and F_dark with our empirically derived PQ and RQ, respectively:

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{light}\mathrm{\_}\mathrm{DIC}}\left(\mathrm{mmol}\mathrm{DIC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}\right)=F_{\text{light}}\times {\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{‾}}{\mathrm{PQ}}}}\text{(10)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{dark}\mathrm{\_}\mathrm{DIC}}\left(\mathrm{mmol}\mathrm{DIC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}\right)=F_{\text{dark}}\times -\stackrel{\mathrm{‾}}{\mathrm{RQ}}.\end{array}$$

We then recalculated daily metabolic DIC fluxes GPP_DIC, CR_DIC and NCP_DIC ($\mathrm{mmol}\mathrm{DIC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) using Eqs. (6)–(8). Due to lack of information on the temporal variability in PQ and RQ, we only interpret DIC fluxes on a daily basis.

2.4.2 Biodiversity

We evaluated biodiversity both from a taxonomic and a functional perspective. For taxonomic diversity, we used the “vegan” package in R (Oksanen et al., 2019) to compute Shannon diversity (H^′) and Pielou's evenness component (J^′). H^′ was converted to effective numbers (H_eff = exp(H^′)) to make it linear and scale to species richness (Jost, 2006). For functional diversity, we first assigned functional traits to each species based on the existing literature (Österling and Pihl, 2001; Törnroos and Bonsdorff, 2012; Queirós et al., 2013; Riera et al., 2020; Remy et al., 2021; Kindeberg et al., 2022) and the databases Biological Traits Information Catalogue (MarLIN, 2023) and Polytraits (Faulwetter et al., 2014). The selection of functional traits was based on direct connections to carbon cycling including feeding mode, bioturbation mode and whether the species is calcifying. Indirect, general traits included movement mode, living habit and environmental position. This selection process resulted in 25 trait modalities from which we constructed a traits-by-species matrix assigning each species to specific trait modalities (refer to Table S4 in the Supplement). Species can exhibit multiple trait modalities, depending on life history and environmental conditions. To address this, and to avoid a disproportionally large influence by generalist species on functional diversity, we used fuzzy coding (Chevenet et al., 1994), whereby species comprising multiple trait modalities were assigned a score between 0 (no association) and 3 (full association), with the total sum of each trait always being 3. Based on this matrix, we calculated community-weighted means (CWM) of trait values and several multivariate components of functional diversity including functional richness (FRic), functional evenness (FEve) and Rao's quadratic entropy (RaoQ). These calculations were performed using the “FD” package in R (Laliberté and Legendre, 2010), and further detailed information on these multivariate components and their taxonomic analogs can be found in Mason et al. (2005) and Villéger et al. (2008). As with H^′, the functional diversity index RaoQ was transformed to effective numbers as FD_eff = $\mathrm{1}/(\mathrm{1}-\mathrm{RaoQ})$.

We measured biomass divided into classes. We obtained wet weight (g) after blotting each specimen on a tissue for 2 s and dry weight (g) after drying at 60 °C for 24 h. Regrettably, due to a computer malfunction, the class division per sample was lost for infauna samples, and only the total pooled biomass per site is available for this group. We combusted pooled samples at 520 °C for 4 h to obtain ash-free dry weight (AFDW; g m⁻²).

2.4.3 Light-use efficiency

We evaluated the relationship between irradiance (PAR) and gross primary productivity (GPP) using a hyperbolic tangent function (Jassby and Platt, 1976; Platt et al., 1980):

$$\begin{array}{}\text{(11)}& \mathrm{GPP}=P_{\mathrm{m}}\times \text{tanh}\left({\displaystyle \frac{\mathit{\alpha }\mathrm{PAR}}{P_{\mathrm{m}}}}\right),\end{array}$$

where P_m is the maximum oxygen flux of gross primary productivity ($\mathrm{mmol}{\mathrm{O}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$), α is the quasi-linear initial slope of the curve and PAR is seabed irradiance as photosynthetic active radiation ($\mathrm{\mu }\mathrm{mol}\mathrm{photons}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$). We performed curve-fitting in OriginPro 2022 using a Levenberg–Marquardt iteration algorithm, and we scaled the standard error of the fitting parameters with the square root of the reduced chi-squared statistic (Attard and Glud 2020).

To examine these relationships further, we calculated the light-use efficiency (LUE) at each site, which indicates the efficiency with which absorbed PAR is converted to primary production, as

$$\begin{array}{}\text{(12)}& \mathrm{LUE}={\displaystyle \frac{\mathrm{GPP}}{\mathrm{PAR}\times f\mathrm{APAR}}},\end{array}$$

where fAPAR is the fraction of absorbed irradiance calculated from the difference between incident and reflected PAR as measured by the upward- and downward-facing PAR sensors (see above). Including fAPAR in the calculation of LUE thereby accounts for any differences owing to meadow characteristics, such as the higher three-dimensional meadow complexity (higher fAPAR) relative to bare sediment (lower fAPAR), and captures the diel differences in seabed reflectance and absorption (Attard and Glud, 2020).

2.4.4 Statistical models

To test the effect of differences in abiotic factors between deployments, and thereby validate the use of the chronosequence approach, we employed linear mixed-effects models (package “lme4” in R (Bates et al., 2015)), testing the effect and interaction of abiotic variables on absolute values of hourly oxygen fluxes ($\mathrm{|}{F}_{{\mathrm{O}}_{\mathrm{2}}}\mathrm{|}$). We used model selection (based on the Akaike information criterion (AIC)) to select the best model, which included sea surface temperature, flow velocity, PAR and turbidity as fixed effects and site as a random factor. We used type-III ANOVA for significance testing of fixed effects and likelihood-ratio tests (LRTs) for the random effect. Assumptions of models were tested using the “performance” package in R (Lüdecke et al., 2020). We assessed differences in oxygen fluxes calculated by the EC and BC using a two-sample t-test and compared oxygen and DIC fluxes in the benthic chambers using linear regression analyses. We tested site differences in biodiversity and sediment parameters using multiple one-way ANOVAs and visually reviewed multivariate community composition using non-metric multidimensional scaling (NMDS) and principal component analyses (PCAs). We set the significance level to α = 0.05 for all statistical tests and performed all analyses in R, version 4.2.3 (Rcoreteam, 2023).

Figure 2 Time series of flow, oxygen and light. Time series of (a) flow velocity, (b) ambient dissolved oxygen and (c) hourly oxygen flux overlaid on photosynthetic active radiation (PAR) in yellow.

2.4.5 Carbon budget

We constructed a carbon budget of daily inorganic carbon fluxes and pools of organic carbon. We based the sediment carbon pool on the POC stock of the top 12 cm of sediment, whereas we inferred seagrass aboveground and belowground carbon from dry weight (DW) and a global average carbon content for Z. marina of 34 % DW (Duarte, 1990). We estimated macroalgal carbon content based on DW and species-specific carbon content of the dominant red and brown algae reported from the area, which ranged from 29.1 %–39.9 % DW (Olsson et al., 2020). For benthic fauna, we converted ash-free dry weight (AFDW) to carbon, assuming a 50 % carbon content (Wijsman et al., 1999; Rodil et al., 2021). We converted organic carbon pools to moles, and they are reported as mol C m⁻². Lastly, we calculated the total pool of organic carbon for each site as the sum of all pool means. We calculated the total propagated uncertainty (SE_total) as

$$\begin{array}{}\text{(13)}& {\mathrm{SE}}_{\mathrm{total}}=\sqrt{{\mathit{\sigma }}_{\text{sediment}}^{\mathrm{2}}+{\mathit{\sigma }}_{\text{AG}}^{\mathrm{2}}+{\mathit{\sigma }}_{\text{BG}}^{\mathrm{2}}}+{\mathit{\sigma }}_{\text{algae}}^{\mathrm{2}}+{\mathit{\sigma }}_{\text{fauna}}^{\mathrm{2}}/\sqrt{n},\end{array}$$

where σ is the standard deviation of each pool mean and n is the number of pools. AG and BG are eelgrass aboveground and belowground biomass, respectively.

3 Results

3.1 Environmental conditions

The weather was sunny and dry during all field deployments with only two minor rain events in between (Fig. S2 in the Supplement). Sea surface temperature ranged from 17.10–19.98 °C, driven mainly by the diel light cycle. Salinity ranged from 24.7 to 28.9 but remained constant (± 0.1) during each individual deployment. Photosynthetic active radiation (PAR) at the seabed was similar between sites and deployments and reached a highest value of 728 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ (Figs. 2 and S2). Flow velocities ranged from 0.9 to 21 cm s⁻¹, averaging 5.6 ± 3.4 cm s⁻¹ across all sites (Figs. 2 and S2).

Ambient seawater chemistry was largely similar between sites, although there was a higher background salinity, TA, DIC and DIN at the 3-year and Nat sites, which were sampled after a weather front passed by, likely exchanging some of the bay water with off-shore fjord water (Table S3; Fig. S2 in the Supplement). The average DO during EC deployments was highest in Bare and lowest in 3-year, averaging (mean ± sd) 302.9 ± 21.8 and 260 ± 37.3 µM, respectively. Turbidity was generally low but increased markedly at the Nat site, following a minor rain event prior to deployment (Fig. S2; Table S2 in the Supplement). However, differences in turbidity did not have any detectable effects on seabed PAR (Fig. S2; Table S2 in the Supplement).

3.2 Hourly oxygen fluxes

Hourly O₂ fluxes followed the diel light cycle and increased both in magnitude and variability going from bare sediments to increasing age of the restored seagrass (Fig. 2).

The largest hourly oxygen fluxes typically occurred in the afternoon, with the highest recorded between 14:30–15:30 CET at the 3-year site (8.96 ± 1.44 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$). The largest oxygen uptake rates were generally observed during hours following midnight, with the most negative hourly flux recorded between 23:30–00:30 CET at Nat (−9.08 ± 5.62 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$).

Flow velocity was on average significantly higher in the 3-year compared to Bare and significantly higher in the 7-year and Nat meadow compared to the 3-year (Table S2 in the Supplement). Although there was a general positive linear relationship between flow velocity and absolute oxygen flux across all deployments, the higher flow velocities in 7-year and Nat generally occurred during short time periods at night and did not correspond to consistent increases in absolute oxygen fluxes for those sites (Fig. S3 in the Supplement). Consequently, site R² values were low, ranging from < 0.001–0.20 (Fig. S3 in the Supplement). Further analysis through linear mixed-effects modeling indicated that, while temperature, PAR, turbidity and flow velocity explained 20 % of the variation in hourly $\mathrm{|}{F}_{{\mathrm{O}}_{\mathrm{2}}}\mathrm{|}$ across all sites, the random effect site was highly significant (LRT = 20.9, p < 0.001) suggesting that some other feature, not included in the model, contributed to the observed differences in oxygen fluxes between sites (Table S5 in the Supplement).

Figure 3 Fluxes and relationships of oxygen, carbon and productivity dynamics. (a) Daily oxygen fluxes as gross primary productivity (GPP), community respiration (CR) and net community productivity (NCP). (b) Linear regression of daily oxygen-based GPP and CR. (c) Daily dissolved inorganic carbon (DIC) fluxes based on eddy covariance fluxes converted using the photosynthetic quotient (PQ) and the respiratory quotient (RQ). (d) Linear regression of oxygen and dissolved inorganic carbon (DIC) hourly flux measured in the benthic chambers used to calculate PQ and RQ. The dashed black line indicates slope equal to 1, and dashed gray lines are the 95 % confidence interval of the fitted slope.

3.3 Daily integrated metabolism

Daily metabolic oxygen fluxes (GPP, CR) as measured by the EC were lowest in the bare sediments and increased with meadow age (Fig. 3a). GPP and CR were tightly coupled, but |CR| was always higher than GPP, amounting to an average GPP:CR ratio of 0.81 (Fig. 3b). Accordingly, we observed net heterotrophy at all sites (NCP < 0; Fig. 3a and b). Oxygen-based NCP decreased 3-fold between the bare and the youngest restored meadow (−5 to −15 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), with a further 40 % decrease in the 7-year-old meadow (−21 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$).

Oxygen fluxes measured by the EC and BC were not significantly different from each other (two-sample t-test: p = 0.69), although there was a tendency to overestimate oxygen fluxes in BC relative to EC by 0.7–4.0 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$. Oxygen and DIC fluxes in the benthic chambers were highly correlated across all incubations (Fig. 3d), irrespective of differences in light conditions. The photosynthetic quotient (PQ) was always less than unity, averaging 0.46 ± 0.10 across the four sites, whereas the respiratory quotient (RQ) averaged 0.93 ± 0.25. Site-specific RQ revealed high variability between sites, ranging from 0.65–1.13.

Estimated DIC fluxes mirrored those of O₂, and the benthic DIC net efflux (NCP_DIC) increased as a function of meadow age from 6 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in the bare sediments to 28 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in the natural meadow, thus confirming the net heterotrophic status of the meadows as determined using oxygen fluxes (Fig. 3b and c).

Table 1 Eelgrass and macroalgal structural diversity. Morphometrics, biomass and diversity across the sites (mean ± SE). “Rep. shoots” represents reproductive shoots with seed spathes present. AG and BG are aboveground and belowground eelgrass biomass, respectively, as captured by sediment cores. Macroalgal biomass represents macroalgae collected from eelgrass canopy samples. Maximum number of species refers to the count of macroalgal species found in a sample. Relative proportion indicates macroalgal biomass relative to total macrophyte biomass. Species richness, diversity (H_eff) and evenness (J^′) refer to macrophytes, including macroalgae and eelgrass. An asterisk indicates statistical significance (p < 0.05).

3.4 Structural and functional diversity

3.4.1 Meadow properties

All three eelgrass sites were characterized by high spatial heterogeneity within each meadow (Table 1). No significant differences were observed in seagrass morphometry such as shoot density, canopy height, or below- or aboveground biomass. The only seagrass parameter that differed between the meadows was the number of reproductive shoots containing seeds, which was significantly higher in the natural meadow (p = 0.004). However, the abundance and biomass of other macrophyte species such as brown and red macroalgae increased markedly with meadow age. In the 3-year meadow, only a small specimen of the brown algae Spermatochnus paradoxus was found in one sample, whereas in the natural meadow large quantities of up to five different macroalgal species were found. However, due to the large variability in biomass between samples within each site (Table 1), the between-site differences in number of species were not statistically significant (ANOVA, p > 0.05). The composition of macrophyte species became more even with meadow age such that while the 3-year meadow was dominated by Z. marina (∼ 99 % of total macrophyte biomass), the 7-year and natural meadows had a more heterogenous and evenly distributed macrophyte community, where Z. marina on average contributed 90 ± 15 % and 64 ± 32 %, respectively, to total macrophyte biomass (Table 1). As a result, the three-dimensional complexity of the canopy increased with meadow age, driven mainly by large-bodied drifting fucoid species (F. serratus and F. vesiculosus) and red algae (Furcellaria lumbricalis) residing unattached within the canopies.

Benthic microalgae, as inferred from chlorophyll a on the sediment surface, showed the opposite trend and decreased with meadow age, and chlorophyll a was significantly lower in sediments underlying 7-year (0.28 ± 0.03 g m⁻²) and Nat (0.26 ± 0.01 g m⁻²) compared to Bare (0.56 ± 0.07 g m⁻²).

3.4.2 Benthic fauna

We collected a total of 1927 individuals representing 43 taxa. Taxonomic diversity parameters (abundance, number of species, Shannon diversity and evenness) exhibited large within-site variability, highlighting the small-scale (< 10 m) heterogeneity of fauna community structure. These parameters consistently showed higher values in vegetated compared to bare sediments but exhibited variable, but generally non-significant, differences between the eelgrass sites (Fig. 4; Table S6 in the Supplement). The abundance of infauna was highest at the 3-year site, primarily dominated by opportunistic polychaetes (e.g., Capitella capitata). However, despite the high abundance at the 3-year site, this site did not show a corresponding spike in infaunal species richness but was reflected by the lowest evenness among all sites (J^′ = 0.47 ± 0.08). By contrast, in the 7-year, the abundance had decreased by one-third, while species diversity (H_eff) and evenness (J^′) nearly doubled, exhibiting similar values to the natural reference meadow (Table S6 in the Supplement). Functional trait metrics revealed that both the functional group richness (FGR) and functional diversity (FD_eff) were significantly higher in the 7-year and Nat compared to the 3-year and Bare sites which exhibited similar values (Table S6 in the Supplement). Functional richness (FRic) was notably low in the bare sediments (0.06 ± 0.05) and tended to increase with meadow age, peaking with the highest mean value in the natural meadow (0.53 ± 0.11). However, due to high within-site variability, FRic did not show statistical differences between sites (Fig. 4c).

When separately targeting the meadows for epifauna, we found that they were species-rich and highly diverse, with number of species ranging from 15–18 and H_eff from 7.4–10.9. However, neither taxonomic nor functional diversity metrics exhibited any significant differences between the meadows, although there were some increasing trends particularly in functional evenness (FEve), which was highest in Nat and lowest in 3-year (Table S6 in the Supplement). Epifaunal biomass increased on average 3-fold at Nat compared to the two restored meadows, displaying the highest within-site variability (15.89 ± 10.48 g m⁻²), primarily driven by gastropods.

Figure 4 Biodiversity patterns in benthic fauna. Panels to the left show species richness (a, b) and functional richness (c, d) of infauna samples (a, c) and epifauna samples (b, d). The large panel to the right (e) shows relative abundance of different classes of all fauna combined (infauna plus epifauna).

Community composition partly shifted as the meadow grew, whereby bare sediments and the youngest restored meadow were dominated by polychaetes, whereas more epifaunal species such as bivalves and crustaceans were found in older meadows (Fig. 4d). Absolute abundances and biomass supported these observations, with bryozoans and gastropods contributing to higher biomass at Nat relative to Bare. However, multivariate analysis of the different communities indicated an overlap in community composition (Fig. S4 in the Supplement).

Our analysis of functional traits highlighted the prevalence of certain bioturbation modes in relation to meadow age. For instance, community-weighted means (CWMs) of biodiffusers displayed a linear increase with meadow age and were significantly higher in the natural meadow and the 7-year compared to the 3-year (F_3,20 = 8.4; p < 0.001). Surficial modifiers among infauna were higher in eelgrass compared to bare sediment, peaking in the oldest restored meadow at a CWM of 0.29 ± 0.10.

Table 2 Sediment properties across sites. Organic matter (OM), particulate organic carbon (POC), particulate inorganic carbon (PIC), total nitrogen (TN) and dry bulk density (DBD) of the top 2 cm of sediment. POC_stock is the depth-integrated carbon stock over 0–12 cm sediment depth. Values are mean ± SE, n = 3 per site.

3.5 Sediment carbon stocks

The sediment within Gåsö bay eelgrass meadows has previously been reported as silty sand, with a median grain size (D₅₀) of the surface sediment between 140–170 µm and a silt–clay content of 26 %–35 % (Infantes et al., 2022; Dr. Martin Dahl, personal communication, 2023). Concentrations of sediment OM, POC and TN were not significantly different between sites (p < 0.05) and did not display any consistent increases or decreases with meadow age (Table 2). However, when integrating the POC density across the top 12 cm, the highest POC stock was found in one natural-meadow core (1529 g m⁻²) and the lowest in a bare-sediment core (209 g m⁻²). However, due to large within-site variability, the sites were not significantly different from each other (F_3,8 = 1.52; p = 0.28; Table 2). Meadow age did not show a clear trend, as demonstrated by the 3-year site, which had an average carbon stock 32 % larger than the 7-year, while the 7-year site was more similar to the bare-sediment site (Table 2). Depth profiles of POC concentration and density down to 20 cm revealed near-constant values down to between 12 and 16 cm, where an increase was observed (Fig. S5 in the Supplement). Natural eelgrass had the highest average POC profile, but average values were highly skewed by one core replicate displaying POC density 4 times higher than the other two replicates within the site.

Figure 5 Light-use efficiency and productivity relationships. Panels (a–c) show different components of light-use efficiency (LUE) as a function of hour of the day: panel (a) shows incident photosynthetic active radiation (PAR); panel (b) shows the fraction of absorbed PAR (fAPAR), and panel (c) shows gross primary productivity (GPP) as a function of time of day. Panel (d) shows the relationship between oxygen flux and PAR as hyperbolic tangent curves estimated for each site. α denotes the quasi-linear initial slope.

3.6 Light-use efficiency

All meadows experienced similar incident light conditions (Fig. 5a). The fraction of absorbed light (fAPAR) was always higher in eelgrass (∼ 97 %) compared to bare sediments (∼ 94 %) but did not differ between eelgrass sites on a daily basis (Fig. 5b). Hourly GPP tracked PAR with a clear hysteresis effect evident in the 7-year and natural meadow but to a lesser extent at the bare site (Figs. 5c and S6 in the Supplement). P–I relationships were best explained by the hyperbolic tangent function yielding R² between 0.45 and 0.74. The irradiance needed for photosynthesis to balance respiration (I_k) was almost 4 times higher at the Bare site compared to the 7-year site, equaling 380 and 97 $\mathrm{\mu }\mathrm{mol}\mathrm{photons}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, respectively (Table S7 in the Supplement; Fig. 5d). Estimated light-use efficiency (LUE) was lowest at Bare (0.001 O₂ photon⁻¹) and increased with meadow age to 0.004 and 0.005 O₂ photon⁻¹ at 3-year and 7-year, respectively. The highest daily LUE was observed at Nat (0.007 O₂ photon⁻¹), coincident with the highest number of macrophyte species and the most diverse community structure (Fig. 6).

Table 3 Daily metabolism as a function of meadow age. Curve-fitting of daily metabolism parameters GPP, CR and NCP to meadow age (SiteAge) converted to logarithmic scale (log₁₀(x+1)). SiteAge was defined as 0 years for the site Bare and 13 years for the natural meadow.

Similar to LUE, GPP and |CR| displayed a positive linear relationship with the number of macrophyte species. There was also a positive trend between these parameters and macrophyte Shannon diversity (H_eff) and the proportion of macroalgal biomass relative to eelgrass biomass (Fig. 6). As such, the model that best explained changes in daily benthic metabolism across the four different stages of seagrass development was a logarithmic model (Table 3).

Figure 6 Biodiversity and productivity relationship. (a) Light-use efficiency (LUE) as a function of the relative biomass of macroalgae to eelgrass ($R_{\mathrm{adj}}^{\mathrm{2}}$ = 0.70). (b) Gross primary productivity (GPP) as a function of macrophyte Shannon diversity index ($R_{\mathrm{adj}}^{\mathrm{2}}$ = 0.46). Black lines represent best fit (log_e(x+1)). Note that the Bare site was not quantitatively sampled for macroalgal proportions or macrophyte diversity and was a posteriori set to 0 % and 1, respectively, for curve-fitting.

Table 4 Carbon pools. Pools of particulate organic carbon (mean ± SE, mol C m⁻²) in the different components of the benthic habitats. AG and BG are above- and belowground eelgrass biomass, respectively. Fauna is total fauna (infauna plus epifauna).

Note: n.d. – not determined.

3.7 Carbon pools

Converting seagrass community components to carbon illustrates the pools of carbon available for export, remineralization or burial. Notably, total carbon pools were higher in eelgrass relative to bare sediment but were similar between restored and natural seagrass (Table 4). Sediment POC stocks were the largest carbon pools followed by eelgrass biomass, which contributed on average 11, 21 and 7 % to the total carbon pool in the 3-year, 7-year and natural meadow, respectively (Table 4).

4 Discussion

We present a comprehensive dataset documenting post-restoration seagrass development that captures several aspects of seagrass metabolism. This dataset enables us to investigate the role of biodiversity and different components of a seagrass ecosystem in carbon cycling. We show that (i) community-integrated photosynthetic (GPP) and respiratory (CR) fluxes increase as a function of meadow age (Fig. 3); (ii) daily |CR| increased more relative to GPP, resulting in net heterotrophy (NCP < 0) on diel timescales during summer; (iii) diversity and biomass of macrophytes other than the restored seagrass could be driving higher primary productivity through increased light-use efficiency (Fig. 5); (iv) faunal communities recover rapidly and attain species richness and functional richness comparable to natural meadows within 7 years since restoration (Fig. 4); and (v) surficial (0–12 cm) sediment carbon stocks are large but are not significantly affected by the presence of seagrass in this sheltered bay.

Based on the above results, we postulate that, while increased macrophyte diversity enhances both GPP and CR, the additional CR stemming from benthic fauna communities, together with labile organic matter input, pushes diverse seagrass meadows toward net heterotrophy during summer. This highlights potential tradeoffs between climate change mitigation and biodiversity conservation as incentives for seagrass restoration. Below we discuss four primary lines of evidence to support this postulation.

4.1 Metabolic fluxes scale to meadow development

We found large daily fluxes of GPP- and CR-derived O₂ and DIC that increased as the system developed from bare sediments to a mature meadow (Table 3; Fig. 3). Our values (mean ± SE) of GPP, CR and NCP across the three seagrass sites were 80 ± 11, −99 ± 13 and −19 ± 2 $\mathrm{mmol}{\mathrm{O}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, respectively, which is relatively low when comparing to the global average GPP, CR and NCP estimated for temperate seagrasses of 166 ± 14, −130 ± 11 and 34 ± 8 $\mathrm{mmol}{\mathrm{O}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, respectively (Duarte et al., 2010). However, it should be noted that discrepancies owing to methodological differences are difficult to account for. An updated assessment of seagrass NCP in temperate areas reported an average of 29 ± 79 $\mathrm{mmol}{\mathrm{O}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, although the study which covered 187 seagrass metabolism studies found that merely 68 % reported net autotrophy (Ward et al., 2022). Accordingly, the notion that seagrass habitats are strongly net autotrophic is being increasingly contested as methods continue to improve. In all our sampled sites, GPP was lower than |CR|, resulting in net heterotrophy (negative NCP), independently established by both EC oxygen fluxes and benthic chamber DIC and oxygen fluxes. Several recent studies have reported instances of sustained net heterotrophy across multiple seagrass species and environments (e.g., Barron et al., 2006; Rheuban et al., 2014a, b; van Dam et al., 2019; Berger et al., 2020; Attard et al., 2019; Berg et al., 2022, Ward et al., 2022). For instance, a recent study of Z. marina using EC reported GPP and CR values similar to our natural meadow (95 and 94 $\mathrm{mmol}{\mathrm{O}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, respectively), resulting in a near-balanced metabolic state across 11 years of monitoring (Berger et al., 2020). The authors reported a generally balanced metabolic state on monthly timescales, but, following a temperature-driven dieback event that diminished seagrass shoot density, GPP and |CR| decreased by 55 % and 48 %, respectively. This shifted the meadow to net heterotrophy during summer (NCP = −26 ± 15 $\mathrm{mmol}{\mathrm{O}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). In the following years, the gradual increase in seagrass shoot density increased primarily GPP, showing clear signs of seagrass recovery (Berger et al., 2020).

Although GPP often substantially increases during summer in temperate seagrass meadows, so does CR to a similar extent (Ward et al., 2022). Consequently, despite large seasonal variability in photosynthesis and respiration, the metabolic state is often relatively stable on an annual basis, granted there are no major ecosystem regime shifts. While interannual variability of NCP has been related to seagrass die-off and recovery episodes (Berger et al., 2020), seagrass phenology linked to abiotic factors such as temperature and light regimes typically dictates the metabolic state on intra-annual timescales (e.g., Champenois and Borges, 2012; Rheuban et al., 2014a, b). Here, we show that biotic components other than the seagrass itself can contribute to both the magnitude of and variability in metabolic fluxes. Irrespective of traditional seagrass metrics such as seagrass shoot density and biomass, GPP and |CR| consistently increased in magnitude with meadow age, which in turn corresponded to higher autotrophic diversity and macroalgal biomass. Further research should address whether these relationships are consistent across seasons and what role differing macrophyte phenologies play.

The chronosequence approach employed in this study utilizes the unique opportunity of assessing contrasting restored seagrass habitats of different ages that exist within a close distance from each other (Fig. 1). This enables comparisons between near-identical geomorphology, bathymetry, hydrodynamics and seawater characteristics. However, due to logistical limitations, we were unable to measure all four sites simultaneously, leading to a temporal mismatch of these comparisons. Consequently, this introduces the risk of potential environmental changes between deployments. Importantly, if the change in environmental conditions is conducive to altered benthic metabolism it can influence the comparison along the chronosequence (i.e., between sites). The combined effect of abiotic variables, including PAR, flow velocity, seawater temperature and turbidity, accounted for 20 % of the variation in O₂ fluxes, as measured by the eddy covariance. Noticeably, PAR reaching the seabed did not differ between sites, despite varying levels of turbidity (Figs. 2 and  S2). Salinity was higher at the 3-year and Nat sites compared to the 7-year and Bare sites (Table S2; Fig. S2 in the Supplement). However, due to missing data, we could not evaluate its impacts on oxygen fluxes within the model. That being said, we found no discernable effects on either oxygen or carbon fluxes during our incubations, suggesting that variability in salinity was not a driving factor of metabolism. Flow velocity peaked at the Nat and 7-year sites, but, while there was a positive relationship between $\mathrm{|}{F}_{{\mathrm{O}}_{\mathrm{2}}}\mathrm{|}$ and flow at the Nat and 3-year sites, no such relationship was evident in the 7-year or Bare sites (Fig. S3 in the Supplement). Nonetheless, we cannot decisively rule out the potential role of varying flow velocities in the observed differences in benthic metabolism between sites.

4.2 Carbon and oxygen balance

As part of this study, we present a methodological approach that estimates in situ DIC fluxes under natural hydrodynamic and light conditions. This is obtained by combining the advantages of aquatic eddy covariance with the ability to constrain the marine carbonate system and oxygen dynamics using benthic chambers. Concurrent deployment of these two methods has been utilized in previous coastal studies (Long et al., 2019; Camillini et al., 2021; Polsenaere et al., 2021) but only for comparing oxygen fluxes.

Assessing the in situ relationship between oxygen ($F_{{\mathrm{O}}_{\mathrm{2}}}$) and DIC fluxes (F_DIC) can provide insights into biogeochemical processes and renders reliable estimates of photosynthetic quotients (PQs) and respiratory quotients (RQs). All else equal, photosynthetic and respiratory quotients are governed by the $\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratio of the fixed and respired organic matter present in the system (Champenois and Borges, 2021). However, considering the various sinks and sources of organic matter present in a seagrass meadow and the multitude of processes affecting $F_{{\mathrm{O}}_{\mathrm{2}}}$ and F_DIC differently, this is not very useful. Deviations from the theoretical 1:1 relationship between $F_{{\mathrm{O}}_{\mathrm{2}}}$ and F_DIC (Fig. 3d) are thus ubiquitous in the literature (e.g., Barron et al., 2006; Turk et al., 2015; Trentman et al., 2023). In fact, the slope we observed is identical to what Pinardi et al. (2009) observed in sediments vegetated with the freshwater macrophyte Vallisneria spiralis using sediment cores. Moreover, our relatively low PQs (0.34–0.52) were similar to what Ribaudo et al. (2011) observed in V. spiralis (0.30–0.68) in microcosms. The authors attributed the low PQ to oxygen transport to roots and subsequent radial oxygen loss (ROL) which fuels aerobic respiration, a process well-documented in Z. marina as well (e.g., Jensen et al., 2005; Frederiksen and Glud, 2006; Borum et al., 2007; Jovanovic et al., 2015). Turk et al. (2015) observed PQs ranging from 0.5–2.6 in seagrass (Thalassia testudinum) and found a temporal component to the variability in PQ with lower values in the morning and higher in the evening (Turk et al., 2015). Similar to our study, Ouisse et al. (2014) obtained a PQ and an RQ of 0.42 ± 0.27 and 0.95 ± 0.22, respectively, using in situ benthic chambers in dwarf eelgrass (Z. noltii) meadows across several seasons. The authors hypothesized that the low PQ could also be due to photorespiration in epiphytic algae on the seagrass leaves which can consume more than 3 mol O₂ for every 1 mol DIC used (Ouisse et al., 2014). We observed large quantities of epiphytic microalgae and biofilm on seagrass leaves across all our studied meadows, albeit only as qualitative observations (Kindeberg, T., pers. obs.). However, seagrass epiphytes are abundant in the area where they can exert detrimental effects on seagrass metabolic performance and positive effects on epifauna distribution (Baden et al., 2010; Gullström et al., 2012; Brodersen et al., 2015). It is important to note that inorganic processes (i.e., CaCO₃ production and dissolution), which can have a large influence on PQ and RQ (Champenois and Borges, 2021), are implicitly accounted for in our F_DIC by subtraction of the 0.5F_TA term in Eq. (3).

While we obtained an average RQ close to unity, it was based on a relatively small sample size compared to PQ due to issues with dark incubations, especially in the natural meadow. It is possible that our acclimation (∼ 30 min) or incubation times (3.0 ± 0.1 h) were too short for accurately capturing dark DIC fluxes, as seen in the temporal lag in DIC fluxes relative to O₂ fluxes in a study by Fenchel and Glud (2000) and a lag in O₂ consumption due to the primary producer cellular machinery (Tang and Kristensen, 2007). Nonetheless, without any ancillary data on other biogeochemical processes, we cannot reconcile the sources of our observed PQ and RQ.

4.3 Macrophyte diversity driving light-use efficiency and higher metabolism

Despite the large research field on the relationship between biodiversity and primary productivity (Tilman et al., 2014), light-use efficiency (LUE) is largely understudied in benthic metabolism studies (Attard and Glud, 2020). Studies have hitherto focused mainly on smaller-scale LUE, such as microalgae in microbial mats and corals (Al-Najjar et al., 2010, 2012; Brodersen et al., 2014). We observed a positive relationship between macrophyte diversity and LUE when controlling for biomass, indicating that mixed meadows consisting of both seagrass and macroalgae utilize light resources more efficiently and are more productive compared to monospecific meadows. Importantly, the restored seagrass meadows became more mixed over time as drifting macroalgae inhabited the meadow. These unattached algae are a common feature in the area, often considered a nuisance that can prevent seagrass recovery (Moksnes et al., 2018). Here it seems they also improve overall LUE of the meadow and contribute to larger metabolic fluxes.

Figure 7 Pools and fluxes of carbon. Schematic illustration of a carbon budget including different benthic pools (mol C m⁻²) of particulate organic carbon in (i) fauna biomass, (ii) sediment, (iii) eelgrass aboveground biomass, (iv) eelgrass belowground biomass and (v) macroalgal biomass. Arrows indicate the daily metabolic fluxes of dissolved inorganic carbon ($\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), where blue arrows are gross primary productivity (GPP_DIC), red arrows are community respiration (CR_DIC) and orange arrows are net community productivity (NCP_DIC).

Whether higher LUE is driven by certain species remains unclear, but the change in canopy structure and increasing three-dimensional complexity can positively influence LUE (Zimmerman, 2003; Binzer et al., 2006). Niche complementarity is common in ecological systems (Loreau and Hector, 2001; Hooper et al., 2005), and it is reasonable to believe that, with increased diversity of autotrophs, pigment complementarity can facilitate optimal resource use, especially as brown and red macroalgae are known to have a wide range of photosynthetic pigments (Enríquez et al., 1994). Additionally, the mere presence of multiple growth morphologies may induce self-shading that further increases LUE (Tait et al., 2014). An increase in photosynthetic pathways (e.g., C3 and C4) with higher macrophyte diversity and differing affinity for forms of inorganic nutrients (e.g., ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$) is also expected. Moreover, both Z. marina and fucoid species are known to utilize both CO₂ and ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ for photosynthesis (Binzer et al., 2006). However, the efficiency differs between species (e.g., Larsson and Axelsson, 1999; Invers et al., 2001), and, considering the large spatiotemporal variability in pH and [${\mathrm{HCO}}_{\mathrm{3}}^{-}$] relative to [CO₂] we observed, this could be another reason for the higher LUE at higher species diversity. Studies from macroalgal canopies have found similar relationships between macrophyte canopy complexity and LUE, attributed to niche complementarity where intact assemblages are more efficient and productive than the sum of their parts (Tait and Schiel, 2011; Tait et al., 2014). For instance, a study by Tait and Schiel (2011) using ex situ incubation chambers found that an intact assemblage of seven species had higher net photosynthesis than the sum of all individual species. The authors observed that different species played different roles at different irradiances. For instance, the fucoid species Cystophora torulosa was exceptionally efficient at photosynthesizing at high irradiance and did not show signs of photoinhibition even at PAR >2000 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ (Tait and Schiel, 2011). Tait et al. (2014) studied P–I relationships in macroalgal assemblages and found that when more sub-canopy species were included (up to four), respiration and photosynthesis increased, thus corroborating our observed trends. However, they found that production did not saturate at incident irradiance of 2000 µmol m² s⁻¹, as opposed to less speciose assemblages (two sub-canopy species) that already reached light saturation of net primary production (NPP) at about 600 µmol m² s⁻¹ (Tait et al., 2014). This is somewhat contrary to what we found for GPP, where P–I curves saturated at lower irradiance with higher macrophyte diversity (Figs. 5d and S6). Albeit not specifically addressing canopy structure or diversity, Rheuban et al. (2014a) found that a younger, 5-year-old restored Z. marina meadow was light-saturated, while an older, 11-year-old meadow did not show any signs of light saturation, and this was consistent across seasons.

Whereas it is rather intuitive that a diverse community of primary producers are better at photosynthesizing (i.e., higher GPP), the relationship is as strong with CR. This is likely explained by the tight coupling between GPP and CR stemming from respiration of labile photosynthates (Penhale and Smith, 1977). However, detritus of macroalgae such as Fucus spp. is also more labile than Z. marina, partly due to a lower C:N ratio and a more bioavailable polysaccharide composition (e.g., Kristensen, 1994; Thomson et al., 2020).

4.4 The role of benthic diversity in seagrass metabolism

The fact that most fauna diversity metrics were not significantly different between the natural meadow and the youngest (3-year) meadow implies that benthic diversity recovers quickly. Similar findings have recently been reported from Z. marina restoration projects in Denmark (Steinfurth et al., 2022) and from the very same sites as in this study (Gagnon et al., 2023). In fact, Gagnon et al. (2023) found that both taxonomic and functional diversity recovered within 15 months after restoration, but after 3 months the abundance was already similar to documented abundances in comparable seagrass meadows in the area. The authors partly attributed this to efficient larval dispersal from the adjacent natural meadow within the bay (Gagnon et al., 2023).

It is generally established that higher diversity yields higher productivity in seagrass meadows (Duffy et al., 2017), although the mechanisms behind the relationship are debated (Hooper et al., 2005; Gamfeldt et al., 2015). Based on our results, it seems that high macrophyte and macrofauna diversity positively influence GPP and CR, respectively, although the relationships with fauna are less clear. Aside from direct cellular respiration, many infauna species indirectly modify metabolic fluxes across the sediment–water interface through bioturbation and sediment reworking (Aller and Aller, 1998; Kristensen et al., 2012). A scrutiny of bioturbation and reworking modes revealed that especially biodiffusers and surficial modifiers increased with meadow age, despite highest total abundance of infauna in the youngest meadow (Tables S4 and S6 in the Supplement). It is possible that these functional modes benefited from larger quantities of macroalgal detritus building up on the sediment surface over the years. Thomson et al. (2020) found the lugworm Arenicola marina, an upward conveyor, contributed to a 37 % higher efflux of DIC in sediments containing F. vesiculosus compared to Z. marina. Macroalgal detritus was to a much higher extent respired or consumed compared to seagrass, which instead was buried in anoxic sediment layers by the lugworm (Thomson et al., 2020). Moreover, the role of bioturbation in oxygenating otherwise anoxic sediment can have large ramifications for sediment–water fluxes of DIC and could hence contribute to our observed CR.

4.5 Implications of seagrass restoration on the carbon budget

The observed net heterotrophy during the productive season implies the system relies either on historic production of autochthonous carbon or on trophic subsidies to sustain metabolism. Albeit only covering a brief period within the summer season, our results suggest that the seagrass in this area receives large amounts of allochthonous carbon that is partly turned over and released as DIC. A large influx and sedimentation of allochthonous carbon was shown in a recent study by Dahl et al. (2023) from the same bay. They reported relatively high carbon accumulation rates (0.91 ± 0.06 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) of which 51 % of sediment carbon originated from eelgrass productivity and 38 % from macroalgae (Dahl et al., 2023). Assuming this rate is constant throughout the year (0.0025 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), this accumulation rate is about 1 order of magnitude lower than our measured summer NCP_DIC, implying that the majority of imported carbon is rapidly remineralized or assimilated by secondary producers (Fig. 7).

Our estimated budget of all carbon pools illustrates that, whereas sediment stocks are the dominant pools, organic carbon is built up in living biomass following restoration (Table 4; Fig 7). Eelgrass and macroalgal biomass in the natural meadow made up 58 % and 27 % of all biomass, respectively, which is on the same order as the relative proportion of sediment POC sources found in Dahl et al. (2023). Accumulation of sediment carbon and production of living biomass can be decoupled on longer timescales, although trophic subsidies (i.e., external inputs) may be required to sustain both (Cebrian et al., 1997; Duarte et al., 2010; Huang et al., 2015). Notably, total fauna biomass also increased with meadow age, despite varying differences in abundance (Figs. 4 and 7; Table S6 in the Supplement).

While we are able to resolve the dominant carbon pools and metabolic fluxes, the import and export of organic carbon over seasonal timescales is required to reconcile the annual carbon cycling at this site. Nevertheless, it is reasonable to infer that NCP and carbon sequestration in these seagrass systems are sustained by lateral import of allochthonous organic carbon.

5 Conclusion

Planting seagrass initiates a profound transformation of the benthic environment that influences biodiversity and carbon cycling. Throughout our field study, we found that, while fauna diversity developed in an anticipated successional pattern, the metabolic fluxes and net release of DIC were always higher in seagrass. These fluxes increased with meadow age, and we observed increasing gross primary productivity and respiration as the seagrass grew and drifting algae and benthic fauna colonized. Collectively, our findings suggest a scenario where higher macrophyte and fauna diversity drives high primary productivity and respiration, respectively. Together with ample input of sestonic organic matter to this sheltered bay, these productive meadows act as effective bioreactors of organic carbon on diel timescales during summer, as evidenced by the net heterotrophic state and net efflux of DIC. These results highlight the intricate connections between carbon cycling and biodiversity that should be taken into consideration when restoring seagrass, especially in sheltered environments with large input of external organic matter. However, identifying the individual mechanisms and constraining the relative importance of fauna and flora diversity for benthic carbon fluxes remains a challenging task and should be a focal point in future research.