Seasonal ecosystem metabolism across shallow benthic habitats measured by aquatic eddy covariance

Shallow benthic habitats are hotspots for carbon cycling and energy flow, but metabolism (primary production and respiration) dynamics and habitat‐specific differences remain poorly understood. We investigated daily, seasonal, and annual metabolism in six key benthic habitats in the Baltic Sea using ~ 2900 h of in situ aquatic eddy covariance oxygen flux measurements. Rocky substrates had the highest metabolism rates. Habitat‐specific annual primary production per m2 was in the order Fucus vesiculosus canopy > Mytilus trossulus reef > Zostera marina canopy > mixed macrophytes canopy > sands, whereas respiration was in the order M. trossulus > F. vesiculosus > Z. marina > mixed macrophytes > sands > aphotic sediments. Winter metabolism contributed 22–31% of annual rates. Spatial upscaling revealed that benthic habitats drive > 90% of ecosystem metabolism in waters ≤5 m depth, highlighting their central role in carbon and nutrient cycling in shallow waters.

Land-ocean transition zones are biodiversity and productivity hotspots. In addition to phytoplankton production occurring in the water column, the submerged coastal landscapes are colonized by a mosaic of distinct communities of seafloor (benthic) vegetation that play a significant role in the oceanic C cycle (Duarte 2017). Large emergent canopies of seagrass and macroalgae, along with biofilms of less conspicuous microscopic algae (microphytobenthos), are estimated to cover a third of the world's continental shelf (Gattuso et al. 2006). Vegetated habitats can synthesize organic matter at remarkably high rates (> 100 mmol C m −2 h −1 ) (Bordeyne et al. 2017), forming extensive canopies and globally significant C stocks (Fourqurean et al. 2012). On an annual basis, vegetated canopies typically generate a surplus of organic C production (net autotrophy), and a substantial proportion (~40%) of their fixed C can be exported (Duarte and Cebrián 1996), enhancing secondary production in terrestrial and aquatic realms through spatial subsidies (Polis et al. 1997).
The important role of shallow benthic habitats for C cycling and energy flow has been known for decades (Mann 1973;Smith 1981;MacIntyre et al. 1996), but we remain widely uninformed about how metabolism (primary productivity and respiration) varies across the coastal seafloor. Obtaining overall estimates of coastal productivity across multiple habitats that span highly contrasting biogenic communities, from rocky algal canopies to soft-sediment environments, has been challenging due to a lack of a common methodology allowing for comparable metabolism estimates. This has limited our understanding of the contribution of biogenic habitats to coastal productivity. Importantly, habitats with contrasting structural biodiversity elements may have very different productivity patterns, despite occurring in close proximity (Eyre and Maher 2011). The first step toward being able to understand productivity in heterogeneous coastal areas is therefore to perform comparative measures of metabolism across the mosaic of habitats that characterize the coastal seascape, and then to extend these measurements over the year to resolve their seasonal dynamics.
In this study, we explore benthic metabolism within the shallow eutrophic waters of the Baltic Sea. We investigate metabolism dynamics in six contrasting shallow subtidal habitats over a year, each representing a major habitat type of the nearshore Baltic (canopies of Fucus vesiculosus, Zostera marina, and mixed macrophytes, a Mytilus trossulus reef, photic sandy sediments, and a deeper aphotic muddy basin). We used the aquatic eddy covariance (AEC) oxygen flux method (Berg et al. 2003), a recent technological development that allows investigating habitat-scale (10s of m 2 ) metabolic rates and their drivers noninvasively at a high temporal resolution (1 h or less) (Berg et al. 2007;Berg et al. 2017). Based on this extensive data set, which also incorporates data on species' abundance and biomass, we elucidate the potential role of the different habitats in C and nutrient cycling, sequestration, and energy transfer within a highly heterogeneous coastal setting.

Study location and sampling
Our six study sites were located nearby the Tvärminne Zoological Station on the Hanko Peninsula in SW Finland (59.844 N, 23.249 E), and consisted of five photic habitats and one deeper aphotic site (Fig. 1), namely: bare sediments (photic: 3 m depth, sandy sediment, and aphotic: 34 m depth, muddy sediment), a mixed macrophyte canopy (3 m depth, sandy), a seagrass canopy (Z. marina, 3.5 m depth, sandy), a macroalgal canopy (F. vesiculosus, 2.0 m depth, rocky), and a mussel reef (M. trossulus, 5 m depth, rocky). In total, flux measurements were performed on 40 occasions between May 2016 and December 2017.

Biological sampling
We developed a sampling protocol to quantify dominant features of seafloor biodiversity within the eddy covariance flux footprint . The extracted samples were used to estimate the standing biomass (g C m −2 ) of the main phototrophic components and macrofauna of each habitat using conversion ratios for dry weight (macrophytes), ash-free dry weight (macrofauna), or chlorophyll a (microphytobenthos) (Supporting Information Section S1).

Benthic oxygen fluxes
Benthic metabolism rates were computed from O 2 fluxes quantified in situ using the AEC technique (Berg et al. 2003). Our AEC systems have two fast-response O 2 microsensors (T 90 ≤ 0.3 s, low stirring sensitivity <1%) (Revsbech 1989) for redundancy and comparison. These sensors are interfaced with an acoustic Doppler velocimeter (Nortek; Fig. 1f) (McGinnis et al. 2011). This setup is capable of resolving very small fluxes down to 1 mmol O 2 m −2 d −1 or less by capturing the entire range of flux-contributing turbulent eddies under typical hydrodynamic conditions within the benthic boundary layer (Berg et al. 2009). Additional sensors located on the AEC frame logged transmitted (seabed) photosynthetically active radiation (PAR) , water temperature and salinity (U24 HOBO), and dissolved O 2 concentration (U26 HOBO) every 15 min.
We followed standard guidelines for instrument setup and data processing (Lorrai et al. 2010;Donis et al. 2015) (Supporting Information Section S2). Computing daily metabolism rates from the 15 min fluxes followed a three-step process (Fig. 2). For each deployment, multiple days of quality-checked O 2 fluxes and PAR were averaged by the time of day to produce a single continuous 24 h time series. The minimum number of flux days per deployment was 1.0, the maximum was 4.9, and the mean was 3.0 (Supporting Information Table S1). The fluxes were then separated into daytime fluxes (FLUX day ; when PAR>0.0 μmol m −2 s −1 ) and nighttime fluxes (FLUX night ; when PAR < 0.0 μmol m −2 s −1 ), and the PAR time series was used to determine the number of daylight hours (h day ). Daily gross primary productivity (GPP, mmol O 2 m −2 d −1 ) was computed as GPP = (FLUX day + |FLUX night |) * h day . Respiration rates (R, mmol O 2 m −2 d −1 ) were calculated as R = |FLUX night | * 24. Net ecosystem metabolism (NEM) was computed as the difference between daily GPP and R. Positive NEM indicates surplus organic C and O 2 production (autotrophy); negative NEM indicates net heterotrophy. Daily GPP, R, and NEM were converted to C equivalents (g C m −2 d −1 ) assuming a quotient of 1.0 for both primary productivity and respiration. Annual rates were estimated as the numerical integration (mathematical area) of daily GPP, R, and NEM over Julian days using linear interpolation to gap-fill between discrete measurements. Linear regression was performed on the accumulated signals for GPP, R, and NEM over the year to investigate whether the slope of the integration was significantly different from zero (95% confidence level for fitting parameters). Numerical integration and statistical analyses were performed in OriginPro 2018 (OriginLab Corporation).

Sampling sites
Bottom-water temperature ranged from 0.2 C in March to 16.5 C in August, and salinity was between 5.2 and 6.7 throughout the year. Daily PAR ranged from 0.4 mol m −2 d −1 at the Z. marina site in December to 30.7 mol m −2 d −1 at the F. vesiculosus canopy (shallowest site) in June (Supporting Information Table S1).

Biomass of plants and macrofauna
Plant biomass at the mixed macrophyte site and at the Z. marina canopy was dominated by the eelgrass Z. marina, with occasional loose-lying accumulations of ephemeral, filamentous algae comprising up to 28% of total biomass. Total plant biomass ranged from 15 to 28 g C m −2 at the mixed macrophyte site (6-50% below-ground) and from 15 to 43 g C m −2 at the Z. marina canopy (8-25% below-ground). Photosynthetic biomass at the F. vesiculosus canopy was an order of magnitude higher, ranging from 195 to 396 g C m −2 , with ephemeral macroalgae comprising < 3%.
Macrofauna biomass generally was similar at the sand, deep mud, mixed macrophyte, and F. vesiculosus sites (annual range from 2 to 10 g C m −2 ). The Z. marina canopy had slightly higher macrofauna biomass within the range of 6-14 g C m −2 , while the mussel reef had the highest macrofauna biomass of the five shallow sites, with values ranging seasonally from 21 to 34 g C m −2 .

Benthic metabolism measurements
The flux data set we present consists of 2926 h of benthic fluxes. Individual data sets used to calculate daily metabolism rates ranged from 24 to 117 h in duration (average = 73 h) (Supporting Information Table S1).
The highest daily rates of GPP were measured at the F. vesiculosus site in June (2.4 g C m −2 d −1 ), coinciding with highest seabed PAR up to 1350 μmol m −2 s −1 and the longest photic period (h day = 20 h). Daily R was highest at the mussel site in August (1.8 g C m −2 d −1 ), coinciding with the warmest water temperatures for the year at this site of 15 C. Daily NEM ranged from strongly net heterotrophic at the mussel  Information Table S1).
In early winter (December), benthic metabolism rates were low at the five photic habitats, with daily GPP ranging from 0.1 to 0.2 g C m −2 d −1 , and daily NEM ranging from net heterotrophic to metabolic balance (−0.2 to 0.0 g C m −2 d −1 ). However, in late winter (March), under colder water temperatures and higher light availability, we measured substantially higher rates of daily GPP at all five photic habitats (0.2--1.1 g C m −2 d −1 ). Altogether, winter GPP (1 st December to 31 st March) accounted for a significant proportion of annual GPP: 28% at the Z. marina site, 22% at F. vesiculosus canopy, 30% at mussel bed, 29% at the bare sediments site, and 26% at the mixed macrophytes site. Winter was similarly important for habitat R, accounting for between 27% and 31% of annual R at the shallow sites (Supporting Information Section S3).
The regression analysis that was performed on the annual integrated rates of GPP, R, and NEM indicated that in all but two cases, the slope differed significantly from zero at the 95% confidence level (R 2 value range from 0.45 to 1.0; global mean R 2 = 0.94). Annual NEM at the Z. marina site and at the sand site was not significantly different from zero (R 2 = −0.10 and −0.16).
Annual GPP was highest at the F. vesiculosus site (0.44 kg C m −2 yr −1 ). Surprisingly, the mussel reef showed remarkably high rates of GPP (0.17 kg C m −2 yr −1 ) that were comparable in magnitude to those that were measured at the Z. marina canopy (0.16 kg C m −2 yr −1 ; Fig. 3). The mixed macrophyte canopy and the sand site had a similar annual GPP (0.10 and 0.07 kg C m −2 yr −1 ). Annual R was lowest at the aphotic muddy site (0.04 kg C m −2 yr −1 ) and highest at the mussel site (0.35 kg C m −2 yr −1 ). Despite being located at greater depth, the aphotic site had surprisingly high annual R rates, comparable to those measured at some shallower sites (Fig. 3, Supporting Information Table S1). The F. vesiculosus canopy was strongly net autotrophic on an annual basis (NEM = 0.29 kg C m −2 yr −1 ), whereas the mussel reef, despite a high annual GPP, was strongly net heterotrophic (NEM = −0.18 kg C m −2 yr −1 ). The mixed macrophyte canopy had a low annual NEM of 0.01 kg C m −2 yr −1 , while annual NEM at the Z. marina site and the sand site were not significantly different from zero. The aphotic site had an annual NEM of −0.04 kg C m −2 yr −1 (Fig. 3).

Seasonal benthic metabolism
Obtaining daily estimates of coastal productivity across multiple habitats that span highly contrasting biogenic communities has been challenging due to a lack of a common methodology allowing for comparable metabolism measurements. For this reason, soft-sediment habitats have been the focus of numerous studies, whereas measurements in rocky substrates have not been performed with the same intensity, despite these habitats being largely characteristic of coastal areas in temperate and high-latitude regions. Importantly, rocky habitats can harbor a massive biomass of autotrophs such as large canopies of macroalgae, or heterotrophs such as dense bivalve reefs. The AEC approach makes it possible to compare these contrasting habitats.
In our study, the hard-bottom substrates represented the two "extremes" in benthic metabolism. The F. vesiculosus canopy was the most autotrophic habitat we investigated, and the mussel bed the most heterotrophic. The metabolism results from the hard-bottom substrates are perhaps also the most intriguing. Whereas high annual GPP at the F. vesiculosus site can be expected due to high autotrophic biomass, with values comparing well with existing GPP estimates for macroalgal canopies worldwide (Krause-Jensen and Duarte 2016), the high annual GPP at the mussel bed was surprising, given the greater depth of this habitat (lower light availability) and the small standing stock of autotrophic biomass, which consisted of biofilms and short tufts of ephemeral filamentous algae (Fig. 1d). The high annual GPP at the mussel bed, which amounted to almost half of that at the F. vesiculosus site, and was comparable to that measured at the Z. marina canopy, suggests efficient recycling of nutrients and C between heterotrophs and autotrophs in the absence of other substantial sources of regenerated nutrients such as sediment deposits. Filter-feeding bivalves such as mussels may play a similar important role in nutrient regeneration through active filtration, deposition of feces and pseudo-feces, and efficient exchange of nutrient-rich waters near the seabed (Kautsky and Evans 1987;Nishizaki and Ackerman 2017). These processes may be important in maintaining high benthic GPP, thereby providing a significant autochthonous organic matter subsidy that was equivalent to~50% of the annual R rate of the mussel reef habitat in our case. Similar direct positive effects between suspension-feeding bivalves and benthic primary producers have been observed in tide-pool communities (Pfister 2007), and remarkably high primary production rates have also been reported for intertidal oyster beds (Volaric et al. 2018).
Despite being net heterotrophic, shallow-water bivalve reefs may synthesize large amounts of organic matter that may significantly offset their annual C demand. Due to the underlying rocky substrate, it is expected that the synthesized autotrophic biomass is either consumed by the mussels, and thus recycled within the habitat itself, or it is exported to surrounding depositional environments along with detrital pools originating from other shallow habitats such as F. vesiculosus beds . At the deeper aphotic site, the annual benthic C demand of the depositional basin sediments (0.04 kg C m −2 yr −1 ) was comparable to the annual water column NEM (0.05 kg C m −2 yr −1 for the 10 m depth photic zone) (Kuparinen et al. 1984). When considering that this site has an estimated annual C burial rate of~50% (data not shown), these results suggest significant connectivity with external sources of organic C likely originating from shallower depths (e.g., export from shallow photic habitats) or from inland waters.
The onset of winter in shallow, high-temperate benthic habitats is often characterized by overwintering of fauna in deeper  Table S1.
waters and burrows, and suppressed biotic interactions and energy flows (Möller et al. 1985). Our measurements of benthic metabolism in early and late winter show contrasting results. Whereas low rates of metabolism (GPP ≤ 0.2 g C m −2 d −1 ) are evident in December, with habitats largely being net heterotrophic (up to −0.2 g C m −2 d −1 ), late winter (March) GPP is mostly above 0.2 g C m −2 d −1 and as high as 1.1 g C m −2 d −1 at the F. vesiculosus site. Some habitats were strongly net autotrophic during this period (NEM up to 0.6 g C m −2 d −1 ). Despite cold water temperature of 1-2 C, the benthic habitats were able to photosynthesize efficiently. It is likely that sufficient light and ample nutrients in late winter following a largely heterotrophic autumn and early winter period make conditions favorable for photosynthetic production, resulting in large amounts of C being produced during winter, up to a third of annual benthic GPP. Similar observations have been made for other high-temperate and Arctic settings (Attard et al. 2014), although measurements of benthic metabolism in winter still remain scarce.

Spatial upscaling
The importance of individual habitats to coastal GPP, R, and NEM was estimated by upscaling our annual measurements to habitat distribution models (20 m grid resolution) that are available for our study area (Virtanen et al. 2018) (Fig. 4). To do this, we defined a 93 km 2 area encompassing all of our study site locations, and limited our analyses to shallow photic depths (0-5 m; 13 km 2 ), where most of the benthic photosynthetic biomass is located (Virtanen et al. 2018). Habitat distribution models for these shallow depths indicate a dominance of F. vesiculosus (5.1 km 2 ≈ 40% habitat coverage), followed by bare seabed (4.4 km 2 ≈ 34% coverage), M. trossulus (1.8 km 2 ≈ 14% coverage), and Z. marina (1.6 km 2 ≈ 13% coverage). Pelagic estimates are from literature values that are available for the study area (GPP = 0.17 kg C m −2 yr −1 , R = 0.12 kg C m −2 yr −1 , NEM = 0.05 kg C m −2 yr −1 ) (Kuparinen et al. 1984;Lignell 1990;Raateoja et al. 2004), and have been scaled linearly to water depth. Using this approach, we identify a remarkably large contribution of benthic habitats for metabolism in shallow waters: the seafloor was responsible for 93% of annual ecosystem GPP, 92% of ecosystem R, and 95% of the ecosystem's surplus autochthonous C production (NEM), which was dominated by the highly autotrophic F. vesiculosus canopies (Fig. 4)  . Furthermore, F. vesiculosus mediates around two-thirds of the total ecosystem GPP in shallow waters, and M. trossulus reefs appear to be of similar importance to Z. marina canopies for GPP in our study area. Benthic habitats are therefore hotspots of production in shallow waters. In the Baltic Sea, broad-scale estimates of benthic GPP are largely based upon microalgal production, with rates of~0.1-0.3 g C m −2 d −1 in summer (Ask et al. 2016). From our measurements, it is evident that shallow benthic habitats, and complex habitats in particular (Z. marina, M. trossulus, F. vesiculosus), can mediate areal rates of GPP that are an order of magnitude higher than these current estimates. Incorporating the contributions of complex benthic habitats into broad-scale estimates of coastal productivity is therefore important. Fig. 4. Spatial upscaling of annual metabolism rates to a 93 km 2 habitat map that is available for Z. marina, F. vesiculosus, and M. trossulus in the study area, considering only shallow photic depths (0-5 m; 13 km 2 ).The solid line on the map in (a) delineates the area considered, and the black dot marks the location of the Tvärminne Zoological Station. Numbers indicate location of the six study sites. Habitat coverage (b) was combined with areal metabolism rates to estimate gross primary production (GPP, c), respiration (R, d) and net ecosystem metabolism (NEM, e). Pelagic estimates are from literature values (see "Discussion" section) and have been scaled linearly to water depth. GPP and R values are positive whereas NEM can be positive or negative.

Biomass turnover rates
Partitioning of biomass production between different phototrophic components within each habitat (e.g., macrophytes, epiphytic microalgae) has important implications for food-web ecology and biogeochemical cycling (Sand-Jensen and Borum 1991;McGlathery et al. 2007). AEC fluxes quantify habitat-scale metabolism that include contributions from all phototrophic components. By considering the ratio between daily benthic GPP (g C m −2 d −1 ) and the estimated standing phototrophic biomass (B, g C m −2 ), we can infer the phototrophic C turnover rates (d −1 ) for different phototrophic habitats and seasons (Fig. 5). Using this analysis, we deduce that habitats with low phototrophic biomass consisting primarily of microalgae and thin ephemeral macroalgae (bare sediments, mussel bed) have rapid phototrophic biomass turnover rates (~0.07 d −1 ), whereas large perennial canopies such as F. vesiculosus turn over their biomass on much longer timescales (~0.003 d −1 ). These results are consistent with the notion that newly produced biomass from microalgae and ephemeral macroalgae is rapidly grazed and decomposed (Sand-Jensen and Borum 1991;Duarte and Cebrián 1996), whereas large perennial macrophytes grow slower (e.g., 0.03 d −1 for F. vesiculosus during the main growth period) (Sand-Jensen and Borum 1991), but they attain much larger standing biomass, providing habitat structure and contributing to storage of atmospheric CO 2 (Fourqurean et al. 2012;Krause-Jensen and Duarte 2016). Over the past century, perennial macrophyte canopies have deteriorated in the Baltic Sea due to eutrophication, favoring instead the proliferation of vast underwater fields of ephemeral macroalgae (Kautsky et al. 1986;Boström et al. 2003). This shift in the structural habitat biodiversity promotes low standing phototrophic biomass, rapid C and nutrient turnover, and potential for a larger local secondary production, at the cost of reducing overall canopy habitat structure and the C and nutrient sequestration capacity of the coastal zone (McGlathery et al. 2007).