Variable contributions of seafloor communities to ecosystem metabolism across a gradient of habitat-forming species

The contributions of habitat-forming species to the biodiversity and ecosystem processes of marine and terrestrial ecosystems are widely recognized. Aquatic plants are considered foundation species in shallow ecosystems, as they maintain biodiversity and sustain many ecosystem functions such as primary production and respiration. Despite the increasing amount of biodiversity-ecosystem functioning experiments in seagrass habitats, the effects of benthic variability on ecosystem functioning are rarely investigated across spatially variable aquatic plant habitats. Here, we quantitatively link seasonal variability in seafloor metabolism (i.e. gross primary production and community respiration) with major benthic community components (i.e. microphytobenthos, aquatic plants and macrofauna) across a structural complexity gradient of habitat-forming species (in terms of shoot density and biomass), ranging from bare sand, to a sparse mixture of plants to a dense monospecific seagrass meadow. The increasing complexity gradient enhanced the magnitude of the relationships between benthic community and seafloor metabolism. The daily average seafloor metabolism per season at the bare site was similar to the sparse site, highlighting the role of microphytobenthos for seafloor metabolism in shallow unvegetated sediments. The contribution of the associated macrofauna to the seafloor respiration was similar to the aquatic plant community contribution. Infauna was the main macrofaunal component significantly explaining the seasonal variability of seafloor respiration. However, benthic community-metabolism relationships were stronger within the plant community than within the macrofauna community (i.e. steepest slopes and lowest p -values). Understanding these relationships are a priority since climate change and biodiversity loss are reducing habitat complexity around the world, jeopardizing valuable ecosystem functions and services.


Introduction
The role of foundation species in organizing the community structure, facilitating biodiversity and adjusting ecological processes has long been recognized in a variety of marine and terrestrial ecosystems (e.g. Ellison et al., 2005;Angelini et al., 2011;Lamy et al., 2020). In marine sedimentary habitats, diverse and active faunal communities significantly contribute to seafloor functioning through complex biogeochemical interactions (Lohrer et al., 2004;Rodil et al., 2011;Villnäs et al., 2013). Seafloor heterogeneity is significantly boosted by emergent biological structures that add habitat complexity and influence biodiversity-function relationships. In coastal areas worldwide, seagrasses are typical habitat-forming species as their presence transforms the seafloor into structurally complex, highly productive and diverse habitats that maintain ecosystem functions and services such as fisheries, biodiversity, primary production and carbon sequestration (Heck et al., 2003;Hovel et al., 2002;Duarte et al., 2010;Fourqueran et al., 2012). Globally, these important ecosystems have declined due to anthropogenic activities (e.g. Waycott et al., 2009;UNEP, 2020). Recent conservation and restoration success (de los Santos, 2019) have stopped the decline in some locations, but to ensure continued and successful conservation of these coastal ecosystems we need to improve our understanding of the ecological factors driving seagrass productivity.
There is a wealth of studies on the role of seagrass biodiversity (e.g. species richness, functional diversity, shoot density or canopy height) on coastal ecosystem functioning (including carbon uptake, nutrient cycling and O 2 production) (e.g. Touchette and Burkholder, 2000;Fourqueran et al., 2012;Rheuban et al., 2014;Gustafsson and Norkko, 2019). The role of benthic macrofauna to the sedimentary oxygen dynamics has also been widely recognized (Glud, 2008), and their contribution to the seafloor metabolism in seagrass meadows has been documented. For instance, some bivalve species are expected to influence seagrass productivity due to sedimentary nutrient enrichment (Peterson and Heck, 2001;Lohrer et al., 2016). While habitat-forming species can be dominant across large spatial scales (Ellison et al., 2005;Boström et al., 2011;Angelini et al., 2011), experiments are typically performed on smaller spatial scales. Most of the experimental studies linking the activity of the benthic communities and the variability in seafloor metabolism have been based on enclosure techniques, either using benthic cores in the laboratory (e.g. Wijsman et al., 1999;Glud, 2008;Delgard et al., 2016) or benthic chambers in the field (e.g. Lohrer et al., 2004;. Hence, there is a need for large-scale studies relating the contribution of benthic biodiversity to the ecosystem metabolism of contrasting coastal habitats. Typically, seascape patterning of coastal meadows shows a mosaic of aquatic plant patches with different complexity (shoot density and biomass), where the interaction of the different elements influences the ecological functioning and biodiversity of the whole mosaic (Boström et al., 2011). Previous studies have often only focused on certain seagrass meadow types (i.e. homogeneous seagrass meadows) (e.g. Rheuban et al., 2014;Drylie et al., 2018), and hence biodiversity-metabolism relationships across the spatial heterogeneity of seagrass beds remain unknown. Bare sand interspersed with patchy and continuous mono-and multispecific seagrass are typical along the Baltic coast (Gustafsson and Norkko, 2019), making this area well suited to study how different benthic components relates to seafloor metabolism across spatially variable coastal vegetation.
New technological developments such as the aquatic eddy covariance (AEC) allows for a detailed habitat-scale measurement of primary production rates of seagrass meadows over time, highlighting the role of vegetated benthic communities as hotspots for metabolic activity in various seagrass ecosystems (Rheuban et al., 2014;Attard et al., 2019;Berger et al., 2020). Recently, we have been able to estimate the seafloor area that contributes the most to the benthic oxygen flux dynamics measured by the AEC (i.e. footprint), and simultaneously characterize the main benthic biodiversity components within the footprint that can influence the seafloor metabolism . We showed that structurally complex coastal habitats of the Baltic Sea drive most of the seasonal ecosystem metabolism, highlighting their central role in primary productivity . Based on seasonal biomass measurements, we estimated that macrofauna communities associated with aquatic plant meadows contributed up to 25% of the overall seafloor respiration, and simultaneously around 10% of the seafloor metabolism translated into secondary production (Rodil et al., 2020). However, we have not been able to quantify the linkages between the temporal variability of the seafloor metabolism with the complexity of contrasting aquatic meadows and the variety of associated benthic organisms. Here, by exploring the variation in the seafloor structural complexity within the AEC footprint area at different time scales we can formally test how benthic biodiversity components influence seafloor metabolism.
Building from previous results, we quantify the links between the seafloor metabolism (i.e. gross primary production and respiration) and the major benthic community components (i.e. microphytobenthos, aquatic plants, macroinfauna and epifauna) across a spatial gradient in meadow complexity based on increasing shoot density and plant biomass (i.e. from bare sand to sparse mixed vegetation to dense monospecific seagrass). The study compared three different seafloor habitat types (i.e. a bare sand site, a sparse vegetated site, and a dense seagrass site) located in different sites representing a natural spatial gradient in increasing complexity. We used existing data from those sites and included new sampling dates to increase the resolution on the temporal variation of the habitat characteristics to understand how seasonal changes in habitat-forming species could modulate biodiversity-ecosystem functioning relationships based on the spatiotemporal variability of benthic metabolism. Elements of seagrass structural complexity such as shoot density, biomass or canopy height have been invoked as important determinants of faunal abundance and seagrass primary production (e.g. Hovel et al., 2002;Rheuban et al., 2014;Lohrer et al., 2016;Gustafsson and Norkko, 2019). We hypothesise that the relationships between seafloor metabolism and benthic community will be stronger following an increase in structural complexity of habitat-forming species.

Study sites
The study sites were located along the Hanko peninsula, Baltic Sea archipelago, SW Finland (59.844 • N, 23.249 • E). We selected three spatially close shallow (3-3.5 m) sedimentary sites ( Fig. S1): (1) a bare sand site with no emergent structures, (2) a sparse vegetated site (composed of a diverse assemblage of limnic and marine plants with a variable and more disperse plant coverage), and (3) a dense seagrass site (mainly Zostera marina Linnaeus, 1756) characterized by different composition and structural complexity (i.e. shoot density and biomass) of the plant stands (Table S1) (see Rodil et al., 2020).
Sampling was conducted by SCUBA divers on five (bare), eight (sparse), and six (dense) occasions from June 2016 to September 2017 ( Fig. 1; Table S1). The present study is based on the combination of two published seasonal data sets and additional unpublished sampling dates to increase the temporal resolution of the habitat characteristics (Tables S1-S3). One dataset (Tables S1-S2) consists of benthic biodiversity information (Rodil et al., 2020), and the other dataset (Tables S3) consists of seafloor metabolism data (i.e. O 2 fluxes) obtained using large-scale, non-invasive aquatic eddy covariance (AEC) measurements .

Benthic oxygen flux measurements
Benthic O 2 fluxes were quantified in situ using the aquatic eddy covariance (AEC) O 2 flux method (Berg et al., 2003), an emerging technological development that allows investigating large-scale (10s of m 2 ) metabolic rates and their drivers non-invasively at a high temporal resolution . The main hardware consisted of an Acoustic Doppler Velocimeter (ADV) and fast-response oxygen microsensors all mounted in the centre of a sturdy lightweight tripod frame with its measurement volume located ~15 cm above the seabed surface ( Fig. 1). A photosynthetic active radiation (PAR) sensor (LI-192, Li-Cor) and a dissolved O 2 optode (U26-001, HOBO) were located on the frame to log at 5 min intervals PAR, dissolved O 2 concentration and temperature. Benthic O 2 fluxes (mmol m − 2 h − 1 ), were extracted from the measured velocity and O 2 microsensor data streams, and benthic metabolism rates were computed as daily gross primary production (GPP) and community respiration rates (CR) following well-established protocols (Berg et al., 2003;Attard et al., 2019) (Table S3). Net ecosystem metabolism (NEM) was computed as the difference between GPP and CR. The AEC was deployed by divers from a small boat onto a central area within each site to cover a major representation of the key biodiversity structures. Individual deployments logged flow velocity and O 2 microsensor output in continuous sampling mode at 32 Hz over 3-4 days. The flux data set we present consists of 1435 h of benthic fluxes. Individual data sets used to calculate daily metabolism rates ranged from 49 to 99 h in duration (average = 75.5 h) (Table S3). Computing daily metabolism rates from the 15min fluxes followed a three-step process (see Attard et al., 2019 for further details). Briefly, 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.

Benthic biodiversity sampling
At the end of each deployment, we conducted a standardized sampling design (see Rodil et al., 2019 for the methodological approach) to quantify dominant features of benthic biodiversity within the AEC footprint area (i.e. the seafloor area that contributes most of the O 2 flux registered). The annual metabolism and AEC footprint of the three study sites has been previously investigated . Hence, we know that the maximum contribution to the seafloor O 2 flux dynamics at these sites can be found for a seafloor circular area of approximately 80 m 2 and within a 5 m upstream distance of the AEC instrument . At each site, we took eight random samples within the AEC area to collect representative components of the benthic community. We harvested all the vegetation and associated epifauna by gently enclosing all plants within a quadrat frame (25 cm × 25 cm) into net-bags. Plant shoot density was measured by counting all shoots from the net-bags (n = 8) and all plant species were identified. Vegetation was divided into aboveground (AG: shoots) and belowground (BG: rhizomes, roots) biomass, and dry weight biomass was estimated in the lab (60 • C, 48 h). All the associated epifauna was sorted, identified and counted. Eight sediment cores (Ø = 5 cm, 15 cm deep) for infauna were randomly taken within the circular area at each site. Sediment was sieved (0.5 mm), animals were sorted, identified and counted (individuals m − 2 ). Finally, sediment surface samples (1 cm) were randomly taken within the sampling area using syringes (Ø = 3.5 cm, n = 3) for organic matter (OM, %) and microphytobenthic pigments (i.e. chlorophyll a, μg g − 1 ) analyses.
We used shoot density and biomass per sampling date as attributes of structural complexity of the plant stands (e.g. Hovel et al., 2002;Gullström et al., 2008) across sites. However, temporal changes in environmental variables (e.g. light and nutrients) are expected to result in changes to the population structure and functioning (i.e. plant physiology and morphology) of plant stands (Weller, 1987;Lonsdale, 1990). Differences in the biomass-density relationships may reflect differing ecology among meadows and, therefore, may be used as a reliable ecological indicator of plant population dynamics between meadows (Weller, 1987;Cabaço et al., 2007). We applied a typical allometric relationship between density (D) and mean plant mass (B) given by log 10 B = log 10 k -0.33 x log 10 D to scale plant growth to size (see Vieira et al., 2018). Here, B refers to aboveground biomass and k is an allometric constant. In addition, we estimated the standing stock (g C m − 2 ) of the main phototrophic (i.e. plant species and microphytobenthos) and macrofauna components at each site using conversion ratios for dry mass (macrophytes), chlorophyll a (microphytobenthos) or ash-free dry mass (macrofauna). Seagrass C content (g C m − 2 ) was estimated by applying species-specific conversion ratios (range from 35 to 41% of dry weight) (Gustafsson and Norkko, 2019). Microphytobenthos organic C (g C m − 2 ) was coarsely estimated from chlorophyll a assuming a multiplication factor of 46 (Dejonge, 1980). Macrofauna AFDM was converted to C content (g C m − 2 ) assuming an organic content of 50% (Wijsman et al., 1999).

Statistical analyses
We examined the distribution of the data using the Draftsman plot routine to assess the degree of skewness, presence of outliers and correlations. Environmental data (i.e. organic matter, pigments, PAR and temperature) was log (x+1)-transformed and normalized, and a principal component analyses (PCA) was used to summarize seasonal patterns that can affect the structure and distribution of the benthic communities across sites. The three sites were seasonally associated with temperature and PAR (Fig. S2). The lowest temperatures (<4 • C) were found from December to March, and the lowest PAR conditions (<4 mol m − 2 d − 1 ) were found from October to March (Table S3; Fig. S3). Hence, temperature and PAR were used as proxies for seasonal variability in the regression-based analysis relating seafloor O 2 fluxes and benthic community. In addition, we examined the seasonal dissimilarity of the macrobenthic assemblages (i.e. species-specific abundance) across the three sites conducting two principal coordinate analysis (PCO) for centroids (plants and macrofauna, separately). The PCOs were based on the Bray-Curtis resemblance measure calculated from 4th-root transformed variables, followed by one-factor non-parametric multivariate analyses of variance (PERMANOVA). We used site as a fixed factor, and we tested any potential seasonal interaction using PAR as a numeric covariable (log (x+1)-transformed). We performed a permutation of residuals I.F. Rodil et al. under a reduced model, and a Type I SS as recommended for analysis involving covariates (Anderson et al., 2008). Significant effects were further investigated through a series of pairwise comparisons. A permutational analysis of multivariate dispersions (PERMDISP) was applied to quantify the degree of sample-point deviation from the respective centroids defined by the seasonal variability across sites. The macrofauna taxa that mostly contributed to the similarity among sites were identified using SIMPER analysis. Daily variability of the O 2 fluxes and the main benthic community characteristics (Tables S1-S3) were analyzed through PERMANOVA (same as above). We calculated distance resemblance matrices using Euclidean dissimilarity measures based on log (x+1)-transformed data, which is equivalent to the univariate ANOVA, but where the p-value for the test is obtained using permutations (4999). Analyses were performed using PRIMER7 (Anderson et al., 2008;Clarke and Gorley, 2015).
We aim to establish links between seasonal AEC seafloor metabolic metrics as response variables (i.e. GPP and CR) and the main benthic community variables (i.e. AGDM, phototrophic biomass, k, macroinfauna and epifauna abundance and biomass) as predictor variables as well as temperature as a covariate. We used separate linear regressions for each predictor variable to estimate the average contribution of the benthic community to the daily-averaged metabolic metrics per season and included temperature or PAR as proxies of seasonal variability. Temperature and PAR were not used simultaneously in the regressions because they were highly correlated (Pearson's R = 0.65) and showed a similar pattern (e.g. lowest temperature and light in winter; Fig. S3). The normality (Shapiro test) and homogeneity (Durbin Watson test) of the residuals were estimated. Regression analyses were fitted in R 3.6.1. (R Development Core Team, 2020).

Characterizing the complexity gradient of habitat-forming species
The PCO of the aquatic plants showed a clear discrimination of the assemblages (81.8% of the variation explained by the two first axis) between vegetated sites (Fig. 1) based on dissimilarities among speciesspecific plant abundances ( Fig. 2A). The PERMDISP analysis detected significantly different multivariate dispersion caused by the low plant assemblage dispersion (high similarity) within the dense meadow compared to the sparse site ( Fig. 2A; Table S4). Thus, the aquatic plant assemblages were significantly different (47.5% dissimilarity) between the vegetated sites (Table S4). The sparse site was composed of a diverse assemblage of aquatic plants (i.e. Ceratophylum demersum, Chara spp., Myriophyllum spicatum, Potamogeton perfoliatus, Ruppia sp., Stuckenia pectinata and Zostera marina), while the dense site was mainly represented by Zostera marina (94 ± 8% of the total macrophyte dry weight).
A second PCO (Fig. 2B) showed a clear segregation of sites according to the macrofauna assemblages (61.6% of the variation). PERMDISP detected significant differences in the dispersion of the macrofauna assemblages caused by the low dispersion within the dense seagrass ( Fig. 2B; Table S4). There were significant seasonal differences in the macrofauna assemblages across sites (i.e. Site x PAR; Table S4). However, the similarity of the macrofauna assemblages within sites was higher than between sites, and there was a significant decrease in the similarity of the assemblages across the gradient ( Fig. 2B; Table S4). The bivalve Macoma balthica was the most common and abundant species across sites and over time (Table S5). The bare sand was the less diverse site, with the invasive polychaete Marenzelleria spp. (especially in winter) and two gastropod species (Peringia ulvae and Potomopyrgus antipodarum) as the most abundant species (Table S5). The sparse site showed a high number of species compared to the bare site, including bivalves (e.g. Cerastoderma glaucum and Mytilus trossulus) and polychaetes (e.g. Hediste diversicolor) as well as a number of epifauna species (Table S5). The dense site showed the highest macrofauna diversity, with a high number of typical epifauna species (e.g. Gammarus spp., Idotea balthica, I. chelipes), as well as mollusc (bivalves and gastropods) and polychaetae species (Table S5).
There were significant differences in the abundance and biomass of the macrobenthic communities following the complexity gradient of habitat-forming species (i.e. dense > sparse > bare; Figs. 3 and 4), and this pattern was consistent over time (i.e. no seasonal × PAR interaction; Table S6). Hence, we included the benthic community characteristics (e. g. phototrophic biomass or infauna biomass) as numeric explanatory variables of the complexity gradient across sites (see supplementary files for a detailed statistical description). Overall, the attributes of complexity of the plant community and related biodiversity metrics (i.e. shoot density, AGDM, k, phototrophic biomass) showed a significant increase across sites ( Fig. 3; Table S6). However, the interaction with PAR was significantly different between sites (Table S6). Thus, the relationship between shoot density and PAR was significantly negative (Shoot = 942.4-74.7 x PAR, R 2 = 0.503; p < 0.001) for the dense meadow (Fig. 3). However, the phototrophic biomass in the dense meadow showed a significant and positive relationship with PAR (Photo = 19.31 + 1.96 x PAR, R 2 = 0.163; p < 0.01) (Fig. 3). Overall, macrofauna abundance and biomass showed a significant increase following and biomass-density allometric relationship (k) between vegetated sites across sampling dates (E) Relationship between the phototrophic biomass (i.e., macrophytes and microphytobenthos; mg C m − 2 ). and PAR across sites (see Table S6).

Ecosystem metabolism across a gradient of structural complexity
Seasonal gross primary production (GPP) and community respiration (CR) computed from daily-averaged metabolism rates (mmol O 2 m − 2 d − 1 ) showed high temporal variation ( Fig. 1; Table S3). Over the year, daily net ecosystem metabolism (NEM, computed as the difference between GPP and CR) of the bare sand was close to zero, while the NEM of the vegetated sites was widely variable ( Fig. 1; Table S3). The NEM of the sparse site was mainly heterotrophic (NEM < 0) over the year (except in June 2016). The NEM of the dense site showed both net heterotrophy (2016) and autotrophy (2017) for the same seasons in different years ( Fig. 1; Table S3). Daily CR was positively related to GPP (log-CR = 0.92 + 0.013GPP, F 1,47 = 38.4; p < 0.001), showing a strong coupling (R 2 adj = 0.44) on seasonal timescales across sites (Fig. 5A). The dense site showed significantly higher rates of CR and GPP compared to the sparse and bare sites, and high seafloor metabolic rates were significantly linked to high PAR values (Table S7), typically associated to summer ( Fig. S2; Table S3), across sites ( Fig. 5B and C).

Benthic community role on the seasonal variability of seafloor metabolism
There was a positive interaction between temperature and seasonally averaged shoot plant density explaining 69% (R 2 adj ) of the seasonal GPP variability (Table 1). This result means that a positive effect of shoot density on GPP coincided with high temperature (Fig. S4). In addition, plant above ground dry mass (AGDM) explained 38% of the GPP seasonal variability across the vegetated sites and sampling dates (i.e. no temperature interaction) (Table 1; Fig. 6A). The allometric biomassdensity relationship (k) and the phototrophic biomass (i.e. microphytobenthos plus macrophytes) showed a positive relationship with GPP across the increasing gradient of structural complexity ( Fig. 6B and C), explaining 28 and 20% of the seasonal variability across sites, respectively (Table 1). Similarly, the seasonal variability of the CR was explained by the positive interaction between shoot density and temperature (Table 1). This relationship suggests, as for GPP, that high shoot density (Table S1) and temperature (Table S3) significantly promoted high rates of CR across sites (Fig. S4). The AGDM ( Table 1) explained 47% of the CR seasonal variability across vegetated sites and sampling dates (Fig. 6D). The allometric k and the phototrophic biomass showed a positive relationship with CR across the gradient of structural complexity ( Fig. 6E and F), explaining 44% and 34% of the seasonal variability across sites, respectively (Table 1).
For the macrofauna community, the abundance was positively related to the CR across the increasing gradient of structural complexity ( Fig. 7A and B). Thus, infauna and epifauna abundance explained 33.3% and 19.2% of the CR seasonal variability, respectively (Table 1). However, the positive relationship between infauna abundance and CR was significantly stronger (51.6%) with increasing temperatures (Fig. 7A). There was a significant increase in the CR with increasing infauna biomass across sites (Fig. 7C). There was a significant and positive interaction between epifauna biomass and temperature (Table 1), meaning that epifauna biomass contributed significantly to the CR with high temperatures (Fig. 7D).

Discussion
In this study we have furthered the understanding of the role of patchiness in aquatic vegetation as a driver of seafloor metabolism.
Improved understanding of such relationships pertains to our ability to address current major environmental challenges in coastal areas, such as the curtailing of biodiversity loss and the cycling and storage of carbon (e.g. Waycott et al., 2009;Duarte et al., 2010).

Ecosystem metabolism across a gradient of structural complexity
The average seafloor O 2 flux variability from the vegetated sites was similar to the high but widely variable metabolism reported across coastal meadows worldwide (e.g. Duarte et al., 2010;Rheuban et al., 2014;Attard et al., 2019). The bare site showed less seasonal variability of the seafloor fluxes compared to the vegetated sites, with a NEM close to zero over the year (Fig. 1). As expected, the relationship between CR and GPP was tightly coupled on seasonal scales across sites (e.g. Duarte et al., 2010;Rheuban et al., 2014). Seasonal changes in light (i.e. PAR) and temperature were the main drivers of ecosystem metabolism across the study sites, with high rates of CR and GPP during high PAR conditions (i.e. summer). Temperature and light are known to drive the structure of benthic communities and seafloor metabolism from coastal habitats (e.g. Delgard et al., 2016;Attard et al., 2019). A stronger metabolic-PAR relationship occurred in the dense meadow compared to the other sites (Fig. 5). Across seasons, seafloor metabolism was on average about two and three times higher in the dense meadow than it was in the sparse and bare sites, respectively (Fig. 5).

The role of primary producers on the seasonal variability of seafloor metabolism
Our study shows that an increase in the complexity of the aquatic plant community promoted an increase in seafloor metabolism across sites (dense > sparse > bare). This is expected given that seagrass density and biomass are related to the amount of photosynthetic biomass present. Hence, the dense site showed a higher range of values and a higher data variability in the coupled relationships between plant community characteristics and seasonal seafloor metabolism compared to the sparse and bare sites (Fig. 6).
Typically, we can expect that an increase in plant shoot density will stimulate productivity (Rheuban et al., 2014;Delgard et al., 2016;Lohrer et al., 2016). However, in our study, shoot density contributed positively to the seasonal seafloor metabolism variability only when interacting with temperature (Fig. S4). Aquatic plant density is directly regulated by environmental variables linked to depth and seasonality such as light and temperature (Ralph et al., 2007). This suggests that density is regulated in a more direct and deterministic manner (i.e. responds faster) by environmental conditions than for instance biomass, Table 1 Significant regressions describing gross primary production (GPP) and community respiration (CR) as the main response variables. Explanatory variables for the aquatic plants (shoot density, AGDM: aboveground dry mass and phototrophic biomass) and macrofauna (infauna abundance and biomass, epifauna abundance and biomass) communities, and temperature (numeric co-variable as proxy of seasonality).  Fig. 6. Relationships between daily-averaged gross primary production (GPP) and the average (n = 16) (A) aboveground dry mass (B) biomass-density relationship (k) and (C) phototrophic biomass, and between daily-averaged community respiration (CR) and (D) aboveground dry mass (E) k allometric constant and (F) phototrophic biomass across sites and dates. All relationships were significant at p < 0.05 (Table 1). and consequently it is the more sensitive of the plant abundance indicators to seasonal and environmental changes (Olesen and Sand-Jansen, 1994). In our study, shoot density was a more variable parameter (in terms of mean and standard deviation) compared to plant biomass. Aboveground plant biomass is the product of shoot density and mean shoot weight; however, both parameters did not necessarily correlate well over time. For instance, the highest shoot densities in the dense meadow coincided with the lowest temperatures and the lowest seafloor metabolism, but also with very low aboveground biomass (Table S1). This biomass-density mismatch could be explained by plant decaying processes or by biomass acclimatization to lower light availability during winter (e.g. Lonsdale, 1990;Olesen and Sand-Jensen, 1994). We showed that the aboveground dry mass and the allometric biomass-density relationship contributed significantly to the seasonal variability of the seafloor metabolism across sites (Fig. 6). Plant biomass has been strongly linked to production in coastal waters (e.g. Duarte and Chiscano, 1999;Barrón et al., 2004). In addition, the allometric scaling of aquatic plants can respond to different environmental conditions (e.g. light) that can reflect differing ecology among plant communities and therefore, it is a useful tool to assess anthropogenic disturbances, plant community dynamics and seagrass productivity (Duarte, 1991;Cabaço et al., 2007). Partitioning of biomass stock of different phototrophic organisms is known to have significant implications for food-web ecology and biogeochemical cycling (McGlathery et al., 2007;Attard et al., 2019). We included sedimentary chlorophyll a as a proxy for photosynthetic microphytobenthos to compare the standing stock (g C m − 2 ) of the main phototrophic components (i.e. microphytobenthos and aquatic plants) across sites. Despite having the lowest phototropic biomass across sites, the bare site showed a comparable seasonal role in the seafloor metabolism to the vegetated sites (Fig. 6). Specifically, the bare and sparse sites showed similar average rates of GPP and CR (i.e. 10-30 mmol O 2 m − 2 d − 1 ). In shallow unvegetated habitats, microphytobenthos and light availability are known to be principal determinants of benthic primary production (Sundbäck et al., 2000;Hope et al., 2019), also during winter (Attard et al., 2014;Delgard et al., 2016). Hence, bare sediments can be net autotrophic due to the high productivity of microphytobenthos (Sundbäck et al., 2000;Rodil et al., 2011), suggesting a significant role of bare sands in maintaining productive and healthy coastal systems beyond vegetated beds. This is of particular relevance in the Baltic Sea where microphytobenthos contributes significantly to the annual coastal primary production (Ask et al., 2016), and seafloor distribution models for coastal habitats indicate a dominance of bare sediments . These results support our hypothesis that the increasing gradient in the complexity of aquatic plants is responsible for some of the significant differences in the seasonal magnitude of the coupled relationships (Fig. 6). Our sites showed no significant presence of associated epiphytes during the sampling dates (Rodil et al., 2020). Therefore, we can relate the macrovegetation presence as the major macrophyte contributing to the overall seafloor productivity.

Macrofauna role on the seasonal variability of seafloor metabolism
The three-dimensional structure of the vegetation combined with the high primary productivity of aquatic plants, enhanced by benthic microalgae, provides a habitat and a supply of organic matter for the macrofauna community (Boström and Bonsdorff, 1997;Heck et al., 2003). Our study shows that the increase in the complexity by habitat-forming organisms following a spatial gradient in meadow composition increased the abundance and biomass of the macrofauna community, consequently enhancing the potential effect of macrofauna activity on seafloor metabolism.
Macrofauna activity (e.g. burrowing or biodeposition) is known to contribute to nutrient dynamics and primary productivity of sedimentary habitats around the world (e.g. Peterson and Heck Jr., 2001;Lohrer et al., 2004;Rodil et al., 2011;Villnäs et al., 2013). Our results support our hypothesis that an increase in structural complexity of habitat-forming species will indirectly influence seafloor metabolic processes such as community respiration (CR) by increasing the Fig. 7. Relationships between daily-averaged community respiration (CR) and the average (n = 16) (A) infauna abundance interacting with temperature, (B) epifauna abundance, (C) infauna biomass and (D) epifauna biomass interacting with temperature across sites and sampling dates. All the relationships were significant at p < 0.05 (Table 1). plant-associated macrofauna community. For example, the dense site showed a higher range of values and a higher data variability in the coupled relationships between macrofauna community characteristics and CR compared to the sparse and bare sites (Fig. 7). Infauna, both in terms of abundance and biomass, was the main macrofauna community component significantly explaining the seasonal variability of seafloor CR (33 and 43%, respectively). Macrofauna community typically contributes ~50% of the CR through respiration and bioturbation (Glud, 2008). In addition, the significant role of the infauna abundance on the CR variability was emphasized with increasing temperature, especially in the dense site. The infauna stimulated CR in the bare site, showing similar average rates of CR to the average rates of the sparse site. Indirectly, the activities of different infauna species are capable of enhancing not only the productivity of macrophytes (Peterson and Heck Jr., 2001;Lohrer et al., 2016), but also benefiting the productivity of microphytobenthos (Lohrer et al., 2004;Rodil et al., 2011) and microbes (Glud, 2008;Michaud et al., 2009) in marine sediments.
Epifauna showed a rather weak role on the seafloor metabolism compared to infauna. In general, infauna contributes more to CR than epifauna, because in addition to their respiration, they also re-oxidize reduced chemical species in the sediment through bioturbation (Glud, 2008). Epifauna abundance was on average an order of magnitude lower than infauna abundance across sites and over time, therefore decreasing the potential contribution to seafloor metabolism. Moreover, epifauna abundance and biomass did not always correlate well. For example, the dense site showed the highest epifauna biomass and the lowest epifauna abundance during winter, coinciding with the highest plant density and the lowest temperature, and with the lowest CR. Epifauna showed a strong positive role in the seasonal CR variability during high temperature sampling dates; coinciding with higher epifauna biomass, and higher activity and respiration rates (Rodil et al., 2020).

Conclusions
Our study confirms the significant role that the structure of aquatic plant communities plays for the seafloor metabolism of coastal habitats and highlights the role of microphytobenthos in the seafloor metabolism of shallow bare sediments. Macrobenthic biodiversity metrics (i.e. abundance and biomass) were positively related to seafloor metabolism (i.e. GPP and CR) across an increasing gradient of structural complexity of habitat-forming species. Macroinfauna contribution (R 2 = 0.43) to the seafloor CR variability was comparable to the aquatic plant community contribution (R 2 = 0.46) in terms of the average variability explained. However, benthic community-metabolism relationships were stronger within the plant community than within the macrofauna community (i. e. steepest slopes and lowest p-values). Aquatic plant meadows rank amongst the most valuable ecosystems in the world in terms of the goods and services they provide (Duarte et al., 2013;UNEP, 2020). The reduction or loss of the structural biodiversity and phototrophic biomass in shallow vegetated habitats due to human-induced activities (e.g. eutrophication and bottom trawling) will affect the productivity, nutrient cycling and carbon sequestration capacities of these heterogeneous environments, and consequently the ecosystem functions and services provided by coastal habitats worldwide (e.g. Duarte et al., 2010;Sala et al., 2021). The detrimental effects will be most likely exacerbated due to climate change (Duarte et al., 2013). The sustainability of coastal ecosystems calls for a strong societal action and formal management commitment to account for habitat conservation and the recovery of disrupted marine populations of habitat-forming species and their associated communities.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.