Towards a sampling design for characterizing habitat-specific benthic biodiversity related to oxygen flux dynamics using Aquatic Eddy Covariance

The Aquatic Eddy Covariance (AEC) technique has emerged as an important method to quantify in situ seafloor metabolism over large areas of heterogeneous benthic communities, enabling cross-habitat comparisons of seafloor productivity. However, the lack of a corresponding sampling protocol to perform biodiversity comparisons across habitats is impeding a full assessment of marine ecosystem metabolism. Here, we study a range of coastal benthic habitats, from rocky-bed communities defined by either perennial macroalgae or blue mussel beds to soft-sediment communities comprised of either seagrass, patches of different macrophyte species or bare sand. We estimated that the maximum contribution to the AEC metabolic flux can be found for a seafloor area of approximately 80 m2 with a 5 meter upstream distance of the instrument across all the habitats. We conducted a sampling approach to characterize and quantify the dominant features of biodiversity (i.e., community biomass) within the main seafloor area of maximum metabolic contribution (i.e., gross primary production and community respiration) measured by the AEC. We documented a high biomass contribution of the macroalgal Fucus vesiculosus, the seagrass Zostera marina and the macroinvertebrate Mytilus edulis to the net ecosystem metabolism of the habitats. We also documented a significant role of the bare sediments for primary productivity compared to vegetated canopies of the soft sediments. The AEC also provided insight into dynamic short-term drivers of productivity such as PAR availability and water flow velocity for the productivity estimate. We regard this study as an important step forward, setting a framework for upcoming research focusing on linking biodiversity metrics and AEC flux measurements across habitats.


Introduction
Coastal benthic zones are diverse and productive environments capable of sustaining vital ecosystem functions and providing valuable societal services [1]. Primary production and respiration are metabolism metrics commonly used for describing ecosystem functions and community health e.g., [2][3][4][5]. The contribution of different elements of the benthic biodiversity such as richness of species, abundance, biomass or range of habitats to the ecosystem metabolism is key for the functioning of coastal ecosystems e.g., [2,3,6,7]. For instance, marine macrophytes, including macroalgae and seagrass, can be dominant sources of primary production and respiration in coastal systems e.g., [7][8][9]. Microphytobenthos, also a key coastal primary producer, has a crucial role for the overall ecosystem metabolism of unvegetated marine habitats e.g., [10,11]. Coastal ecosystems also harbour a high diversity of macrobenthic consumers that contribute to the ecosystem metabolism reflected by elevated benthic O 2 uptake e.g., [3,6,12]. However, the majority of the studies on benthic biodiversity and/or ecosystem metabolism have been conducted within the same type of habitat e.g., [8,[9][10][11][12][13], and cross-site information from different habitats with complex emergent structures (e.g., seagrass, seaweeds, and blue mussels) is largely lacking.
Coastal systems are complex environments characterized by heterogeneous benthic communities [1,3]. As benthic biodiversity changes across communities and habitats, it is essential to understand the role of the habitat-specific components in the variability of the O 2 fluxes [14]. Lately, the aquatic eddy covariance (AEC) technique has been implemented to quantify rates of primary production and respiration from continuous measurements of vertical O 2 fluxes within the benthic boundary layer [15]. The AEC enables extended temporal measurements and integrates large areas of the heterogeneous seafloor under the naturally varying characteristics of the dynamic coastal environments, overcoming the limitations of previous methods for measuring benthic metabolism [15][16][17]. This rapidly developing technique provides great promise in resolving differences in benthic metabolism across contrasting habitats. The AEC has been successfully adapted to quantify O 2 dynamics in a wide variety of seafloor environments including high latitude rocky embayment [17], coral reefs [18], seagrass meadows [19], and oyster beds [20]. These pioneering studies are transforming the research on ecosystem functioning. A few studies have established some preliminary links between habitat biodiversity and AEC fluxes e.g., [18,19,[21][22][23][24]. However, there are no biodiversity surveys specifically developed within the multidirectional context of the seafloor metabolic flux measurements obtained by the AEC technique.
We need to develop simple but efficient methodologies that can relate habitat biodiversity and heterogeneity information to the seafloor area of maximum metabolic contribution measured by the AEC from multiple directions. The size and shape of the seafloor surface area that contributes most of the flux registered (i.e., 90%) at the AEC measuring point is known as the 'flux footprint' [16]. The non-invasive nature of the AEC allows for quantitative footprint flux comparisons between different habitats under true in situ conditions [15,16]. However, there is no equivalent sampling protocol developed to comparably quantify the abundance and biomass of the biota across different benthic habitats. Therefore, our goal is to propose a sampling method based on established protocols for routine biodiversity sampling of subtidal benthic communities to reliably characterize the dominant biomass elements of different coastal habitats within the main metabolic area of influence measured by the AEC. As a first step, our study attempts to provide a framework which enables comparisons to be made across benthic habitats before further studies can start constraining spatial and temporal benthic biodiversity and flux dynamics in more detail. Here, we measured the main seafloor metabolic and biodiversity metrics across a range of benthic habitats in a temperate coastal setting of the Baltic Sea archipelago. We determined the multidirectional AEC footprint area with the maximum contribution to the flux signal, and then we superimposed a sampling design to characterize the habitat-specific biodiversity contributing to the seafloor-specific footprint. We define the benthic biodiversity footprint, a term associated to the AEC footprint, as the main benthic biodiversity area of influence that contains the largest metabolic contribution measured by the AEC. Eventually, this information will be pivotal for future research, when comparing and linking benthic biodiversity measures with AEC metrics of productivity and respiration across coastal habitats.

Study sites
The archipelago coastline presents a complex structure of islands and a mosaic of soft and rocky bottom shores. We selected five shallow (� 5 m) habitats of the Baltic Sea archipelago (i.e., one site per habitat), on the Hanko Peninsula, SW Finland (Fig 1, Table 1). Three habitats were selected as main representatives of the soft-sediment communities: (1) a vegetated site comprised of mixed macrophyte species (henceforth, MM), (2) an adjacent bare sand site (henceforth, BS), and (3) a seagrass meadow (henceforth, SG), mainly comprised by Zostera Mixed macrophyte (MM) and bare sand (BS) habitats, exposed seagrass meadow (SG), littoral bladder-wrack belt (FV), and blue mussel reef (BM). A selection of photographs showing the biodiversity sampling area at the MM and BS, the AEC deployed at the SG, and a closer look to the Fucus vesiculosus and blue mussels at the BM and FV rocky beds, respectively.

Aquatic eddy covariance system, deployment, flux measurements and AEC footprint
Benthic O 2 fluxes were quantified in situ using the AEC technique [15]. Our AEC systems ( Fig  1) are very similar to the original design [25], and have been described in detail in previous studies [26,27]. The main hardware consisted of an acoustic Doppler velocimeter (Nortek) mounted vertically in the centre of a tripod frame and 2 Clark-type O 2 microsensors connected to the velocimeter via submersible picoamplifiers (response time < 0.3 s). A photosynthetic active radiation (PAR) sensor (LI-192, Li-Cor), a dissolved O 2 optode (U26-001, HOBO), and a saltwater conductivity sensor (U24-002-C, HOBO) located on the frame logged transmitted (seabed) PAR, dissolved O 2 concentration, temperature and salinity throughout each deployment at 5 min intervals. The AEC was deployed from a small boat, and a lift-bag was used by divers to gently lower the tripod onto the seabed, approximately in a central and representative area of the habitat. Once deployed, the instrument provided an accurate measure of the actual distance between the sensor measurement height and the seabed, and Table 1. Aquatic eddy covariance (AEC) deployment information and environmental characteristics (mean±SE). This table includes information on the AEC footprint area and productivity flux estimates (Gross primary production, GPP, community respiration, R, and net ecosystem metabolism, NEM) at the five habitats. Photosynthetic active radiation (PAR) is calculated as an integration of continuous measurements over 24 h. X max is the integrated area with the maximum flux contribution (see [16]). Benthic O 2 fluxes (mmol m -2 h -1 ), were extracted from the measured velocity and O 2 microsensor data streams following well-established protocols [15,27,28]. Fluxes for consecutive 15 min intervals were computed using the open-source software SOHFEA [29]. Multiple days of flux measurements were bin-averaged to produce a single, continuous 24 h time series of fluxes for each site (see Fig 2). The fluxes were aligned with PAR and separated into daytime (PAR > 0.0) and night-time (PAR = 0.0) fluxes (μmol m -2 s -1 ). These fluxes were averaged to provide a mean ± SD Flux light and Flux dark , and the PAR time series was used to determine the number of daytime hours (h light ). Daily gross primary productivity was calculated as GPP = (|Flux light |+|Flux dark |) x h light . Daily respiration rates were calculated as R = |Flux dark | x 24. The net ecosystem metabolism (NEM), representing the daily balance between GPP and R, was computed (mmol m -2 d -1 ). Positive NEM indicates a surplus of organic C and O 2 production (net autotrophy), whereas negative NEM indicates net heterotrophic communities.

AEC info
For each habitat, the size and shape of the seafloor surface area included in the AEC flux measurements were examined from measurements of seafloor hydraulic roughness (z 0 ) and sensor measurement height (Table 1), as described by [16]. The length, width, and region of maximum contribution (X max ) within the AEC footprint were calculated for each site using the equations provided by [16]. The length and width of the flux footprint were used to estimate the surface area by assuming that the footprint was elliptical in shape (Table 1).

Flow direction effects on seafloor productivity estimates
We examined the effects of the water flow direction on the AEC estimates of seafloor productivity (i.e., GPP, R, and NEM), and the corresponding changes in the location of the flux footprint (see Fig 3A). Combined with specific diversity surveys of the footprint areas, analysis of directional O 2 flux data would provide valuable information about any spatial distribution of productivity on the seafloor, its variability within-habitat, and its drivers. To do this, the software package SOHFEA [29] was used to transform the AEC orthogonal coordinate system (instrument coordinates) to Earth coordinates (ENU) by correcting for the instrument heading. This allowed aligning the AEC coordinate system to the seafloor specific diversity measurements. SOHFEA then generated a corresponding mean flow direction and flow velocity magnitude value for each 15 min O 2 flux. Directional data were grouped into 8 direction sectors, each representing 45˚increments. Data processing and windrose plots were performed in OriginPro 8.5 (OriginLab).

Defining the benthic community area of contribution: The benthic footprint
Characteristics of the AEC footprint (i.e., size and shape) can be estimated with knowledge of site-dependent factors such as water depth, seafloor roughness length scale (z 0 ), and sensor measurement height [16]. Conservative tracer tracking simulations [16] provided an estimate for the length and width of an ellipse that contributes to 90% of the flux signal, and identified the region of maximum flux contribution (X max ) within this area. The benthic conditions across various coastal habitats around the world, including our study habitats, show a range of seafloor roughness length scales (z 0 , range from 0.0005 to 0.04 m) that will affect the average of the AEC flux measurements (Fig 4A). For typical coastal seafloor settings, the estimated length of the footprint ranges from < 10 m for habitats with rough benthic surfaces (e.g., oyster beds, [20] or vegetated canopies, [30]) to > 80 m for habitats with smoother surfaces such as bare sediments in unidirectional flows [31] (Fig 4B). The AEC flux signal is not distributed evenly within the footprint area, but most of the flux signal originates from a smaller region (X max ) typically located < 5 m upstream from the instrument ( Fig 4B). In our study, we estimated that a 5 m upstream distance will capture 50 to 80% of the total flux-contributing seafloor area, whether we consider the simple bare sand site or the more complex sites. Therefore, a circularshaped benthic footprint (Ø = 10 m), covering multiple flow directions ( Fig 3A) was superimposed on the typically elliptical-shaped AEC footprint ( Fig 3B) to characterize quantitatively the main biodiversity elements (in terms of biomass) within the AEC footprint of maximum contribution. We defined the benthic biodiversity footprint as the minimum benthic area of influence that contained the largest AEC footprint contribution.

Sampling the benthic footprint: Quantifying benthic communities across habitats
We developed a simple but realistic sampling protocol based on well-established protocols for routine biodiversity sampling of subtidal benthic communities to methodically conduct comparable biodiversity characterization of the main community drivers (e.g., macrophytes and  Table 1. The circular biodiversity footprint (80 m 2 , r = 5 m) was superimposed to the AEC footprint ensuring that the maximum flux contribution measured in all the habitats (X max = 5 m) was covered. The grey scale shows a theoretical gradient of the benthic contribution to the flux registered within the footprint area (darker grey implies higher contribution closer to the device). (B) Illustration of the biodiversity footprint area displaying 8 equal wedge-sections for community characterization. Quantitative random samples (n = 8, one replicate sample per each of the 8 wedgesections) of the main benthic organisms (v macrofauna, macrophytes and macroalgae) were collected (see Methods).
https://doi.org/10.1371/journal.pone.0211673.g003 macrofauna) across different benthic habitats. We used a sampling device composed of an octagonal metallic ring as the centre point and eight diving reels attached at each side of the ring (Fig 1). From each reel, a 5-m long transect line was extended and four equal wedge-sections (+A, +B, +C, +D) plus four mirrored sections (-A, -B, -C, -D) were used to delineate the sampling area ( Fig 3B). The +A guideline was always orientated to the North to align the AEC and biodiversity measurements. Thus, we obtained a circular sampling area of approximately 80 m 2 divided into eight equal wedge-sections ( Fig 3B). Each transect line was pre-labelled at 1, 3 and 5 m from the central ring for diving orientation and photographic characterization of the seafloor (Fig 3B). Following this configuration, photographs of 25 cm x 25 cm quadrats were taken at 1, 3 and 5 meters along the guidelines (total of 24 pictures per habitat, Fig 3B) to estimate the cover of the main structural components (i.e., seagrass, macroalgae, blue mussel, and bare sediment) in each habitat. Multiple issues arise within the framework of the fundamental question of where, how, and with what frequency to take samples. The answers depend upon the kind of information explored in the habitat and the substrate types to be sampled, although resource availability usually determines these boundaries. All the sampling was conducted by divers directly after the AEC deployments.
Hard-bottom habitats. Benthic samples were collected (n = 8, one random replicate per wedge-section) at the BM reef using a modified Kautsky sampler [35] for hard bottoms (i.e., a 20 cm x 20 cm square with three metal sides and a sampling bag attached to the fourth side) (Fig 3B). A spoon was used to scrape the communities from the rock into the bag. In the lab, the mussels were distributed evenly on a water-filled tray sectioned in eight sectors. Four sectors were randomly chosen, and mussels within the sector counted and abundance calculated (ind m -2 ). Due to the large abundance of mussels, all the individuals per sample were categorised into size classes, using mesh sizes of 1, 2, 4 and 9.5 mm [36]. Then, approximately 100 from every size class were randomly selected and the length (maximum anterior-posterior axis) was measured with a Vernier calliper to an accuracy of 0.1 mm. All the invertebrate species associated to the mussels were counted, identified, and the biomass was calculated (see section 2.6).
The number and height of all the F. vesiculosus individuals (n = 8) enclosed in randomly selected PVC-frames (50 cm x 50 cm, one replicate per wedge-section, Fig 3B) were nondestructively counted (ind m -2 ) and measured (height in cm) in situ at the FV site. Eight complete F. vesiculosus mature individuals were randomly taken (one per wedge-section) cutting off the holdfast and enclosed in 0.5 mm mesh-bags (Fig 3B). In the lab, all F. vesiculosus individuals were dried to constant weight at 60˚C (accuracy of 0.1 g). The average individual dry weight (in grams, n = 8) and the average abundance of F. vesiculosus (ind m -2 , in 50 cm x 50 cm) were used to calculate the average F. vesiculosus biomass (dwt, g m -2 ) of the mature individuals. We sorted all the invertebrates associated to the individual samples of F. vesiculosus (i.e., epifauna). The epifauna was counted (ind m -2 ), identified, and the biomass was calculated (section 2.6). Random replicates per wedge-section (n = 8) were collected for samples of the macroinvertebrates directly associated to the rock (i.e., seabed fauna) using a Kautsky sampler, as performed on the BM reef. In the lab, we counted (ind m -2 ), identified and calculated the biomass of all the seabed macrofauna (section 2.6).

Macroinvertebrate biomass calculation
From all the habitats, all shell-less invertebrates (i.e., polychaetes, crustaceans and others) were blot wet weighed (accuracy of 0.0001 g). The length (maximum anterior-posterior axis) of molluscs and gastropods was measured with a Vernier calliper to an accuracy of 0.01 mm. Finally, biometric conversion factors for invertebrates of the Baltic Sea [37,38] were used to obtain biomass information for all the species (dry weight, g m -2 ). Dry weight biomass for the mussels was estimated) as: Meat weight = shell length 2.307 x 10 −4.744 [36].

Coverage of the main benthic structural components
We applied a supervised image classification technique to map all the pictures (n = 24; Fig 1  and S1 Fig), performing a maximum likelihood classification on a set of raster bands (ArcGIS 10.1 geoprocessing tool). The benthic structural components were selected, and the area covered (%) was classified according to their colour (pixel counts) using the spectral signature previously defined in a training set (ArcMap). For instance, at the BS site the microphytobenthos layer can be distinguished from bare sediment using the processed pictures and their different colour density criteria, the former being dark to light shades of green and the latter brown or black.

AEC footprint, water flow and productivity estimates
The characteristics of the footprint, calculated for a water depth of < 5 m and a measuring height of 25-30 cm, were habitat dependent ( Table 1). The length of the footprint ranged between 8.8 m (SG) and 42.7 (BS) m, while the width ranged between 1.6 (FV) and 2.0 (SG) m ( Table 1). As a result of this distribution, the footprint area was relatively large in all the habitats, ranging between 13.6 (SG) and 57 (BS) m 2 . However, the maximum contribution to the flux (X max ) at the measuring point was found for a spot at the seafloor surface located � 2 m upstream of the instrument, ranging from 0 (SG) to 2 (BS) m ( Table 1). The AEC flux data provided variable rates of benthic productivity depending on the habitat. Thus, the BS site was a net autotrophic habitat (6.5±2.0, mmol O 2 m -2 d -1 ), while the MM site was a metabolically balanced habitat (i.e., NEM � zero; 0.2±1.6, mmol O 2 m -2 d -1 ) at the time of the study (Table 1, Fig 5A). The SG site was a net heterotrophic habitat (i.e., -2.4±1.4, mmol O 2 m -2 d -1 ). The FV site showed the largest O 2 production (NEM was almost twice as high as BS, i.e., 11.2±6.3 mmol O 2 m -2 d -1 ), whereas the BM site showed the largest net O 2 uptake compared to the other habitats (NEM was 17 times more negative than in the SG site, -42.1±3.4 mmol O 2 m -2 d -1 ) ( Table 1, Fig 5A).
Directional polar plots show the water flow direction as a function of velocity and frequency and how it relates to the habitat community composition. The BS site showed a relatively uniform cover of microphytobenthos ( Fig 6A) and a unidirectional flow for most of the dataset (Fig 6B). The mean flow velocity (Fig 6B) and the PAR were 2-fold higher in day 2 than day 3 (Fig 6C). GPP and R were also higher in day 2 than day 3, and NEM was negative in day 2 versus positive in day 3 (Fig 6C). Z. marina was the most abundant plant species in the MM site (S1 File, S1 Fig), and its coverage distribution (24 ± 18 and 69 ± 19%) was variable (Fig 6A). There were two predominant flow directions, Dir A and Dir B. PAR, GPP and R were 2-fold higher in Dir B than in Dir A (Fig 6C). NEM was positive in Dir A and negative in Dir B ( Fig  6C). The SG site showed a more homogenous coverage (~60% in all directions) of Z. marina compared to the MM site (Fig 6A) within the three main flow directions (A, B, C) (Fig 6A and  6B). There were similar rates of GPP, R, and NEM for Dir A and B under similar PAR ( Fig  6C). Dir C showed a similar PAR, but higher flow velocities (33% of the time � 8 cm s -1 ) compared to other directions (Fig 6B and 6C). GPP and R were higher in Dir C than in other directions (Fig 6C), but NEM was similar and negative in all three directions (Fig 6C). F. vesiculosus clearly dominated the FV site (Fig 7A, S1 Fig) and the unidirectional flow (Fig 7B) showed a positive NEM (Fig 7C). No comparisons between days and sectors are possible since we only have one complete day of data for this site. The BM site was homogenously covered in all the sampled sections by M. edulis, ranging from 60 ± 2 to 82 ± 10% (Fig 7A). The mean flow velocity for the direction sectors ranged from 0.5 ± 0.4 to 3.2 ± 1.5 cm s -1 (Fig 7B). There was no light response in the eddy fluxes (i.e., no measurable primary productivity), thus the O 2 fluxes represent R only. The mean R for the direction sectors ranged from 19 ± 2 (n = 55) to 67 ± 4 mmol O 2 m -2 d -1 (n = 104) (Fig 7C).

Habitat biodiversity characterization: The role of benthic consumers and producers
There were significant differences in the biodiversity characteristics of the dominant benthic consumers (see Fig 5B, S1 File). The lowest number of macroinvertebrate species was found at the BS site (6), while the rest of the sites presented a similar number of species (10-13) ( Table 2, S1 File). The abundance (Pseudo-F 4,35 = 61.1; p < 0.001) and biomass (Pseudo-F 4,35 = 4.1; p < 0.01) was significantly different across habitats ( Fig 5B, Table 2, S1 File). The highest and the lowest biomass were related to the BM and FV sites, respectively ( Fig 5B,  Table 2). The biomass of the mussels was significantly different (i.e. pairwise test: FV = MM = BS = SG < BM; p < 0.01) compared to the other habitats ( Fig 5B, Table 2).

Discussion
Coastal systems are connected by highly heterogeneous and productive benthic habitats. AEC measurements integrate over a large and heterogeneous seafloor surface area, and allow comparisons to be made between specific sites with contrasting structural biodiversity elements. However, the lack of a standardized sampling protocol to perform a realistic characterization of the dominant biodiversity features when using AEC measurements is impeding a comprehensive assessment of marine ecosystem metabolism across coastal habitats. To our knowledge, this is the first time that benthic biodiversity surveys across habitats have been performed exclusively upon the multidirectional AEC flux footprint characteristics (i.e., length, width, and area of the footprint and maximum contribution to the flux, X max ). The  Table 2. Mean (±SE) abundance, biomass and the number of species observed at the five habitats. All the macrofauna was divided between those species associated to the macrovegetation or macroalgae (i.e., epifauna), and those species directly related to the seafloor substrate (i.e., seabed fauna), either buried in the sediment or attached to the rocky bed. Abundance (individuals m -2 ), biomass (dry weight, g m -2 ). AEC technique improves the estimates of whole-system benthic metabolism incorporating different structural biodiversity components at an unprecedented level [15,16]. Whereas it is often unpractical to survey the benthic biodiversity on the scale of the entire AEC footprint, by focusing on the area closest to the instrument it is possible to obtain significant information about the biodiversity components that contribute the most to the AEC flux signal. In this study, we defined a seafloor area of approximately 80 m 2 with a 5 meter upstream distance as the largest biodiversity area of influence that contains the largest AEC footprint contribution measured (i.e., 50-80%, see section 2.4.). At the bare sand site, the selected upstream distance captured the lowest AEC footprint contribution, i.e., approximately 50% of the total flux-contributing seafloor area. However, bare sand sites are homogeneous habitats with less variable sedimentary characteristics, benthic biodiversity and environmental conditions compared to a heterogeneous landscape e.g., [42]. Therefore, we are confident that the biodiversity components captured in the bare sand site, following our sampling protocol, are representative of the main benthic biodiversity contribution to the AEC flux signal in this homogenous habitat. We have measured five major shallow habitats from the Baltic Sea coast with very specific biodiversity features. We determined the biomass of the habitat-specific benthic organisms to understand the role of the dominant benthic components in the system, and to discern the broad effects of biodiversity on the AEC measurements, specifically when sampling seafloor habitats with contrasting structural components (e.g., seagrass meadows or blue mussel reefs). However, there are other key coastal habitats supported by different and also important benthic biodiversity elements (e.g., oyster reefs or corals) whose idiosyncrasies must be considered when establishing potential links between biodiversity contributors and ecosystem metabolism. While our protocol presents support for its specific use in quickly estimating community biomass and species composition during the sampling of the main benthic habitats of the Baltic Sea, the procedure possesses inherent versatility for implementation in a wide variety of benthic habitats with other structural biodiversity elements and with different seafloor roughness lengths. The scope of this study was not to draw direct relationships between benthic biodiversity elements and seafloor O 2 eddy fluxes, but to provide a realistic sampling design to quantify dominant features of biodiversity within the seafloor area of maximum AEC metabolic contribution across a range of representative habitats.

Benthic biodiversity metrics and productivity estimates across-habitats
Habitats with different properties of structural biodiversity (i.e., macroinvertebrates, macrophytes and microphytobenthos) also functioned very differently. In general, the net ecosystem metabolism of the study habitats followed a productivity gradient from net autotrophic to net heterotrophic related to the main structural biodiversity components. The bladder-wrack (FV) and bare sand (BS) sites were net autotrophic (i.e., positive NEM), the mixed macrophyte (MM) site was at metabolic balance (i.e., NEM � zero), and both the seagrass (SG) meadow and the blue mussel (BM) reef were net heterotrophic (i.e., negative NEM). There were no significant amounts of epiphytes, ephemeral algae or detritus accumulations in any of the study sites (3-6% cover calculated from pictures and ArcMap). The ephemeral species basically responsible for high biomass in this region occur in spring an early summer [43], and the maximum occurrence of mobile detritus accumulations coincides in summer (see [44]). Therefore, the presence of macrophytes (i.e., F. vesiculosus, Z. marina and other angiosperms) suggests an overall high contribution of the macrovegetation biomass to the productivity of the vegetated habitats. Seafloor communities made of canopy-forming macroalgae such as Fucus spp. are regarded as autotrophic entities [9,30,[45][46][47]. Thus, the bladder-wrack community showed the highest average GPP and the lowest average R compared to the rest of the study habitats, and the NEM rendered a net autotrophic metabolism. It has been suggested that a low bladderwrack respiration during night-time could be related to a gas storage effect within gas spaces of F. vesiculosus [30]. The presence of aquatic angiosperms, mainly Z. marina, is known to influence the productivity estimates, increasing the metabolism of the seagrass meadows [8,23,24,48]. The differences in the net ecosystem metabolism between the SG and the MM sites could be related to the different shoot densities between sites, thus seagrass shoot density in the SG site was almost three times higher than in the MM site. When shoot density is very high, although the upper parts of the leaves may be light-saturated, the understory can show high respiratory losses because of the shelf-shading effect [8] causing light attenuation. Interestingly, we documented a high productivity role of the bare sediments compared to the vegetated-canopies of the soft-sediment sites (i.e., SG and MM). The productivity estimates did not confirm the general hypothesis that canopy-forming vegetation are far more productive than soft sediments without macrovegetation [13]. Alternative primary producers represented by a uniform sedimentary coverage of microphytobenthos (approximately a 60% cover of the total average estimated through ArcMap) can contribute to the metabolism of the unvegetated BS site through the benthic microalgal role on system productivity and trophic dynamics e.g., [11,49,50]. The relatively high community R recorded in all the soft-sediment sites may also indicate an increase in the sediment heterotrophic O 2 demand through the respiration and bioturbation of the benthic macroinfauna [4,48]. The high respiration rates of the filter-feeder M. edulis, although seasonally variable due to reproduction periods and food availability, are expected to play a major role on the net heterotrophy of hard bottoms [51,52]. The BM reef shaped a dense heterotrophic community expected to exert an increase in the biological O 2 demand through respiration compared to other habitats with less abundant macroinvertebrates or with a larger presence of primary producers. We acknowledge that our data on the metabolic budget for the different habitats only represents one sampling season. Thus, we need to generate a more representative data set by means of sampling across different spatial and temporal scales if we want to formally establish meaningful comparisons of oxygen fluxes between habitats, and to build significant relationships between biodiversity metrics and productivity estimates (see section 4.3.).

Biodiversity and oxygen flux measurements within habitats
Marine benthic habitats are normally systems of intense, but variable metabolism. Thus, the regulation of metabolic processes is related to a combination of environmental drivers, including light, temperature and flow velocity [23,50,53,54]. Understanding the variability of these environmental drivers is key for establishing the role of benthic communities on the ecosystem metabolism of coastal habitats. For instance, the short-term response of the benthic metabolism to the environmental drivers in a given seagrass meadow can be altered by the water flow conditions [55,56]. It is evident from our measurements that the average productivity estimates (i.e., GPP and R) in the soft-sediment sites were affected by changes in the irradiance and the water flow. For instance, the single water flow direction recorded in the BS site rendered variable productivity estimates over time, thus the highest productivity rates were related to those days with higher flow velocity and PAR conditions. During prolonged AEC deployments (i.e.,. 3-5 days) we often observed prompt changes in velocity magnitude and direction. Both MM and SG sites showed multiple flow directions with variable velocities and frequencies that affected the productivity estimate. At MM, macrophyte coverage was heterogeneous and patchy, and the highest productivity rates were recorded across the direction sector with the highest PAR conditions and the largest flow velocities (i.e., Dir B). In the SG site, light conditions were similar and the seagrass was relatively homogeneous across all the sectors, consequently the NEM was similar (~-3 mmol O 2 m -2 d -1 ) in all three main directions. However, the maximum productivity rates corresponded to the direction sector with the highest flow velocity (i.e., Dir C).
Light is considered a critical factor controlling the growth and function of benthic producers (i.e., macrophytes and microphytobenthos), thus governing the ecosystem metabolism response of coastal habitats [48,49,57]. There are a number of physical (e.g., sediment resuspension and advection, shelf-shading) and water-quality processes (e.g., nutrient loads and algal proliferation) that can affect light penetration to the seabed at different spatiotemporal scales [8,58]. Water flow has been considered a main driver of photosynthesis rates by modifying the responses of light, nutrient availability and temperature of vegetated beds [55,58]. For instance, moderately reduced flow and turbulence within vegetated beds may indirectly increase light availability by reducing self-shading on macrophyte canopies and by increasing water transparency due to sediment deposition (see [59]). However, vegetated beds can benefit from increased flow and turbulence due to faster removal of undesired substances and intensified exchange between the overlying water column and that within the meadow [55,58].
The hard-bottom sites were the most contrasting habitats in terms of the biodiversity, productivity estimates and water flow measured. The FV site is an exposed community with high annual averaged wave amplitudes and water flow conditions [30]. The unidirectional flow at the FV site was based on only one complete day of data due to sensor damage incurred~35 h into the measurements, and thus no integrated comparisons between days and flow directions are possible. The light availability in the FV site was similar to the autotrophic (BS) and metabolically balanced (MM) soft-sediment sites at the time of the study. However, it has been suggested that a high proportion of the F. vesiculosus biomass is photosynthetically-active when compared to other vegetation species such as seagrass, and thus are more efficient at photosynthesizing [60]. The large and photosynthetically-active canopy standing biomass from the FV site is present year-round (see [30]), turning this site into a net autotrophic system for most of the year, with higher annual NEM values than some of the most productive seagrass beds worldwide [30]. The O 2 fluxes recorded in the BM reef represent only respiration that can be quantified within the different direction sectors. The blue mussel cover was homogenous and abundant in all the sectors, and the multiple flow velocity and frequency was variable, and thus not directly related to the metabolic rates. However, physical episodic drivers such as upwelling occurrences or high flow velocity-induced O 2 uptakes due to e.g., advective flushing of anoxic pore-waters and resuspension might have a major role on the variability of the eddy fluxes measured in the BM reef.

Promising leads and ways forward: A future perspective
An improved understanding of seafloor metabolism requires an appropriate characterization of the community biodiversity. The relationship between benthic community biomass and seafloor metabolism measured by AEC (i.e., GPP and R) provides a compelling approach for meeting this requirement. The biomass of marine benthic organisms can be easily measured, enabling direct comparison across benthic habitats regardless of the specific structural biodiversity components. Importantly, when trying to establish links between biodiversity metrics and AEC flux measurements, it is key to have a broad ecological knowledge of the system under study, e.g., through identification of the prevailing structural features and the phenology of the dominant organism that may dictate the productivity rates of the system.
We acknowledge that the biodiversity metrics and the AEC flux measurements in this study were only taken once per habitat. The high natural variability of most of the benthic ecosystem processes, and the potential importance of scaling issues, imply that the relationship between biodiversity metrics and productivity rates in the seafloor will change at different spatiotemporal scales [1]. Thus, the coastal habitats from the Baltic Sea exhibit highly heterogeneous seafloors bound to constant changes in the environmental drivers and in the seafloor features. Regular changes in the environmental conditions such as current velocity, temperature, nutrients, or light can modify the O 2 flux dynamics of the seafloor e.g., [27,55,61]. We regard the present paper as an important step forward, setting the stage for upcoming research focusing on spatial biodiversity-productivity relationships and seasonal across-habitat productivity. Therefore, future research trying to compare and link benthic biodiversity measures with AEC metrics of productivity and respiration needs to increase the replication of the different habitats to cope with increasing spatial heterogeneity. Moreover, it is crucial to conduct seasonal sampling to build confidence in the AEC metabolic measurements and to establish links to biodiversity across different benthic habitats, spanning periods of high productivity in summer to winter conditions with low productivity [24,30]. Also, seasonal sampling is key to detect and understand concurrent changes in biodiversity metrics (e.g., biomass) and sediment-water O 2 fluxes. Long-term (2-3 months) observatory platforms (e.g., landers, sensors, video) for coastal application studies deployed to record physical, biogeochemical and/or biological activity can help to relate major environmental changes in the water (e.g., algal blooms) with the AEC measurements. Other sampling methods and non-destructive techniques such as photogrammetry, multibeam sonar, or ROV surveys for biodiversity characterization of other specific habitats could be incorporated depending on the habitat-specific requirements [62][63][64]. The natural heterogeneity of the benthic communities needs to be considered in planning AEC measurements because seafloor heterogeneity increases with scale, affecting the benthic biomass and consequently the NEM of the habitats [3,21,48]. Few empirical studies have documented how the dominant drivers of ecosystem metabolism shift across multiple scales in heterogeneous environments (see [3]). Following our sampling design, coupled relationships between biodiversity metrics (e.g., biomass) and AEC measurements (e.g., GPP) could be properly established by increasing the precision of sampling (i.e., the number of replicates), consequently improving the estimates of variation (i.e., statistical power), and the confidence of the conclusions drawn from potential biological links. Analysing the strength and variability of such coupled biological relationships following steep environmental and biodiversity gradients within habitats can provide useful ecological insight on the role of biodiversity for ecosystem functioning and benthic community health [3,14,65]. For temporal studies on the same habitat (i.e., revisiting sites), we recommend not to oversample to avoid destruction of the habitat biodiversity contributors that could potentially affect the outcome of the productivity rates. Given the novelty of the AEC measurements in the aquatic environment, there is limited experience with the approach and several aspects of the technique still need to be investigated. However, the use of the AEC approach combined with a standardized biodiversity sampling protocol opens new avenues for marine ecology research using various metrics of ecosystem functioning such as light efficiency, standing biomass, macrofauna respiration, nutrient retention or turnover rates as common currencies for potential across-habitat comparisons. Also, the implementation of the AEC technique combining cross-scale interactions provides a new perspective on the role of biodiversity and heterogeneous habitats for ecosystem metabolism, functioning and services at scales of societal relevance and appropriate for coastal management.
Supporting information S1 File. Habitat biodiversity characterization. Detailed information on the benthic biodiversity (producers and consumers) across all the coastal habitats. (PDF)