Sediment Bulk Density Effects on Benthic Macrofauna Burrowing and Bioturbation Behavior

Introduction

Estuarine intertidal systems are among the most productive ecosystems globally and are responsible for important ecosystem services such as providing coastal protection, carbon sequestration, food production, recreation areas, as well as habitat and nursery grounds for fish and birds (Meynecke et al., 2008; Koch et al., 2009; Barbier et al., 2011; Seitz et al., 2014). At the margins of the land and sea, many factors, both biotic and abiotic, drive the development of these ecosystems (Gray and Elliott, 2009). Benthic macrofauna are a key component of intertidal systems: not only are these animals important prey for birds and fish (Piersma et al., 1993; Zwarts and Wanink, 1993; Bocher et al., 2014), as ecosystem engineers they modify their sedimentary environment (Bouma et al., 2009). Benthic macrofauna–sediment interactions vary considerably in space and time. They depend on species distribution and behavior as well as local environmental conditions such as sediment composition. The strength of these interactions, in turn, underpins the supply of ecosystem services delivered by intertidal ecosystems. Some sediment characteristics, such as grain size and mud content, have been highly investigated for correlation with species occurrence (Thrush et al., 2003; Pratt et al., 2014) and effects on benthic macrofauna behavior (Dorgan, 2015; McCartain et al., 2017). One sediment characteristic that does not have a well-known effect on benthic macrofauna and their behavior is sediment bulk density, an indicator of both sediment compaction and water content (Grabowski et al., 2011), which is the focus of this paper.

By being inversely related to sediment porosity, or the amount of water retained in a waterlogged sediment, bulk density influences the sediment oxygen content, chemistry, and organic matter (Gray and Elliott, 2009; Dowd et al., 2014). Thus sediment bulk density (dry sediment weight per sediment wet volume) is an important characteristic of intertidal geomorphology and describes a measure of sediment compaction which is missing from grain size gradients. The range of bulk density values vary on a regional and local scale. Bulk density tends to increase with grain size (Ysebaert et al., 2005) and sediment strength (Lucking et al., 2017), while it tends to decrease with silt content, erodibility, and organic matter content (Grabowski et al., 2011; Stringer et al., 2016; Joensuu et al., 2018). More compact sandy sediments typically have a bulk density of around 1–2 g cm^–3 and softer muddier sediments typically have a bulk density of around 0.2–1.5 g cm^–3 (Andersen et al., 2005; Grabowski et al., 2011; Stringer et al., 2016). However, sediments with similar grain size composition can have a range of bulk densities due to different water contents (e.g., Widdows et al., 2007; Soares and Sobral, 2009). At very low bulk densities, the sediment is viscous and more akin to a fluid (Grabowski et al., 2011). These conditions are found in many systems worldwide, in particular mangroves (Stringer et al., 2016) and intertidal mud flats (Walles et al., 2017). Furthermore, sediment bulk density tends to increase with erosion and decrease with sediment deposition (Dyer et al., 2000). Besides long-term erosion and sedimentation trends, intertidal flats experience short-term bed-level variation which differs between sites of contrasting wave exposure (Hu et al., 2017). These bed-level trends may contribute to the variation of sediment bulk density over several temporal and spatial scales in estuarine areas. Though bulk density is an important biogeomorphological characteristic of intertidal sediments, we still do not understand how significant this characteristic is for determining benthic macrofauna organism-sediment interactions.

Macrofauna’s mixing of sediment, or bioturbation, typically affects parameters like sediment permeability, grain size, and erodibility (Volkenborn et al., 2009; Kristensen et al., 2013; Harris et al., 2016), and also drives many biogeochemical processes of tidal flats (Kristensen, 1988; Gray and Elliott, 2009). Sediment mixing by benthic macrofauna increases sediment permeability, which plays a vital role in nutrient cycling by driving the circulation of oxygen and nutrients below the sediment surface (Aller, 1994; Thrush et al., 2006). Though benthic macrofauna can move laterally by crawling through the sediment [e.g., H. diversicolor (Aberson et al., 2011) and C. edule (Richardson et al., 1993)] animals usually avoid crawling at the sediment surface since it exposes them to predation (Ens et al., 1997). In this study we focus on vertical sediment mixing which is produced mainly by burrowing and deposit feeding (Kristensen et al., 2012). The mode of mixing is important for determining what kind of ecosystem services the organisms provide. Local mixing, or biodiffusion (Boudreau, 1986a), such as by C. edule, can discourage the build-up of fine mud (Montserrat et al., 2008) and oxygenate the top sediment layers (Mermillod-Blondin et al., 2004). Non-local mixing, where surface particles are transported from depth to the surface and vice versa (Boudreau, 1986b), such as by the polychaete H. diversicolor, transports particles faster and deeper in the sediment, which creates new habitat for deeper living organisms (Gray and Elliott, 2009), thus increasing biodiversity (Sturdivant and Shimizu, 2017). Moreover, non-local mixing also facilitates biogeochemical processes such as denitrification and permanent burial of pollutants (Mermillod-Blondin et al., 2004).

While macrofauna drive tidal flat biogeomorphology and biogeochemistry through their movement behavior, they also change their behavior under different sediment conditions. Many studies have shown that sediment properties determine macrofauna behavior, for example, Macomona liliana decreases its siphon movement in cohesive muddy sediment compared to sandier sediment (McCartain et al., 2017), which could have consequences for sediment permeability and oxygen flow. While the effects of some sediment characteristics, like sediment grain size, on behavior have been well-studied, bulk density effects have not. Studying how macrofauna change their behavior under different bulk densities will help better understand how sediment conditions affect the ecosystem functioning of benthic macrofauna and thus the biogeochemistry and biogeomorphology of tidal flats.

The objective of this study was to investigate how five dominant and functionally different bivalve and polychaete species modify their burrowing and bioturbation behavior under different bulk densities. We hypothesized (1) that both bivalves and polychaetes would have lower mobility (slower burrowing speed, shallower burrowing depth, less active bioturbation) in more compact sediments, (2) that bulk density would have a similar effect on benthic macrofauna behavior in both sandy and muddy sediment, and (3) that larger bivalves would be slower and less active in higher bulk density sediments than the smaller and younger ones, due to greater biomass hampering their mobility. To test our hypotheses, we performed a mesocosm experiment where we subjected three bivalve species (Limecola balthica, Scrobicularia plana, and C. edule) and two polychaete species (Hediste diversicolor and Arenicola marina) to four sediment treatments of varying bulk density and mud content. We measured the following indicators of movement behavior: burrowing speed, burrowing depth, bioturbation activity and respiration rates. We selected the species based on their dominant prevalence in the Scheldt intertidal (Cozzoli et al., 2013) as well as their functional differences. In particular, the five species covered the different bioturbation types (Table 1).

Materials and Methods

Experimental Design

We designed a mesocosm experiment in which three bivalve species (Limecola balthica, Scrobicularia plana, and Cerastoderma edule) and two polychaete species (Hediste diversicolor and Arenicola marina) were subjected to two crossed bulk density (dry sediment weight per sediment wet volume) and sediment grain size treatments (Figure 1b). Our four treatments were: compact-sandy (CS), soft-sandy (SS), compact-muddy (CM), and soft-muddy (SM) (see Table 2 for sediment characteristics). We used two sediment grain sizes to cover a larger part of the range of bulk densities observed in the Scheldt intertidal, as well as to test the interaction between sediment bulk density and grain size on animal behavior. We prepared the treatments from two sediments collected from intertidal flats in the Eastern Scheldt estuary (Figure 1a). The sandier sediment was collected at the Oesterdam (51.466700 and 4.221389) and the muddier sediment was collected from Prosperhaven (51.490305 and 4.259167) nearby. Both were collected several months before the start of the experiment and were stored in outdoor closed bins. The sandier sediment was drained and passed through a large (5 mm) sieve before use. The sticky consistency of the muddy sediment did not allow it to be passed through a sieve, therefore we removed the larger fragments, such as shells, by hand. We incorporated seawater (amount equaling 7% of the sediment weight for compact high bulk density treatments and 15% for soft low bulk density treatments) into the sifted sediments with a standing industrial mixer the day before we added the animals (see Table 2 for resulting sediment characteristics). The sediment was mixed in batches of 5 kg for at least 5 min to ensure the mixture was homogenous. The sediment was placed in pots made from a sawed-off PVC pipe (height 12 cm, diameter 11.5 cm). Each PVC pipe was capped off on the bottom with a removable plastic cap and lined with a plastic bag to prevent water loss. We encircled the brims of the pots with mesh (1 mm mesh size) that extended 2 cm above the high tide height to keep the animals contained to their unit. The pots were filled to 0.5 cm below the brim and placed in the mesocosm tanks.

The experiment was conducted in ten tidal tanks in a climate-controlled room where the water was kept at 18°C, the same temperature as the Eastern Scheldt in July 2018 when the experiment was carried out. Each tidal tank was composed of two 1.2 m by 0.8 m tanks stacked on top of each other (see Figure 1c for example). Unfiltered water from the Eastern Scheldt estuary was pumped from the bottom tank up to the top tank to simulate tidal conditions. High tide conditions (5 cm water above experimental units) lasted 6 h and occurred twice a day. We changed the water once a week, and in addition to the nutrients contained in the raw Eastern Scheldt water, we fed the animals with an algal concentrate (Shellfish Diet from Reed Mariculture) 5 mL per tank twice per week.

We had eight experimental blocks composed of two polychaete species, three bivalve species with two size classes for S. plana and C. edule, and controls without animals (Table 1). We collected the animals at three sites in the Eastern Scheldt: A. marina were collected at the Oesterdam (51.46670 and 4.22139), C. edule and L. balthica were collected at Dortsman (51.543678 and 4.055841), and S. plana and H. diversicolor were collected in Yerseke (51.489245 and 4.057288), 3 days to 24 h before the start of the experiment and stored in a tidal tank in the mesocosm (Figure 1a). Each block was replicated six times for each sediment treatment. Thus, we had 24 units for each species/size block, resulting in a total experiment of 8 × 24 = 192 units (see details in Table 1).

All the animals used occurred naturally in the two sediments we used for the experiment (see Cozzoli et al., 2013 for biomass probability distribution in the Scheldt depending on grain size for C. edule, L. balthica, H. diversicolor, and A. marina; see probability distribution for S. plana depending on grain size in Van Colen et al. (2014). A. marina was the only species placed in unnatural conditions (muddy sediment) not found in the field. Not unexpectedly, A. marina had low survival in the muddy sediment treatment, especially the compact-muddy treatment. There was some low mortality for H. diversicolor (see Results section “Survival Rate and Growth”) but this was not linked to any sediment treatment. We recorded no mortality in the bivalves.

The units were divided randomly into the ten tidal tanks (24 per tank) to minimize a tank effect. We used realistic but low animal densities to avoid competition for space and food resources. The realistic densities also allowed us to obtain realistic behavioral responses from the animals. Thus, the number of animals per core varied between species, which means that responses and animal abundances were colinear. The number of animals per pot were one S. plana, two C. edule, five L. balthica, five A. marina, and three H. diversicolor. We only compared the magnitude of the species effect between size class treatments (i.e., small and large C. edule and S. plana). Though the abundances varied per species, the biomass per pot was roughly similar.

Sediment Characterization

We prepared 24 extra cores to take bulk density and sediment penetration resistance measurements to determine whether the sediment treatments changed over the course of the experiment. The cores were placed in a separate tidal tank in the same climate room as the experiment. We used a universal testing machine (Instron, Baltimore, MD, United States) to measure sediment penetration resistance, or the amount of pressure necessary to compress the sediment. The universal testing machine measures the load as a crosshead penetrates the sediment and extends to the bottom of the sediment core at a constant speed (Bokuniewicz et al., 1975). To measure bulk density (dry sediment weight per sediment wet volume), we determined the water content of a fixed volume of sediment which was collected using a cut-off syringe and weighed before and after freeze-drying for 48 h. We took penetration resistance and bulk density measurements of the different sediment treatments after 12 h, 1 week, 2 weeks, and at the end of the experiment. The four sets of cores were only used for each set of measurements once because the sampling for bulk density destroys the core.

While the penetration resistance changed over time, there was no overlap in penetration resistance profiles between the four treatments by the end of the experiment (Figure 2). The bulk densities of the four treatments remained significantly distinct from each other over the course of the experiment, F (3,80) = 649, p < 0.001, and post hoc Tukey test p < 0.008 for all treatment comparisons. In addition, the inclusion of the measurement date in the linear model evaluating treatment effects on sediment bulk density did not improve the model fit, F (1,79) = 0.001, p = 0.97, indicating that the sediment bulk densities did not change between the beginning and end of the experiment.

Process Measurements: Burrowing Rates

Burrowing rates of bivalves were monitored by counting the animals still present at the sediment surface shortly after the start of the experiment and then every 12 h until they were all buried, after 2 days. This method worked well for the bivalves, but not for the fast burrowing polychaetes as they burrowed almost immediately. Hence, we performed an additional short experiment to investigate how bulk density impacted the burrowing ability of polychaetes. We prepared four sediment treatments with different bulk densities from the sandy sediment by adding in seawater at 2, 10, 15, and 25% of the sediment weight and homogenizing with the standing cement mixer. The four sediment treatments had average bulk densities of: 1.35, 1.31, 1.26, and 1.22 g cm^–3. We then added the sediment to empty pots in the same way as we did for the cores in the mesocosm. We placed six A. marina or six H. diversicolor on the sediment surface and recorded how many animals were still visible at thirty second intervals. We repeated the experiment twice for each species/sediment combination.

Process Measurements: Metabolic Activity

We performed oxygen incubation experiments to test whether different sediment characteristics affected macrofaunal metabolic activity. The incubation experiments were run 24 days after the experiment had begun so that the animals had acclimated to the different sediment conditions. We performed the oxygen consumption experiments with three replicates each for: the controls, H. diversicolor, L. balthica, both size classes of S. plana, and the large C. edule. We could not perform the oxygen consumption experiments on all species blocks because of equipment limitations. The experiment consisted of a water bath containing capped PVC tubes filled with Eastern Scheldt water. We placed each core into the larger PVC tube. The cores descended into the tube slowly and the surfaces of the sediment were not disturbed. We oxygenated the water to raise the oxygen content until it was near saturation and sealed the cores with a custom cap. The cap included a stirrer to homogenize the entire water column’s oxygen content. We then measured the amount of oxygen in the water over time using a FireSting sensor (see Braeckman et al., 2014 for oxygen incubation method). The experiments ran for at least 5 h and we did not allow the oxygen level to descend below 60%. There was a possibility that differences in respiration rates between cores could have been due to animal death, however, all the bivalves used in the oxygen incubation experiments were alive at the end of the 29 days. H. diversicolor experienced low mortality, but without a sediment treatment effect.

Process Measurements: Sediment Mixing

To better understand how the sediment treatments would affect animal bioturbation activity and the resulting sediment mixing, we added luminophores (Environmental Tracing Systems, United Kingdom) to the sediment cores, which are inert natural sediment particles dyed with luminescent paint used to track bioturbation (Gerino et al., 2003; Solan et al., 2004). Luminophores are widely used as a non-toxic tracer to study sediment mixing by benthic animals as small amounts do not affect the animal behavior and do not harm the animals (Maire et al., 2008). We made frozen sediment disks following our standard procedure for mixing the sediment treatments but replaced ten percent of the sifted sediment with luminophores. We added one 0.5 cm thick disk to the top of all cores after the animal had burrowed in to avoid recording the initial burrowing movement as bioturbation activity. Some of the large S. plana never burrowed and so we added the luminophores on top of them.

After 29 days, we processed the experiment. We sliced the cores once lengthwise to obtain a vertical profile of luminophore incorporation into the sediment. These profiles show the total amount of mixing that was done by the animal over the course of the experiment, thus representing the time-integrated outcome. As we were not using the luminophores to model bioturbation intensity, a single lengthwise slice was sufficient for our purposes of investigating the sediment treatment effect on bioturbation activity and mixing mode. In addition, the lengthwise slicing allowed us to recover the animals intact for subsequent physiological measurements. As a consequence of the vertical slicing, we may in some cases have underestimated bioturbation activity if a burrow did not traverse our slice. However, we compensated for possible underestimation by averaging the profiles over six cores for each species-sediment combination. This allowed us to obtain a general estimate of the bioturbation patterns.

To count the luminophores as a function of depth we photographed the two halves of the core under a blacklight using a digital mirror-reflex camera (Canon EOS 1100d) attached to a tripod. The pictures (2848 × 4247 pixels) were analyzed using a custom ImageJ (Fiji) script. In this process, the red layer of an RGB filter was used to highlight the magenta colored luminophores. The brightness value (125–255) was used to select luminophores from other pixels. We estimated bioturbation activity by counting the luminophore pixels in the pictures of the halved cores by 0.5 cm layer, or “bin.” The bins had an area of 10.5 × 0.5 cm, which corresponded to 1722 × 82 pixels. Because one luminophore pixel does not necessarily correspond to a single luminophore grain, we use “luminophore pixels” instead of “luminophore count” as the unit for luminophore quantity throughout the text. The edges of the cores (0.5 cm from the sides and the bottom 1.5 cm) were excluded to avoid skewing our estimates because the plastic lining caused luminophores to accumulate at the sides and bottom of the cores. We used an pixel count value averaged between the two halves of the cores for each bin to get a more robust measure of the luminophore distribution through the sediment. We smoothed the average profiles by species and sediment treatment using local polynomials [R package KernSmooth (Wand, 2019)] to better show general patterns in mixing.

When examining the luminophore profiles we looked for subsurface peaks to differentiate between local (biodiffusion) and non-local (advective transport) mixing. Local mixing causes the number of luminophores to decrease exponentially with depth, whereas non-local mixing causes a peak in luminophores at depth. We determined whether there was non-local mixing by examining the depth of maximum values of luminophore pixel counts over a moving depth window. We were interested in how mixing activity might differ close to the surface and at depth, so we isolated the peaks in the diffusion layer (0–4 or 0–6 cm depending on species/sediment treatment) and in the lower part of the core. We compared the depth and intensity of peaks between species and sediment treatments.

Process Measurements: Burying Depth and Survival

To measure animal burrowing depth for the bivalves, we recorded the lowest point of the bivalve’s position. For the polychaetes, we separated the cores containing polychaetes into 0 to 2 cm, 2 to 5 cm, and 5 to 12 cm layers and sieved each layer and recorded which contained animals. We considered all unrecovered animals as dead.

Statistical Analysis

We wished to determine the effect of sediment bulk density on the behavioral response of our species. Because we wished to know whether the response to bulk density would be different depending on grain size, we evaluated models with only main effects of the sediment treatment (Response ∼ Species + grain size + bulk density) and models that included an interaction between sediment grain size and bulk density (Response ∼ Species + grain size × bulk density). Because we wanted to compare the direction of the response between species, we also tested for an interaction between the species and sediment treatments (Response ∼ Species × grain size × bulk density). We selected the best models based on Akaike Information Criterion (AIC). We mainly used ANOVAs and logistic regression to evaluate the animal behavioral responses to sediment bulk density. See Supplementary Table 2 for a summary of the logistic regression model and Supplementary Table 3 for a summary of the best ANOVA models.

The model fit for all linear regressions and ANOVAs were evaluated using residual plots and QQ-plots. The logarithm of the dependent variable (i.e., depth and luminophore pixel count) was used when appropriate. For all models, the species treatments included separate categories for the two bivalve size classes. After performing ANOVAs, we used post hoc Tukey tests to compare relevant sediment treatment effects, as well as differences in responses between size classes for C. edule and S. plana.

We modeled the initial burrowing time of the bivalves using logistic regression. The response variable was the presence of animals at the surface in each unit at each time step. Our model was of the form of: Response ∼ Species + grain size + bulk density + hour, with interactions tested as detailed above. We used Wald tests to compare sediment treatment effects and investigate the effect of bivalve size on burrowing speed. We also used logistic regression to examine whether sediment treatment had an effect on polychaete survival and burrowing time.

The oxygen data were analyzed in R [R package turbo (K. Soetaert and Provoost, 2017)]. We corrected for the respiration of the microbial community by subtracting the mean control value for each sediment treatment from the values for cores containing animals. We evaluated whether oxygen consumption significantly varied among species, bulk density, and grain size using a three-way ANOVA.

We performed a three-way ANOVA to test whether there were differences in burial depth among bivalve species and sediment bulk density, and grain size. To satisfy the condition of normality, we excluded zeros (non-buried S. plana). To compare polychaete depth strata, we performed a multinomial logistic regression using functions from R package nnet (W. Venables and Ripley, 2002).

Luminophore dispersal was calculated as an indicator for bioturbation intensity by summing the number of luminophore pixels counted in each 0.5 cm depth bin for each luminophore profile. The values of the luminophore dispersal increased when more luminophores were transported into the sediment. The final dispersal values were corrected for the average dispersal values of the controls. This correction allowed us to account for the greater permeability of the sandy sediments. We performed a three-way ANOVA to evaluate whether there were differences in luminophore dispersal among species, sediment bulk density, and grain size.

We also compared the depth and luminophore pixel counts of non-local mixing peaks for species that showed non-local mixing in the luminophore profiles: L. balthica, and H. diversicolor, and both size classes of S. plana. The peaks at depth for L. balthica were not as prominent in the profiles as for the other two species, however, there is evidence from literature that L. balthica induces non-local mixing (Morys et al., 2017). We evaluated whether there were differences in the depth and the value of the peak luminophore pixel counts between 4.5 and 10.5 cm among species and sediment treatments using a three-way ANOVA.

Results

Survival Rate and Growth

Survival was very high for the bivalves (99 ± 7%). We recovered the polychaetes at a lower rate (71 ± 33% of H. diversicolor and 20 ± 31% for A. marina), but it’s unclear to what extent this was due to mortality or due to the polychaetes escaping the cores. If we assume that the treatment effect was due to mortality, then A. marina survived at a higher rate in the soft sediments vs. the compact sediments (30 ± 39% vs. 10 ± 13% survival), Wald’s X² (1) = 7, p = 0.008. There was no sediment treatment effect on H. diversicolor survival.

Burrowing Rates of Bivalves

The burrowing rate was significantly slower for the compact, high bulk density treatments than for the soft, low bulk density treatments, Wald’s X² (1) p < 0.001 (Figure 3 and Supplementary Table 2). Burrowing was around 200 times slower in the compact-sandy treatment and 80 times slower in the compact-muddy treatment than the corresponding soft treatments.

The small S. plana had significantly faster burrowing rates than the larger S. plana, Wald’s X² (1) = 11.7, p = 0.001, but the effect was sediment dependent. The small S. plana burrowed 2 times faster than the large ones in the soft-sandy treatment and they burrowed 5 times faster than the large ones in the muddy-compact treatment, Wald’s X² (1) = 9.9, p = 0.002. The small C. edule burrowed 0.7 times faster than the large C. edule in the soft-sandy treatment, Wald’s X² (1) = 5.7, p = 0.017. Size had no effect on burrowing rate in the soft-muddy treatment (fastest burrowing rate overall) and compact-sandy treatment (slowest burrowing rate overall) for both species. And while all the small S. plana successfully burrowed, 60% of the large S. plana in the compact-sandy treatment and 30% of the large S. plana in the compact-muddy treatment failed to burrow at all.

Burrowing Rates of Polychaetes

The effect of bulk density on the burrowing of polychaetes followed a similar pattern to the bivalves, though at a different scale (i.e., note difference in the x-axis in Figure 3 - hours vs. Figure 4 - seconds). The polychaete burrowing rate significantly increased with the softness of the treatments (Figure 4), Wald’s X² (4) = 357.6, p < 0.001. The biggest difference in burrowing rates was between the two sediment treatments with the highest bulk densities: burrowing rate was 150 times slower in the 1.35 g cm^–3 treatment compared to the 1.31 g cm^–3 treatment, Wald’s X² (1) = 18.4, p < 0.001, suggesting that there might be an absolute threshold that prevents burrowing. Burrowing rates did not significantly differ between 1.26 and 1.22 g cm^–3 treatments, Wald’s X² (1) = 0.31, p = 0.58. H. diversicolor burrowed four times slower than A. marina, Wald’s X² (1) = 5.7, p = 0.017.

Oxygen Consumption as an Indicator of Metabolic Activity

Overall, oxygen consumption significantly differed between species, F (3, 53) = 32.480, p < 0.001, but not sediment treatments (Figure 5). The mean oxygen consumption in the control muddy sediments was 1.6 times greater than the oxygen consumption in the control sandy sediments. That is, the average muddy sediment oxygen consumption was 20.7 + 2.6 SE mmol O₂ m^–2 d^–1 vs. an average oxygen consumption of 13.3 + 2.3 mmol O₂ m^–2 d^–1 for sandy sediment, which might be explained by the greater amount of organic material available for microbial consumption in muddy sediment. There was no significant difference between the oxygen consumption of the large and small S. plana.

All the animals tended to be most active in the sediments that were closest to their natural habitat (see Table 1 for sediment preferences). C. edule and L. balthica tended to have the greatest oxygen consumption in the soft-sandy sediments whereas the H. diversicolor and S. plana tended to have the greatest oxygen consumption in the compact-muddy sediment (Figure 5). While the sediment treatment effect on the species’ oxygen consumption was not significant, the higher oxygen consumption in these sediments could indicate greater activity.

Macrofauna Burrowing Depth

Bivalve burrowing depth was significantly different between sediment bulk densities, grain size, and species, with interactions between the sediment characteristics and species (Figure 6 and Supplementary Tables 1, 3). A post hoc Tukey test showed that the sediment treatments affected the burrowing depth of two out of the five bivalve species: the L. balthica and small C. edule. The L. balthica burrowed deeper in the soft sediments and the small C. edule burrowed deeper in the muddy sediments (Figure 6). The lack of a visible treatment effect on depth for S. plana could be due to the short length (12 cm) of the cores. Many of the S. plana were found near the bottom of the cores and perhaps would have burrowed deeper had they been given more space. While there was a statistically significant difference in the burrowing depth for the two size classes of C. edule, this was not the case for S. plana (Tukey post hoc p = 0.0001 and 0.16, respectively). The best multinomial model of polychaete depth distribution included sediment treatment. Including species did not improve the model, X² (2) = 2.76, p = 0.25. There is small evidence that polychaetes burrowed deeper in compact muddy treatment than the soft muddy treatment (p = 0.057).

Sediment Mixing as Measured by Luminophores

Overall, the benthic macrofauna induced significantly greater luminophore dispersal in the soft sediments than in the compact sediments [Figure 7; F (1,154) = 32.7, p < 0.001]. The greatest luminophore dispersal occurred in the soft-muddy sediment (average = 4.02 × 10⁵ luminophore pixels ± 0.78 × 10⁵ luminophore pixels). The luminophore dispersal was greater in soft-sandy sediment (average 3.68 × 10⁵ luminophore pixels ± 1.75 × 10⁵ luminophore pixels) than in the compact sandy-sediment (average 2.97 × 10⁵ luminophore pixels ± 1.10 × 10⁵ luminophore pixels).

Species induced different amounts of sediment mixing, F (6,154) = 9.4, p < 0.001 (Figure 7). Both size classes of the S. plana and the polychaetes mixed greater amounts of luminophores into the sediment than the L. balthica and C. edule (Figure 7). There was no significant difference between the luminophore dispersal of either C. edule or S. plana size classes, though the small S. plana tended to induce greater sediment mixing than large S. plana (Figures 7, 8).

We further compared the bioturbation mode of the macrofauna by examining the shapes of the luminophore profiles, especially by looking at the presence and depth of non-local mixing (Figures 8, 9). C. edule performed local mixing in the top three centimeters, with the greatest activity occurring in the low bulk density sediment which is visible by the greater quantity of luminophores mixed into the low bulk density sediments than the high bulk density sediments at the same depth (Figure 8). The luminophore profiles for S. plana, L. balthica, and H. diversicolor, showed non-local mixing peaks at depth (Figures 8, 9 white depth interval). The L. balthica and H. diversicolor peaks seemed most prominent in the soft sediments, and the S. plana appeared to produce greater peaks in the soft-sandy sediment (Figure 8).

For the species that showed the clearest non-local mixing (S. plana and H. diversicolor), the maximum luminophore pixel counts at depth were greater in soft sediment than in compact sediment, F (1, 64) = 8.7, p = 0.004 (Supplementary Table 3). Across all species that exhibited non-local mixing, the peaks were greatest in the soft-sandy sediment (average = 0.57 × 10⁵ luminophore pixels ± 0.31 × 10⁵ luminophore pixels).

For all species, greater mixing occurred in the soft treatments vs. the compact treatments. The polychaetes transported more luminophore in the muddy treatments than the sandy ones. Of all the species, S. plana mixed the greatest amount of luminophores into the sediment and was the only species that tended to transport more luminophores in the soft-sandy treatment than the soft-muddy treatment. H. diversicolor and S. plana had greater non-local mixing peak values at depth in the soft treatments than the compact treatments.

Discussion

The aim of our study was to investigate how bulk density affects benthic macrofauna behavior, as a step toward understanding how biogeomorphological and biogeochemical processes on intertidal flats might change under different bulk density conditions. We found clear effects of bulk density on benthic macrofauna burrowing behavior and bioturbation activity, which are summarized in Figure 10. In line with our first hypothesis that the mobility of the benthic macrofauna would be lower in compact (i.e., high bulk density) vs. soft (i.e., low bulk density) sediment, we found that in compact (vs. soft) sediment all animals burrowed slower, all animals transported fewer luminophores, and L. balthica burrowed shallower. Our second hypothesis was that bulk density would have similar effects on the behavior of benthic macrofauna in both sandy and muddy sediment. We did not detect any significant interactions between sediment grain size and sediment bulk density, indicating that animal responses had similar directions in both muddy and sandy sediments. We also did not detect a sediment treatment effect on respiration rate, however, the animals tended to have the greatest respiration in the sediment that were closer to their natural habitat. In line with our third hypothesis that larger bivalves would be more sensitive to differences in bulk density than smaller bivalves of the same species, we found that smaller C. edule and S. plana burrowed significantly faster than the larger ones in the compact muddy treatment (S. plana) and soft sandy treatment (C. edule). In addition, the small S. plana tended to transport more luminophores than the larger ones, but this appeared to be independent of sediment treatment.

Importance of Sediment Grain Size vs. Bulk Density Effects on Burrowing

We found it surprising that the burrowing rates of bivalves were significantly faster in the soft treatments compared to the compact treatments, rather than being largely driven by sediment grain size. Indeed, we expected the burrowing rates to be driven by penetration resistance and the muddy treatments had much lower penetration resistance than the sandy treatments (Figure 2 and Table 2). The observed burrowing patterns might be explained through the difference in sediment cohesiveness between the sandy and muddy treatments. Cohesiveness is largely governed by clay content (Joensuu et al., 2018). At high water contents, muddy sediment becomes akin to a viscous liquid and easy to entrain, whereas at low water contents, muddy sediments have much greater cohesion and a higher erosion threshold (Grabowski et al., 2011). In our study, although the compact-muddy sediment was more penetrable and had a lower absolute bulk density than the soft-sandy sediment (see Table 2), it was noticeably more cohesive and stickier than any of the other sediment treatments. The biomechanics of burrowing are different depending on the sediment type (Crane and Merz, 2017): in cohesive mud, animals burrow through crack propagation, whereas in coarse sand they may burrow through local fluidization or excavation (Dorgan, 2015). The cohesiveness of the compact-muddy sediment could have presented an obstacle to the bivalves’ burrowing of similar magnitude to the high penetration resistance of the compact-sandy sediment, and probably affected the biomechanics of the burrowing animals.

Though we found that within our experiment behavioral differences between the treatments were mainly driven by sediment bulk density rather than grain size, in nature the sediment grain size provides important constraints for species habitat. Many studies have described species assemblages to vary along a sediment grain size gradient (e.g., Ysebaert et al., 2002; Thrush et al., 2003; Compton et al., 2013; Pratt et al., 2014). However, most of these studies are correlative and the mechanics that underpin the habitat-animal associations prove to be elusive (Snelgrove and Butman, 1994). Other factors, like hydrodynamics, may be equally important. For example, the sandy areas in the Western Scheldt have a high degree of hydrodynamic stress and have impoverished benthos communities compared with the species-rich sandy areas in the Eastern Scheldt which have low hydrodynamic stress (Cozzoli et al., 2013). Bulk-density effects are generally not included in these kind of field studies. Our study highlights that including bulk density measurements may add an extra level of understanding to benthic macrofauna distribution and especially macrofauna activity in terms of sediment mixing.

High burrowing ability is essential for bivalves to survive in unstable sediments (Alexander et al., 1993; Takeuchi et al., 2015). All animals in our study burrowed faster in the soft sediments, and the small C. edule and S. plana burrowed faster than the larger adult ones in certain treatments. However, soft sediment is more easily eroded than compact sediment (Grabowski et al., 2011). In most bivalve species, like C. edule and L. balthica, juveniles live closer to the surface than adults and may hence be more easily dislodged by erosion during storms (Tallqvist, 2001; St-Onge and Miron, 2007). Yet, smaller and younger bivalves may compensate for their shallower living depth during erosion events by being able to burrow faster than adults in high bulk density conditions. It would be interesting to further investigate whether there is a tradeoff between sediment compaction effects on the erodibility and burrowing rate, and how this might impact the overwintering survival of juvenile versus adult bivalves.

Sediment Bulk Density Effects on Benthic Macrofauna Survival and Predation

Extreme bulk density sediments may present difficult living conditions for benthic macrofauna and affect their survival due to physiological constraints. At very low bulk densities, sediments might become so soft that animals have to expend a great amount of energy to keep their position in the sediment or unclog their feeding apparatus of small mud particles (Lohrer et al., 2006; Mestdagh et al., 2018). High bulk densities would present different challenges. For example, our high bulk density treatments most likely inhibited A. marina’s ability to ventilate their burrows which is energetically costly in low sediment permeability (Meysman et al., 2005), thus greatly reducing their survival.

Other ecological mechanisms, like predation risk and growth efficiency, might be affected by sediment bulk density as well. At shallower depths, the feeding area of deposit-feeding bivalves is increased as the siphon can be extended onto a larger surface area (Zwarts et al., 1994), whereas bivalves respond to predator presence by burrowing deeper (Griffiths and Richardson, 2006; Flynn and Smee, 2010). In high bulk density sediment, L. balthica may have increased energy expenditure during burrowing or feeding, which would reduce L. balthica growth efficiency. In addition, L. balthica would also stay closer to the surface in high bulk density sediment, thereby increasing its vulnerability to predation. Similarly, when H. diversicolor feeds at the surface, it is vulnerable to predation and its escape depends on speed (Ens et al., 1997). Considering that H. diversicolor burrowed significantly faster in softer low bulk density sediments, we can conclude that H. diversicolor and other surface deposit feeding polychaetes might also be more vulnerable to predators in high bulk density sediments because of decreased possibility of escape. Thus benthic animals may have lower survival in high bulk density conditions due to increased risk of predation.

Implications for Animal-Sediment Interactions and Ecosystem-Scale Impacts

Benthic macrofauna change their behavior under different bulk densities which can have consequences for the biogeomorphology and biogeochemistry of tidal flats. Under higher bulk densities, a reduction in the depth at which infauna burrow, which we observed in particular for L. balthica (Figure 6), may lead to a shallower apparent redox potential discontinuity (Gerwing et al., 2017), which could decrease the depth of the biologically active zone depth (Sturdivant and Shimizu, 2017). Furthermore, a decrease in bioturbation and especially non-local mixing activity in high bulk density sediment, which we observed for important gallery diffusor H. diversicolor (Figure 8), could lead to reduced sediment permeability and oxygen penetration of the sediment (Aller and Aller, 1998; Michaud et al., 2006) as well as a build-up of a mud layer due to decreased resuspension or incorporation of mud into the sediment matrix (Montserrat et al., 2008; McCartain et al., 2017). In addition, a reduction in non-local mixing of H. diversicolor would decrease microbial processing of organic material which would reduce nutrient release into the porewater (Mermillod-Blondin et al., 2004). Changes in bivalve behavior may also affect tidal flat biogeochemistry, for example, reduced siphon movement and pumping rate due to more compact high bulk density sediment would increase the time between bouts of oxygenation of the sediment which might lead to short-term anoxic conditions and decreased denitrification (Volkenborn et al., 2012). Thus higher sediment bulk densities may have negative consequences for the ecosystem function of benthic macrofauna. Indeed, a reduction in their burrowing depth and bioturbation activity could lead to a shallower and less well-oxygenated surface sediment layer which would impact microorganisms and eventually nutrient cycling. Low bulk density sediments may have the opposite effect and stimulate nutrient cycling due to increased bioturbation.

Because sediment bulk density tends to increase with sediment erosion and decrease with sediment deposition (Dyer et al., 2000), we speculate that the macrofauna in an eroding tidal flat are typically less mobile than macrofauna in a depositing tidal flat of a similar sediment grain size. Animal-sediment interactions between bioturbation and bulk density are likely to create positive feedback loops for tidal flat biogeomorphology. Bioturbation destabilizes sediment (Widdows et al., 2000) which decreases sediment bulk density and, as we found in our study, a soft sediment encourages animal movement. These interactions could create a positive feedback loop between low sediment bulk density conditions and elevated benthic macrofauna movement. The opposite feedback loop would occur under high bulk density sediment where sediment conditions discourage animal movement which in turn may lead to further sediment compaction. However, non-cohesive soft sediments are generally more vulnerable to erosion than compact sediments (Grabowski et al., 2011), and greater animal activity may cause these sediments to be even more easily eroded. Further investigation on animal-sediment interactions, particularly on whether animal activity under different bulk densities affects tidal flat biogeochemistry and generates positive biogeomorphology feedback loops, would be necessary to tease out the importance of sediment bulk densities for tidal flat functioning and evolution.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

LW, TB, TY, and CM contributed to the study’s conception. LW and NS contributed to the study’s data acquisition and data analysis. LW wrote the first draft. All authors contributed to revising the manuscript and approved the final manuscript.

Funding

This study was funded by the Buitendijkse project of Rijkswaterstaat (Netherlands Ministry of Infrastructure and Water Management) and supported by the Netherlands Organisation for Scientific Research (NWO) via the project “EMERGO – Ecomorphological functioning and management of tidal flats” (850.13.020, 850.13.022, and 850.13.023).

Conflict of Interest

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.

Conflict of Interest

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 would like to thank Werna Werna for help in carrying out the experiment. We would also like to thank Jeroen van Dalen, Lennart van Ijzerloo, Daniel Blok, Anton Tramper, and Justin Tiano for technical assistance.

Supplementary Material

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

References