In Search of a Field-Based Relationship Between Benthic Macrofauna and Biogeochemistry in a Modern Brackish Coastal Sea

Gogina, Mayya; Lipka, Marko; Woelfel, Jana; Liu, Bo; Morys, Claudia; Böttcher, Michael E.; Zettler, Michael L.

doi:10.3389/fmars.2018.00489

ORIGINAL RESEARCH article

Front. Mar. Sci., 19 December 2018

Sec. Coastal Ocean Processes

Volume 5 - 2018 | https://doi.org/10.3389/fmars.2018.00489

This article is part of the Research Topic Living Along Gradients: Past, Present, Future View all 28 articles

In Search of a Field-Based Relationship Between Benthic Macrofauna and Biogeochemistry in a Modern Brackish Coastal Sea

$\r\nMayya Gogina*$ Mayya Gogina¹^*

Marko Lipka²

Jana Woelfel³

Bo Liu²^†

Claudia Morys⁴^†

Michael E. Böttcher²

Michael L. Zettler¹

¹Biological Oceanography, Leibniz Institute for Baltic Sea Research (IOW), Warnemünde, Germany
²Geochemistry & Isotope Biogeochemistry, Leibniz Institute for Baltic Sea Research (IOW), Warnemünde, Germany
³Chemical Oceanography, Leibniz Institute for Baltic Sea Research (IOW), Warnemünde, Germany
⁴Institute for Biosciences - Marine Biology, University of Rostock, Rostock, Germany

During several cruises in the southern Baltic Sea conducted in different seasons from 2014 to 2016, sediment cores were collected for the investigation of pore-water biogeochemistry and associated nutrient fluxes across the sediment-water interface. Six stations were positioned along a salinity gradient (ranging from 22 to 8) and covered various sedimentary habitats ranging from mud to sand. Integrated fluxes of nutrients in the supernatant water and sediment oxygen consumption were additionally derived from incubations of intact sediment cores. Subsequently, sediment from the pore-water and incubation cores was sieved for taxonomic identification and estimation of benthic macrofauna density. This combined dataset was used to determine the dominant factors influencing the vertical distribution of geochemical parameters in the pore-waters of the studied habitats and to find similarities and patterns explaining significant variations of solute fluxes across the sediment-water interface. A statistical relationship between the thickness of sulfide-free surface sediments, solute fluxes of sulfide, ammonium, and phosphate as well as oxygen consumption and taxonomic and functional characteristics of macrobenthic communities were tested. Our data and modeling results indicate that bioturbation and bioirrigation alter near-surface pore-water nutrient concentrations toward bottom water values. Besides sediment properties and microbial activity, the biogeochemical fluxes can further be explained by the functional structure of benthic macrofauna. Community bioturbation potential, species richness, and biomass of biodiffusers were the best proxies among the tested set of biotic and abiotic parameters and could explain 63% of multivariate total benthic flux variations. The effects of macrobenthos on ecosystem functioning differ between sediment types, specific locations and seasons. Both, species distribution and nutrient fluxes are temporally dynamic. Those natural patterns, as well as potential anthropogenic and natural disturbances (e.g., fishery, storm events), may cause impacts on field data in a way beyond our present capability of quantitative prediction, and require more detailed seasonal studies. The data presented here adds to our understanding of the complexity of natural ecosystem functioning under anthropogenic pressure.

Introduction

Ecosystem functioning, defined as stocks and fluxes of energy and material within and across the boundaries of a given system and its relative stability over time (Paterson et al., 2012), is mediated by diverse biological communities (Snelgrove et al., 2018). Among all the interconnected processes and ecosystem components, marine benthic macrofauna is involved in the provision of ecosystem services as sediment (de)stabilization, food webs, nutrients cycling, and transport. It has long been demonstrated that the biogeochemistry of coastal sediments is strongly influenced by the presence and activity of benthic invertebrates (Gray, 1974; Rhoads, 1974; Haese, 2006; Snelgrove et al., 2018). Though this is recognized, the actual quantitative estimates linking major parameters describing macrofaunal community inhabiting the sediment (such as its taxonomic and functional diversity and density) with nutrient fluxes within the sediment-water interface (SWI) are still sporadic and scarce. Most of the published data is based on laboratory or mesocosm experiments (e.g., Mermillod-Blondin et al., 2004; Bonaglia et al., 2013). In respect of the Baltic Sea, the results of several field studies have now been published, e.g., for the Danish Sounds and for the Archipelago Sea (Hansen and Kristensen, 1997; Karlson et al., 2007; Delefosse and Kristensen, 2012; Norkko et al., 2012). Gammal et al. (2017) emphasized the increase of solute effluxes of phosphate and ammonium from sediments with decreased macrofaunal activity induced by coastal hypoxia in the western Gulf of Finland. The recent study in the Bay of Gdansk by Thoms et al. (2018) determined the greatest impact of macrofauna on sedimentary fluxes at stations where communities were dominated by deep-burrowing polychaetes. The authors suggested that these communities can serve as a descriptive indicator to estimate the extent of the coastal filter, where organic matter is reworked and relatively large amounts of nutrients ( ${NH}_{4}^{+}$ , ${NO}_{3}^{-}$ , ${PO}_{4}^{3 -}$ , SiO₂) are released into the bottom water. Organic matter and nutrients cannot accumulate in the sediment because they are immediately used. They are buried by deep-burrowing fauna, but are also made available again to microbial degradation, therefore high turnover rates are induced due to particle mixing and surface increment (Rusch et al., 2006).

The effects of changes in benthic macrofaunal activity (both bioturbation and bioirrigation) on biogeochemical cycling are not firmly established, and few quantitative rules exist to predict their role (Mermillod-Blondin et al., 2004; Braeckman et al., 2010). Several reviews have attempted to consolidate the recent findings related to diverse ecosystems [for instance with a focus on the role of invasive ecosystem engineers by Guy-Haim et al. (2017)]. There is a certain lack of consensus even with regard to the direction of the effects. Increased faunal activity imply more intense redistribution of sediment particles, oxygenation of the sediment, decline of hydrogen sulfide, reduction of phosphate release into the bottom water (as phosphates are accumulated in the sediment bound to iron minerals), general decrease of the ammonium outflux, and an increase in the dinitrogen gas efflux due to coupled nitrification/denitrification in bioturbated sediments (Karlson et al., 2007; Norkko et al., 2012; Janas et al., 2017). Contradicting with the above, other authors (e.g., Pelegri and Blackburn, 1994; Mermillod-Blondin et al., 2004; Jordan et al., 2009; Gilbertson et al., 2012) describe the increased ammonium flux from the sediment to the overlying water column observed in the presence of burrowing invertebrates. In the Gulf of Finland, Maximov et al. (2014) found evidence for an effect of an increased population of invasive polychaete Marenzelleria spp. on binding of phosphate in sediments and reducing eutrophication conditions. This all suggests not only natural variability in nutrient cycling between regions, but also high context-dependency due to the complex interplay of multiple processes even on a very local scale (e.g., within few kilometers or less, Gammal et al., 2017), making basin scale extrapolations difficult.

There is a lack of such estimates, especially based on field data (in situ derived), in respect of the German part of the south-western Baltic Sea. This underlines the importance of further field studies and quantitative estimates, crucial for understanding and parameterization of the biogeochemical processes in order to support scientifically sound ecosystem management, to ensure the health of the Baltic Sea ecosystem and to mitigate unfavorable changes.

Overall physical transports between the free water column and sediment interstitial water, low in diffusion-dominated and high in advection-dominated habitats, control microbial activity inside sediments (Aller, 1994). According to Mermillod-Blondin (2011) bioturbation (particle mixing) and biodeposition (settling of feces and pseudofeces) produced by benthic invertebrates may play a bigger role in diffusion-dominated habitats where they can significantly modify water and particle fluxes at the water–sediment interface, whereas only slight influence of ecosystem engineers is expected in advection dominated habitats predominantly controlled by hydrological processes. Here it is important to specify that generally diffusion is the transport of dissolved material along a concentration gradient due to random movement. Advection is the transport of both dissolved material and its solvent. This includes hydrodynamics, bioturbation, bioirrigation, bottom trawling, and every process that is moving the pore water itself together with chemical species dissolved in it. However, Mermillod-Blondin (2011) referred “diffusion-dominated habitats” to marine coastal benthic environments with fine sediment texture and low physical exchange of water between bottom water and interstitial water, characterized by low interstitial flow rates. In contrast, “advection-dominated habitats” were defined as benthic environments (in shallow areas) characterized by coarser permeable sandy sediments, where high rates of movement of solutes with water are driven by hydraulic gradients. With this definition Mermillod-Blondin (2011) suggested that with high interstitial flows arising due to high rates of hydrological exchanges, bioturbation by invertebrates can only moderately modify the water circulation patterns in sediments. This study highlighted the need to test this general framework on a wide range of marine (and freshwater) habitats by quantifying the interactions between the functional traits of species and the fluxes of water and particles at the SWI.

Linking the detailed biological and chemical information is a crucial step toward a mechanistic understanding of the role of macrofauna in the marine sediments biogeochemistry and nutrient cycles. In this context, the present study aims to interpret field data obtained within a 21 month period during several cruises from sandy and muddy sediments of the southern Baltic Sea. Research was carried out within the interdisciplinary research project SECOS (“The Service of Sediments in German Coastal Seas”), aiming to determine and understand the performance of marine sediments in the German part of the Baltic Sea. The main objectives of this case study were (1) to determine the (regional and context-dependent) quantitative relationships between biogeochemical solute fluxes and macrofaunal data derived from the sediment cores, and (2) to analyze the relevance of descriptors and indicators of benthic communities functional structure that significantly influence sediment biogeochemistry. Our study is centered around the key research question: what interactions (if any) occur between the macrobenthos composition and quantitative changes in solute fluxes in muddy and sandy habitats (under modern anthropogenic impacts)? In order to address it, the following hypotheses were set off and tested based on the available data: (H1) the signals of intensive faunal bioturbation activity are reflected in pore-water profiles, (H2) there is significant correlation between biogenic mixing depth and sulfide appearance depths, (H3) the variability of biogeochemical patterns can at least partly be explained by functional structure of macrofauna.

Materials and Methods

Study Area

In the southern Baltic Sea (Figure 1), the shallow areas along the shore and on top of the offshore glacial elevations are characterized by a mosaic of coarse sediments and sands, whereas increasing water depths and finer and more organic-rich muddy sediments generally dominate in the deeper basins and bays. Fine sands cover about 32% of the German EEZ, mud up to 25% (Tauber, 2012). Near-bottom salinity and oxygen concentration are the major factors influencing species richness and composition of macrozoobenthic communities in the area (Zettler et al., 2017). Salinity declines from 20–25 in the western part of Kiel Bay to 7–8 in the Oder Bank. The water exchange between the western Baltic and the Baltic Proper is inhibited by the Darss and Drodgen Sills, usually causing the highest temporal variability of salinity in the western part of the study area. However, due to the influence of major Baltic Inflow event in December 2014 (Mohrholz et al., 2015) during SECOS cruises the variability of salinity was clearly higher in AB than in LB and MB (Table 1). Aperiodic seasonal oxygen depletion events occur in the deeper areas of the Kiel Bay, the Bay of Mecklenburg and in the Arkona Basin, and have negative effects on the diversity and density of soft-bottom fauna (Arntz, 1981). Such oxygen depletion event was observed in our study during the cruise EMB111 at stations MB, LB, ST, and AB, when only OB stations remained oxygenated (see Table 1).

FIGURE 1

Figure 1. Map of the study area in the south-western Baltic Sea with locations of sampling stations. Background colors show the distribution of sediment types after Tauber (2012). The number of cores obtained during each sampling event is given in Table 1. Geographic locations: LB, Bay of Lübeck; MB, Bay of Mecklenburg; ST, Stoltera; TW, Tromper Wiek; AB, Arkona Basin; OB, Oder Bank.

TABLE 1

Table 1. The number of cores obtained for porewater analysis and incubation experiment during each sampling event as well as values (mean ± standard deviation) of physical, chemical, and sedimentary variables measured at each station.

It is important to mention that the presence of anthropogenic effects (bottom trawling traces) was clearly observed at station MB during side-scan sonar survey (Lipka et al., 2018; see Bunke, 2018). Further, at the shallow OB station, sediments rearrangements during short-term natural physical events (e.g., storms) are known to occur (Mohrholz et al., 2015; Lipka et al., 2018).

Biogeochemical Data

Several research cruises were carried out between April 2014 and January 2016 with research vessels “Alkor” (ALK434), “Poseidon” (POS475), “Elisabeth Mann Borgese” (EMB100, EMB111), and “Maria S. Merian” (MSM50), see Table 1. During each cruise and at each visited location, sediment sampling was carried out by collecting up to eight parallel short sediment cores with a multicorer in acrylic tubes (60 cm long, 10 cm internal diameter). Cores disturbed during sampling (e.g., with shattered or inclined sediment-water interface) were discarded.

Pore-water (PW) samples were taken immediately after recovery of cores from 1, 2, 3, 4, 5, 7, 9, 11, 15, and 20 cm below and from bottom water right above the SWI using rhizons (with 0.2 μm pore size of the filter membrane) attached to clean syringes (Seeberg-Elverfeldt et al., 2005). About 10 ml pore-water were extracted from each depth (up to 20 cm below the sediment surface) after discarding the first 1–2 ml. A detailed description of the sampling procedure and pore-water analysis is given in Lipka et al. (2018). Following the methods described in Kowalski et al. (2012) and Reyes et al. (2016), the concentrations of the major, minor and trace elements (including total P, Fe, and Si) dissolved in pore-water, were measured by inductively coupled plasma optical emission spectroscopy (ICP-OES; Thermo, iCAP 6300 Duo, Thermo Fisher Scientific GmbH, Dreieich, Germany), ${NO}_{3}^{-}$ and ${NH}_{4}^{+}$ concentrations—using QuAAtro nutrient analyser (SEAL analytical) as described in Winde et al. (2014) and H₂S (sum of the dissolved sulfide species) concentrations—photometrically according to Cline (1969). Benthic solute reservoirs were calculated as depth integrated concentrations in the pore-waters of the top 10 cm surface sediments considering porosity that was determined as outlined in Lipka et al. (2018).

Incubation (INC) of intact short sediment cores (10–20 top centimeter of the sediment column) with naturally inhabiting macrofaunal communities was carried out on board of the research vessel to estimate the total oxygen uptake and total benthic fluxes of NH₄, PO₄, and Si. The volume of each “incubation chamber” (height of sediment as well as of the water phase) was determined and integrated into the calculations. Cores were incubated at constant temperature (4°C, except for during EMB111 in autumn 2015–10°C, see Table 1 for more details) in darkness for 5–10 days. Stirring mechanisms were used to homogenize the water column for analyses without resuspension of sediment particles, keeping the bottom water movement close to natural conditions. In the closed incubation chambers, overlaying water oxygen was measured using in-core attached optode spots (Fibox 4 sensor system, PreSens Precision Sensing GmbH, Regensburg, Germany): in each core, optical phase values were logged with a 5 s rate for 5–10 min three to five times a day. These values were converted into oxygen concentrations and adjusted to prevailing temperature and salinity using the mathematical conversion matrix of the manufacturer PreSens. Incubation periods under oxic conditions ranged from 3 to 10 days during EMB100, from 1.5 to 6.5 days during EMB111 when hypoxic conditions were observed at all stations with the exception of OB, and from 3 to 6 days during MSM50 cruise.

Sensor calibrations of both O₂ procedures were performed in ambient sea water (100% atmospheric saturation) aerated for 30 min, and in saturated seawater sodium dithionite solution (0% oxygen), cross checked by Winkler titration (Winkler, 1888). The deviations of oxygen concentrations measured several seconds apart were 0.2%.

For nutrient analysis, supernatant water sampling was performed every 12 to 24 h. The stirrer was switched off during the sampling to prevent intermixing of the sampled water with injected reservoir water for volume compensation. About 10 ml was filled into various sample tubes for subsequent analysis of main and trace element concentrations, sulfide/sulfate concentrations, dissolved inorganic carbon (DIC) concentration, and on-board nutrient analysis (see above). For more details on the incubation setup see Lipka et al. (2018).

Concentration changes of the bottom water solutes over time were evaluated via linear regression. Net flux rates (e.g., total oxygen uptake, TOU) were calculated using the slope (R² > 0.85, p < 0.1), separately for oxic (>89 μM O₂), hypoxic (< 89 μM O₂) and anoxic phases, corrected for the surface area to volume ration of each chamber. Because the effects of macrofauna activity under natural conditions were of interest, only oxic phase fluxes of O₂, NH₄, PO₄, and Si derived for the 21 cores incubated under oxic conditions were used in this study. Nutrient fluxes were calculated from linear increase or decrease of concentration vs. time, taking into account the volume of water above the incubated sediment.

Macrofaunal Data and Macrobenthic Parameters Tested

Macrozoobenthic samples were collected from short cores after completing the pore-water analysis or incubation. Thus, our data represents 2 sets: PW-dataset with 36 cores from 6 stations (Figure 1) and 5 sampling campaigns (March/April 2014, September/October 2014, April 2015, August/September 2015, and January 2016) analyzed for pore-water nutrients and diffusive fluxes; and INC-dataset with 41 cores incubated on board to estimate the total oxygen uptake and total benthic fluxes of NH₄, PO₄, and Si (of which only 21 are included in the evaluation) from 5 stations sampled during 3 latter cruises (see Table 1). The area covered by each core (10 cm inner diameter) was 0.00785 m². For macrofauna sampling sediment from the cores was sieved using a 1.0 mm sieve mesh size and samples were preserved in 4% buffered formaldehyde–seawater solution. In the laboratory, the organisms were sorted, identified to species level (with the exception of genus Phoronis and family Naididae), counted and weighted. The nomenclature was cross-checked following the World Register of Marine Species (WoRMS Editorial Board, 2018). Abundance and biomass data were standardized to an area of 1 m². Ash-free dry weight (AFDW) biomass (used here in the analysis) was estimated from the wet weight (determined including calcareous structures, without tubes) using species-specific conversion factors from the in-house list of the Leibniz Institute for Baltic Sea Research, Warnemünde (includes over 500 species, based on 15,000 measurements after drying at 60°C for 24 h and ignition at 500°C for 5 h).

Additionally, at every station and sampling event, at least 3 replicate samples were collected by the Van Veen grab (sampling area 0.1 m²). The samples were processed in the same way as described above for the short cores. Considering the patchy distribution of benthic assemblages (especially at finer scales, even in homogenous habitats) and in order to keep the amount of analysis feasible, this study is mainly focused on the collation of macrofaunal and biogeochemical data obtained directly from the same short cores. Therefore, the data from grab samples (referred to as VV) was used here only occasionally for comparative purposes (some more details on this data can be found in Gogina et al., 2017).

General descriptors of macrofaunal community considered in the analysis were species richness, total macrofaunal biomass and abundance. The metric of bioturbation potential (Solan et al., 2004) was used as proxy for sediment mixing activity:

B P_{c} = \sum_{i = 1}^{n} B P_{p i}, where B P_{p i} = {(B_{i} / A_{i})}^{0.5} * A_{i} * M_{i} * R_{i}

Where B_i is AFDW biomass in g/m² and A_i is abundance (in ind/m²) of the ith species, respectively. M_i and R_i are categorical scores reflecting increasing mobility and sediment reworking, respectively (indicating for M_i: 1—living in a fixed tube, 2—limited movement, 3—slow, free movement through the sediment, 4—free three dimensional movement via burrow system; for R_i: 1—epifauna that bioturbate at the sediment–water interface, 2—surficial modifiers, 3—upward and downward conveyors, 4—biodiffusors, 5—regenerators that excavate holes, transferring sediment at depth to the surface). M_i and R_i scores applied in this study are published in Gogina et al. (2017). BP_pi summed across all species in a sample gives an estimate of community-level bioturbation potential, BP_c

Quantitative services for ecosystem functioning provided by macrofauna vary mainly due to differences in biomass and burial depth between communities. Vertically resolved information on the distribution of organisms at various sediment depths in the PW and INC cores was lacking (but see Morys et al., 2017 for high resolution macrofauna depth distribution in spring 2014). In order to estimate the differences in the potential mixing depth for analyzed cores, an empirical equation that provides an estimate of biogenic mixing depth (BMD) based on the BPc published in Solan et al. (2004), i.e., BMD = 4.55 + 0.719*logBPc, was used. It should be noted that this metric is only a proxy, and is therefore not a precise match with the real sediment depth at which individuals were occurring.

To express the community functional structure in the analyzed cores, we used the biological traits approach (BTA) with a selected set of traits relevant for actuating flows of chemical substances, particularly fluxes of solutes, across the SWI and generally important for species-environmental relationships (Bremner et al., 2006; van der Linden et al., 2017). Table 2 provides the full list of 8 traits and 32 modalities considered. Based on the affinity of species to a particular modality of each trait, partial BTA scores, ranging from 0 (no affinity) to 1 (exclusive affinity), are assigned for each species so that the sum of scores of all modalities of each trait would sum to one. This “fuzzy coding” approach allows each taxon and each biological trait to have the same weight in further analysis (Chevenet et al., 1994; Bremner et al., 2006; van der Linden et al., 2017). To code the taxa by their functional traits we have used scores published in Gogina et al. (2014).

TABLE 2

Table 2. Biological traits used to describe functional structure of macrofauna community.

Accompanying Environmental Data

Environmental characteristics (Table 1) were measured at each sampling event parallel to the collection of cores. Water depth, near bottom water salinity and temperature were obtained from ship CTD. Near bottom oxygen concentration was derived using Winkler titration. One sample was taken at each sampling event for the determination of sediment characteristics (organic content, fraction finer 63 μm, median grain size, sorting, skewness). Chlorophyll a content of top 0.5 cm of sediment was determined as described in Morys et al. (2016). Sediment water contents of the freezed (under −20°C) and vacuum-dried samples were calculated gravimetrically by Bunke (2018).

Data Interpretation and Statistical Analysis

Pore-water profiles for separate cores were visually compared against the BPc values based on animals found in particular PW-cores, which represented a proxy for potential faunal activity.

Results of Shapiro-Wilk test of normality indicated that the null hypothesis of data originating from a normal distribution could be rejected for all but three variables (total NH₄ flux, species richness, and near-bottom oxygen concentrations). Therefore, only non-parametric methods could be allied to test for significant differences in geochemical, environmental and macrozoobenthic parameters between sediment types, sampling sites, and seasons.

Correlations between macrofaunal parameters such as the number of species, abundance, biomass, BPc, and pore-water reservoirs and solute fluxes at SWI were tested by Spearman correlation coefficients with p-value approximated by permutation. A significance level threshold of p ≤ 0.05 was assumed for all the analyses. Software PRIMER 7 (Clarke et al., 2014) and PERMANOVA+ (Anderson et al., 2008) were used to test the differences of functional traits structure between stations (settings: unrestricted permutation of raw data, number of permutations = 9999, estimation of P-values based on Monte Carlo simulations to overcome the problem of low number of replicates).

Fuzzy Correspondence Analysis (FCA, Chevenet et al., 1994) was applied to ordinate and explore biological traits structure (see section Macrofaunal data and macrobenthic parameters tested) in PW cores, using R-package “ade4” (Dray and Dufour, 2007). FCA is able to incorporate discrete fuzzy coded biological traits variables and apply corrections to account for the fact that each biological trait has multiple trait modalities that sum to no more than one. FCA ordination of cores and traits was obtained using Euclidean distances from relative frequencies of traits weighted by AFDW biomass. The result is an ordination plot where each point represents the biomass-weighted trait composition of each core, and proximity of points on the plot corresponds to structural similarity. FCA allows extraction of variability covered by each ordination axis and correlation ratios of each trait with principal axes.

To determine the set of environmental abiotic and biological variables (that reflect faunal metabolic and bioturbation activity) that best explain the multivariate variation of oxygen consumption and total benthic fluxes measured in incubation experiments, we performed a distance-based linear model permutation test (DistLM with dbRDA, McArdle and Anderson, 2001) employing the routine from the PRIMER 7 with PERMANOVA+ add-on. A similarity matrix based on the Euclidean distance of normalized total multivariate fluxes (total oxygen uptake and total benthic fluxes of NH₄, PO₄, and Si) data was obtained for 21 incubation cores. Highly correlated predictor variables were removed from the analysis (correlation coefficient ≥ 0.90), retaining: sediment grain size fraction below 63 μm, CTD derived near-bottom water salinity, oxygen concentrations, temperature, BPc, biomass (AFDW) of upward and downward conveyors, biomass of biodiffusers and biomass of suspension feeders and species richness (Nsp) of fauna colonizing the particular core (see Appendix 4). Step-wise selection criterion procedure based on 9,999 random permutations and the adjusted R² was used to determine the best set of predictors (similar to the approach used in Belley and Snelgrove, 2016). Variables showing significance levels above the threshold (p > 0.05) in the marginal test were excluded from the final model. Scatter plots of macrofaunal BPc against total solute fluxes were produced to aid interpretation of results.

To provide clear definitions, “diffusive fluxes” here are referred to the (theoretical) proportion of the flux that is only driven by molecular movements without any rearrangement of the solvent. This flux can only be reliably measured if advective processes can be largely excluded because advective fluxes have a great influence on the concentration gradients and therefore on the calculated diffusive fluxes. “Total fluxes,” estimated from incubation experiments, and denote the fluxes of dissolved chemical species by both advection of water (e.g., by organisms, via hydro-irrigation or due to relocation of pore-water and particles in surface sediments by human activities) and simple molecular diffusion across sediment-water interface. Total fluxes are the fluxes measured by incubation or estimated from process rates (like oxygen consumption) under the assumption that a steady state prevails. “Advective fluxes” in this study denote the transport of dissolved materials across sediment-water interface due to the transport of water by non-diffusive processes such as bioirrigation and bioturbation by burrowing organisms, hydro-irrigation and antropogenic activities like bottom trawling. Here “advective fluxes” are estimated as difference of total and diffusive fluxes. However, this estimate of the advective flux component using the difference of total and diffusive fluxes is only feasible for parameters that are not much involved in secondary reactions such as redox-reactions. An example of a reaction that excludes estimates of advective component of phosphate fluxes in our data is the reoxidation of Fe at the SWI resulting in the diffusive flux directed not toward the water column, but into the oxic zone of Fe-oxide precipitation, located somewhere in the top millimeters of sediment (Lipka et al., 2018). Similar secondary processes might occur with NH₄ that seems to be consumed in the oxic bottom waters during incubation experiments (see e.g., Figure 6B). Thus, out of the parameters covered with total fluxes measurements in our study arguably only Si seems suitable to separately estimate the advective flux. Our intention was to check if the effect of macrofaunal activity (expressed by BPc, abundance and biomass) on advective component of the Si flux is more distinct compared to that on the total or diffusive Si fluxes. As diffusive and total fluxes were derived from different cores (with high variability between the parallels), for this estimate the only option was the use of averages of the two methods for each site and sampling campaign.

Results

Pore-Water Profiles and Faunal Bioturbation Activity (PW Cores)

In total, 1,151 individuals representing 42 taxa (excluding the records of uncountable species—Bryozoa and Hydrozoa—from two OB, two LB, and one ST core from spring 2014 and autumn 2015) were identified in 36 PW cores (Appendix 1). Overall, Gastropoda were the most abundant, particularly at sandy OB and silty TW, followed by Bivalvia (most abundant at sandy ST, muddy AB, and MB sites) and Polychaeta (most abundant at muddy LB and second-most abundant at TW and ST). Polychaeta was also the most diverse group (represented by 20 taxa), followed by Bivalvia (8 species) and Gastropoda (3 species only).

To explore the effect of macrofaunal activity on the dynamic of nutrients in pore-waters we have compared the vertical concentration profiles of the cores with the lowest and highest BPc values within our PW-dataset. Pore-water concentration profiles of NH₄, P, H₂S, Fe, and H₄SiO₄ were plotted against the corresponding BPc values in Figures 2, 3. Generally, cores with higher BPc values display pore water concentration profiles with declining gradients toward SWI and concentrations that are closer to bottom water values, and therefore lower benthic solute reservoirs in the top centimeters. In typical nutrient-poor sandy sediments of OB site (Figure 3) cores indicating higher BPc showed deeper throughout mixing and more peaks, suggesting more intensive mixing, and active transport mechanism. This is visible in the two parallel cores OB2 and OB3 from EMB100 (Figure 3). Core OB3 indicated an overall second highest BPc (3,122) in PW-dataset (exceeded only slightly by the BPc value for core OB5 collected in January 2016, 3,288), with larger (>5 mm) individuals of deeper-dwelling bivalve Mya arenaria contributing most to this value.

FIGURE 2

Figure 2. Pore-water profiles of ${NH}_{4}^{+}$ , P, H₂S, Fe²⁺, and H₄SiO₄ for sediment cores collected from mud stations plotted against the BPc values defined by the colonizing benthic macrofauna community. Note the differences in axes scales. Cores are labeled according to sampling station (first 2 characters), followed by one character indicating core ID number, whereas the last 6 characters represent the cruise ID (respective seasons are indicated in Table 1).

FIGURE 3

Figure 3. Pore-water profiles of ${NH}_{4}^{+}$ , P, H₂S, Fe²⁺, and H₄SiO₄ for sediment cores collected from silt and sand stations plotted against the BPc values defined by the colonizing benthic macrofauna community. Note the differences in axes scales. Cores are labeled according to sampling station (first 2 characters), followed by one character indicating core ID number, whereas the last 6 characters represent the cruise ID (respective seasons are indicated in Table 1).

A similar pattern was evidenced for pore-water cores from muddy stations. From three parallel AB cores from ALK434 cruise AB3 had the lowest of three BPc values (110, determined by only 3 single individuals of Diastylis rathkei, Scoloplos armiger, and Nephtys hombergii with total biomass of 2.5 mg) and showed near-surface peaks in P, Fe²⁺,and H₄SiO₄ profiles (Figure 2). On the contrary, AB1 and AB2 displayed the highest BPc values recorded at this station (991, with the highest contribution of 10 individuals of Limecola balthica and 10 individuals of Scoloplos armiger, and 744 mostly shaped by large single individual of Nephtys ciliata and 5 L. balthica, respectively), corresponding with lower P, H₄SiO₄ and NH₄ reservoirs. AB1 showed no peaks in P, Fe²⁺, and H₄SiO₄ profiles, whereas peaks in AB2 occurred in deeper sediment layers and were less pronounced compared to those in AB3. In the example of two parallel cores from the sandy ST site collected in spring 2015 (ST4_EMB100 and ST5_EMB100, Figure 3), despite the higher BPc in core ST4 (where surficial modifiers Arctica islandica and L. balthica contributed most), the presence of one individual of N. hombergii in core ST5 seemed to have a greater effect on the Si and P profiles.

However, this relation is not always preserved. Generally, at muddy stations MB and LB biomass was dominated by Cumacea followed by Polychaeta. Sediments from two parallel MB cores from EMB100 cruise were nearly devoid of macrofauna (only one organism of surficial biodiffusor Diastylis rathkei was found in core MB2), and thus biogenic mixing by macrofauna is not expected here (Figure 2). Yet, despite the occurrence of peak and high concentrations of Fe in MB3 profile close to the sediment surface, detection of H₂S only below 10 cm sediment depth suggests possible impact of disturbance of sediment surface (for example by bottom trawling, known to be common for this area).

Sulfide Appearance Depth vs. Biogenic Mixing Depth

Biogenic mixing depth BMD, calculated based on BPc values of the community inhabiting each core, was used as proxy for macrofaunal burial depth. It was highest in sands and lowest in muds, with silt cores showing intermediate values (Figure 4). BMD significantly and positively correlated to sulfide appearance depths (R = 0.65, p < 0.01), suggesting the direct impact of faunal mixing. Cores from sandy stations also showed the largest variability in observed sulfide appearance depths values.

FIGURE 4

Figure 4. Biogenic mixing depth (BMD) values derived based on bioturbation potential (BPc) of the community inhabiting each core after Solan et al. (2004) plotted against sulfide appearance depths values (determined from pore-water profiles).

Functional Structure and Benthic Solute Reservoirs

The functional structure of macrofauna inhabiting the cores from the PW-dataset corresponded to changes in sediment type with sandy stations clustering on the left part of the ordination plot and muddy and silt stations located more to the right along the first FAC axis (Figure 5, Table 3). PERMANOVA analysis confirmed that differences in biomass-weighted trait data were significant between mud and sand cores (t = 4.02, P_(perm) = 0.001) as well as between sand and silt cores (t = 2.07, P_(perm) = 0.029), but not between mud and silt cores (t = 0.18, P_(perm) = 0.957). The co-structure between species biomass weighted trait-by-station arrays and pore-water solute reservoirs data set was not random [Monte-Carlo test: coinertia = 0.238 (coefficient of correlation between the two tables); p = 0.004].

FIGURE 5

Figure 5. Results of Fuzzy Correspondence Analysis on macrofauna data for pore-water (PW) cores. 2D ordination of samples based on trait categories exhibited by the communities (weighted by the untransformed biomass); samples (points) are grouped (circles) for stations, arrows point to the center of gravity of the samples presenting that category modality. Ellipses show the 95% confidence intervals. First axis explained 35% of total variation of the trait data, second axis−16%. See Table 2 for traits abbreviations.

TABLE 3

Table 3. Correlation ratios of biological traits for the FCA ordination, that corresponds to the proportion of variance on each axis that can be explained by the trait modalities (abbreviation as in Table 2) within each trait.

The sediment gradient from sand to mud was associated with a change from suspension to deposit feeding and a declining role of species penetrating below the oxygenated zone as well as species with large individual sizes. It corresponded to the increase of benthic solute reservoirs of ammonium, phosphate and silica in the top 10 cm of sediment (confirmed by significant positive correlation with FCA1, Table 3, see also Appendix 3), and the respective diffusive fluxes out of the sediments (supported by significant negative correlations with FCA1, Table 3). The FCA2 axis displayed only a weak significant correlation with Fe reservoir significantly declined along with FCA2 together with biomass of semipelagic taxa, taxa creating permanent burrows, and those contributing to deep-to-surface sediment transport.

Incubation (INC) Cores

A total of 21 incubations under oxic conditions spanned 5 sites and 3 seasons. In total, we identified 1,584 individuals representing 30 taxa (Appendix 1, Data Sheet 2). Overall in INC cores the most diverse class was Polychaeta, Bivalvia dominated the biomass, whereas Gastropoda were the most abundant. Peringia ulvae followed by Mya arenaria and Cerastoderma glaucum occurred in the highest densities in OB cores. Biomass was dominated by C. glaucum, L. balthica, and M. arenaria. At our sandy stations, functional richness, species richness, and density were significantly higher than in the adjacent muddy sediments (independent 2-group Mann-Whitney U-Test based on 33 PW cores, excluding TW cores: W = 24, W = 22.5 and W = 11, respectively, p-value ≤ 0.001).

PERMANOVA performed on Bray-Curtis similarity matrices of log (x+1) transformed AFDW data with pacifier variables (to account for cores without macrofauna), indicated significant differences in macrofauna composition among the 5 sampling sites and thus significantly greater variability in assemblages among sites than within sites (P_(perm) < 0.01). Pair-wise comparisons showed significantly (P_(perm) < 0.05) different benthic communities for all but three pairs of sites (AB and MB, P_(perm) = 0.421; LB and MB, P_(perm) = 0.196; LB and ST, P_(perm) = 0.177). Results of the Monte-Carlo permutation test further indicated distinct statistically significant differences of macrofauna composition expressed in AFDW biomass only between the OB and the other 4 stations.

PERMANOVA performed on Euclidean similarity matrix of normalized benthic total fluxes data indicated significant differences in benthic fluxes among the sampling sites and also significantly greater variability in flux among sites than within sites [P_(perm) < 0.05, Table 1]. Pair-wise comparisons with the Monte-Carlo test showed that benthic fluxes at LB differed significantly from fluxes measured at AB and ST [AB, P_(perm) = 0.019; ST, P_(perm) = 0.039], but no significant difference was found with fluxes at MB or OB. Pair-wise comparisons also showed that benthic fluxes at MB differed significantly from fluxes measured at AB and ST, that were also significantly different from each other [MB and AB, P_(perm) = 0.025; MB and ST, P_(perm) = 0.03; AB and ST, P_(perm) = 0.007] (data displayed in Figure 6).

FIGURE 6

Figure 6. Total oxygen uptake (A), NH₄ (B), Si (C) and PO₄ (D) total solute fluxes [μmol m²/day] across the SWI (negative values indicate release from the sediment) derived from incubation of natural cores (Lipka et al., 2018), plotted against community bioturbation potential BPc. Different colors indicate different stations, different symbols show seasonal variability.

Analysis of outliers performed for each station separately highlighted that core OB1_spr15 showed extremely high BPc, total AFDW as well as AFDW of bivalves and species attributed to biodiffusers (R4), which in turn associated with extremely high values of oxygen uptake and release of PO₄, NH₄, and Si. Results of DistLM obtained on a dataset excluding outliers showed that log-transformed BPc, species richness (Nsp) and biomass of biodiffusers were the best proxies among the tested set of biotic and abiotic parameters and could explain 63.6% of the multivariate total benthic flux variations (Appendix 4).

For a set of 21 incubated cores, the only significant values of Spearman rank correlation found for benthic fluxes and macrofaunal parameters were between total Si flux and BPc (r = 0.470, p = 0.032) and between Si flux and species richness Nsp (r = −0.442, p = 0.045). Correlation between NH₄ flux and AFDW of biodiffusors (R4afdw) was negative and close to significance threshold (r = −0.423, p = 0.056). No consistent relationship between ammonium fluxes measured in INC-dataset and the presence and biomass of infauna creating permanent burrows was found. For two parallel cores from the OB station sampled in spring 2015, the two overall highest values of NH₄ release were measured. In this case, increase of ammonium release corresponded with increased BPc observed in the cores (Figure 6B). Despite the presence of burrowing invertebrates and the second-high BPc values in both OB cores from winter 2016, only one of them displayed significant release of ammonium from the sediment, whereas for the other incubated core, no significant flux across the SWI could be detected. Moreover, even though burrowing macrofauna was present at OB in autumn 2015, sediment uptake of ammonium from the overlaying water was measured during the incubation. This was also the case at AB station in spring 2015 and winter 2016, though biomass of burrowing infauna (or species scored as biodiffusors) in AB cores was comparatively low (Appendix 5).

Estimate of Advective Flux Component and Effect of Macrofaunal Activity

The results of comparing the effects of macrofaunal activity on the total, diffusive and advective components of the Si flux are presented in Figure 7 and Table 4. As we expected, the most pronounced effect was the positive correlation between advective flux component and infauna activity expressed by bioturbation potential (and, to a smaller degree, abundance and biomass), while the diffusive flux negatively correlated with benthos activity. When muddy and sandy sites were considered separately, no significant effects were found for muddy stations. Looking into more detail, at muddy LB site in autumn 2015 only 65% of Si flux was due to advection and it corresponded with very low bioturbation potential of infauna enclosed in INC cores. In contrast, during winter 2016, diffusive Si fluxes were close to 0, therefore advection was nearly 100%, and BPc in INC cores was comparatively high (mostly due to biomass of bivalve A. islandica that could mix surficial sediments). At MB station in spring 2015 high total and low diffusive Si fluxes resulted in rather high advective flux component (75%), but the lowest macrozoobenthos activity was recorded in INC and PW cores (however, in Van Veen garb samples from MB site collected on same occasion AFDW biomass was rather high, mainly again due to A. islandica exceeding 10 g/m²). For sandy stations alone, advective component displayed correlation only with biomass averaged across incubated cores, and this relation was negative.

FIGURE 7

Figure 7. Estimate of advective component (in %) of Si total fluxes at different sites and seasons (EMB100–spring 2015, EMB111–autumn 2015, MSM50–winter 2016) plotted against log-transformed major macrofaunal parameters (BPc_VV, BPc_PW, BPc_INC—community bioturbation potential estimated per m² based on Van Veen Grab (0.1 m²) samples (VV), pore-water (PW) cores set and incubation (INC) cores set, respectively, Abu refers to respective abundance [ind/m²] and AFDW—to biomass [AFDW g/m²].

TABLE 4

Table 4. Spearman rank correlation coefficients between the estimated total, diffusive, and advective fluxes, as well as advective component of Si, and general macrofaunal parameters (averaged per cruise and site).

Discussion

Previous studies have established that macrofaunal reworking, including bioturbation, can cause much higher nutrient flux from sediments into the water column than molecular diffusion in coastal areas (Mortimer et al., 1999; Christensen et al., 2000). Bioirrigation generates transport of pore-water especially in the presence of tube-dwelling animals or gallery-building biodiffusers (Aller, 1982; Kristensen, 2001; Shull et al., 2009; Braeckman et al., 2010). Where macrofauna is abundant, oxygen is introduced punctually deeper into the sediment surface, therefore shifting the sulfidic zone downwards. Although this is recognized, actual quantitative field estimates are scarce.

It is suggested that functional group composition of biota is the most critical aspect of biodiversity for ecosystem functioning in the context of global biogeochemical cycles (Snelgrove et al., 2018). We attempted to identify the main players responsible for altered geochemical profiles and looked for similarities and patterns explaining significant alterations of calculated fluxes across the SWI in the south-western Baltic Sea. Taxonomic and functional identities, estimates of abundance, biomass and community bioturbation potential were used to assess the role of macrofauna in variation of pore-water reservoirs, oxygen uptake and total nutrient efflux. The first general observation is that data scatters a lot over time, differences in scale and from repetitions.

Among the possible reasons explaining why the interpretation of data reveals no clear strong dependencies are three major points: (1) the activity of organisms has stochastic nature, it is not constant and mixing, pumping, and respiration work of individuals captured in the analyzed cores may be reduced in intensity or ceased, especially under stressed conditions associated with sampling procedures; (2) human induced impacts on the seafloor such as trawling fishery evidenced e.g., in the Bay of Mecklenburg (that were not quantitatively captured by our analysis) may overwhelm the coupling of macrofauna activity and sedimentary biogeochemistry; (3) marine microorganisms are a major component of global nutrient cycles (Arrigo, 2005), whereas macrozoobenthos is only a mediator of bacterial-induced processes in benthic-pelagic coupling (Yazdani Foshtomi et al., 2015). Leuzinger et al. (2011) found that the magnitude of responses decline in higher order interactions and suggest that global scale impacts on the ecosystem may be dampened as a result of positive and negative drivers of change.

Parameters like salinity, temperature and oxygen affect both redox-reactions and animal distribution, density and behavior (Gogina et al., 2014; Gammal et al., 2017). During autumn 2015, our study observed oxygen depletion at all stations with the exception of OB, and near-bottom water temperature was much warmer compared to the other cruises (Table 1). Hypoxia could be stimulated by an inflow of North Sea waters into the Baltic Sea in December 2014 that terminated a decade of continuously increasing oxygen deficiency in the central deep basins (Mohrholz et al., 2015). According to findings of Braeckman et al. (2010) faunal simulation can double benthic microbial respiration and mineralization especially in summer, whereas in winter, before the spring phytoplankton bloom, faunal effects can be negligible. The authors explained the strong summer effect by reference to higher macrobenthic activity owing to the elevated temperature and better animals' conditions (enhanced metabolism). Perhaps due to the complexity of interactive effects, we were unable to find distinct seasonal patterns in our data. Indeed, at LB station we observed the highest total oxygen uptake in autumn 2015 among muddy sites, corresponding with one of the lowest BPc values, whereas at MB site intermediate total oxygen uptake was found in the core with the highest BPc (Figure 6). The absence of abiotic variables in the best solution of our DistLM model (to explain the multivariate patterns of variation of oxygen consumption and total benthic fluxes measured in incubation experiments, see Data sheet 2) could be due to the fact that in our case information on macrofauna was resolved to the core level whereas abiotic variables were measured at coarser spatial scale (per station and sampling event), that can explain the low predictive power of the latter in this limited dataset. On the other hand, it can also be due to the close linkage between species and their environment, that the former well mirroring the latter (in agreement with Belley and Snelgrove, 2016).

Estuarine ecosystems are subjected to a high degree of variability in their environmental conditions due to usually large areas with shallow water depth, freshwater and seawater inputs, as well as climatic and anthropogenic impacts. van der Linden et al. (2017) found polychaetes to display higher levels of functional diversity compared to mollusc group in tropical estuaries. Whereas infaunal activity inherently contributes to increased solute exchange between the sediment and the overlying water, the contribution of the underlying drivers varies among functional groups through increased diffusive fluxes, advective pore water bioirrigation (Kristensen et al., 2012), and animal excretion. Mermillod-Blondin et al. (2004) showed that in diffusion dominated systems Cerastoderma edule acted as a biodiffuser in the top 2 cm of the sediment, but had little effect on O₂ consumption, water exchange across SWI, microbial characteristics, and release of nutrients from the sediment. In contrast, Corophium volutator and Hediste diversicolor build burrows in the sediment that allows transport of surface particles into biogenic structures and enhanced solute exchange. H. diversicolor typically bury down to 10 cm, maximal—down to 30 cm sediment depth (Zettler et al., 1994). C. glaucum (having traits similar to C. edule) and H. diversicolor were also typically present in INC cores from our OB station (see Appendix 2), together with another burrowing polychaetes Marenzelleria viridis. Particularly in the core OB1, where overall highest nutrient fluxes released from the sediments were observed, H. diversicolor contributed up to 11% to total biomass (with 3.8 g AFDW/m²) and dominated the BPc (with a share of 24%), whereas C. glaucum had the share of 26% in biomass (with 9.1 g AFDW/m²) and 18% in BPc.

Norkko et al. (2013) depicted the number of large bivalves as the best predictor of ecosystem functioning. In our data, bivalves that can reach individual sizes exceeding 5 mm (at their longest dimension) are represented by A. islandica, C. glaucum, L. balthica and deep-dwelling (down to 25 cm sediment depth) M. arenaria. Increased biomass of A. islandica (e.g., ST cores from spring 2015 and winter 2016) was observed along with a declining release of ammonium and phosphate, whereas higher abundance of M. arenaria coincided with a significant increase of ammonium and phosphate total fluxes out of the sediment at OB site. Also, sulfide appearance depth showed a strong positive correlation with abundance of M. arenaria (see Appendix 6). At muddy stations, total phosphate release was recorded only in cores from the MB site where no bivalves were present (with the exception of 1 small individual of Corbula gibba having AFDW biomass of 0.22 mg). The only incubated core from MB station that displayed no phosphate release was inhabited by one relatively large (15 mg) individual of N. hombergii. This all suggests the importance of both tube dwelling and burrowing polychaetes and bivalves acting as biodiffusors.

The significant positive relationship between biogenic mixing depth sensu Solan et al. (2004) and sulfide appearance depth found in our study is in agreement with findings of Sturdivant and Shimizu (2017) who showed a match between the burrow depth and the aRPD values (apparent color redox-potential discontinuity). Our results also support the conclusions of the previous studies that found BPc to be a good proxy for ecosystem efficiency (e.g., Baldrighi et al., 2017). According to Cozzoli (2016), metabolic energy will be spent on sediment reworking independent of the type of bioturbation, in order to obtain food, resulting in sediment mixing and the exposure of fine particles to the surface, where they can be resuspended. The ability to spend this energy will scale with metabolism, and therefore it is natural to expect a scaling of ecosystem engineering with body size and abundance.

We have tested the effect of macrofaunal activity (expressed by abundance, biomass and BPc that includes both) on the advective component of the Si flux and compared it to the effect on the total and diffusive Si fluxes. In accordance with our expectations, the most pronounced effect was the positive correlation between the advective flux component and infauna activity expressed by bioturbation potential (and, to a smaller degree, abundance and biomass), while the diffusive flux negatively correlated with benthos activity, suggesting that gradients were leveled out by advective processes (see Appendix 7). In line with Huettel et al. (1998), advective transport processes were the most important in controlling total fluxes in permeable sediments. However, the significance of the relationship between advective component and macrofaunal activity in our data did not preserve when sandy and muddy sites were treated separately (for sandy cores, only the negative correlation with biomass was significant).

Based on this study and the high variability of data, it is difficult to predict implications that can be expected if the invertebrate community of those typical habitats of the south-western Baltic Sea would change (e.g., due to climatic variability or human-induced pressures). Generally, increased biodiversity enhances ecosystem functioning due to increased probability of the presence of species fulfilling particular functions (Magurran, 2012). Our study was ultimately constrained spatially (to our set of sampling stations) and temporally (by the onset of cruises with respective conditions of preceding inflow, seasonal hypoxia in autumn 2015 and/or storms). Both species distributions as well as nutrient fluxes are temporally dynamic (Figure 6) and those natural patterns inhibit clear interpretation of presented field data, i.e., this is an example of the need for more repeated samplings and replicates at these areas of high spatio-temporal variability there is a need of more repeated samplings and replicates (but see Bolam et al., 2002). Nevertheless, such data free from experimental manipulations (such as de-oxygenation or defaunation) presents valuable information on the true complexity of ecosystem functioning in the natural system (Carpenter, 1996). The description of the variability itself has a value as a source of documented information for the future generalizations (Coble et al., 2016). Until a great deal of empirical experience is collected, we do not know whether a given example of structure is extraordinary, or it is a mere trivial expression of something which we may learn to expect all the time (Hutchinson, 1953).

Overall, from the parameters tested, bioturbation potential was identified as the best proxy for variation in pore-water reservoirs of biogeochemical components compared to available abiotic proxies, abundance, or biomass parameters (in total or considered separately for different functional groups). This was anticipated due to high scoring of deep-burrowing and bioturbating taxa (particularly large individuals) that have a major influence on all oxygen-dependent biogeochemical processes, supporting the importance of functional groups (Bolam et al., 2002).

Small-scale disturbances such as hypoxia (observed during EMB111 cruise in autumn 2015), sediment transport due to currents, trawling or storms, may play a driving, often diminishing, role for biodiversity in general and bioturbation in particular (Villnäs et al., 2013; Gogina et al., 2017). These disturbances also affect fluxes in and out of the sediment (Villnäs et al., 2013). Whereas the mixing effect of physical disturbance can be similar to bioturbation, the effects of species respiration or enlarged pool of their metabolic waste products, such as ammonium (Thoms et al., 2018), is likely to be different. This highlights the complexity of ecosystem function relationships in the soft sediment benthos.

Another, more important factor that may impact the biogeochemical sediment services that has to be considered in future studies is the abundance, community structure and activity of microorganisms, that are particularly active in muddy and sandy surface sediments and cause the liberation of metabolites to the pore-waters (Llobet-Brossa et al., 1998, 2002; Sahm et al., 1999; de Beer et al., 2005), therefore creating the disequilibrium physicochemical environment that promotes element fluxes within the sediments and across the sediment-water interface (Boudreau, 1997). A further factor that may impact the transport processes of electrons in the surface sediment and may have to be addressed in future studies is the activity of cable bacteria (van de Velde et al., 2017).

Despite of our initial expectation, no clear relationships and consistent patterns were found between major macrobenthic parameters and biogeochemical fluxes at stations of the same sediment type. However, we should not ignore this absence of correlations even if it does not meet the general expectation, and therefore we find reporting of these negative results here particularly important. For example, Jacquot et al. (2018) also reported the lack of the relationship between macrofauna and denitrification they were aiming to find, and suggested that macrofaunal density and/or burrowing activity should be above some threshold to increase bioturbation or irrigation to sufficient levels to affect denitrification (see also Appendix 8).

Conclusions and Perspectives

Our data indicates that bioturbation and bioirrigation alter near-surface pore-water profiles of nutrients toward bottom water signals. Biogeochemical fluxes can partly be explained by functional structure of macrofauna, with community bioturbation potential performing as the best proxy among the tested set of biotic and abiotic variables. The effects of macrobenthos on benthic solute reservoirs and total fluxes differed between sediment types, specific locations, and seasons. High complexity and variability in spatial and temporal distribution, anthropogenic influence (e.g., fishery) and short-term physical events (e.g., storms) hinder interpretability of collected data. Thus, the outcomes of this case study are context dependent, but the search and exploration of similar parallel datasets combining characterization of biogeochemical processes across the sediment-water interface should be continued at different scales: in the study area as well as other regions. Ideally, microbiological data should be included in the analysis, to reveal significant relationships, address their generality and better capture the natural variability and various hidden mechanisms.

Author Contributions

The contributions of the authors are: MG, ML, JW, MB, and MZ designed the study. MG did macrofauna sampling, most of the analysis and calculations and wrote the paper. ML carried out sediment sampling and pore-water analysis. ML and JW conducted the incubation experiments, and ML performed the fluxes calculations. BL contributed to modeling the impact of physical disturbances on pore-water profiles. CM contributed the chlorophyll a data and was involved in sampling and analysis. All authors contributed to writing, review and editing, read the text and agreed about the final version of the manuscript.

Funding

This work was supported by the German Federal Ministry for Education and Research: KÜNO Projects SECOS (03F0666) and SECOS-Synthese (03F0738).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We greatly acknowledge the valuable support of all IOW colleagues deployed in field sampling and laboratory analyses, in particular, Anne Köhler and Iris Schmiedinger, and the captains and crews of RVs Elisabeth Mann Borgese, Poseidon, Alkor, and Maria S. Merian. We greatly acknowledge the comments by two reviewers and the editor Teresa Radziejewska that helped to considerably improve this manuscript. We express our thanks to Irina Steinberg for correcting the English.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2018.00489/full#supplementary-material

References

Aller, R. C. (1982). “The effects of macrobenthos on chemical properties of marine sediment and overlying water,” in Animal–Sediment Relations. Topics in Geobiology, Vol 100, eds P.L. McCall and M. J. S. Tevesz (Boston, MA: Springer), 53–102.

Google Scholar

Aller, R. C. (1994). Bioturbation and remineralization of sedimentary organic matter: effects of redox oscillation. Chem. Geol. 114, 331–345. doi: 10.1016/0009-2541(94)90062-0