The combined influence of body size and density on cohesive sediment resuspension by bioturbators

We propose an empirical framework to scale the effects of bioturbation on sediment resuspension to population bioturbation activity, approximated as population metabolic rate. Individual metabolic rates have been estimated as functions of body size and extrapolated to population level. We used experimental flumes to test this approach across different types of marine, soft-sediment bioturbators. We observed that a large part of the variance in biota-mediated sediment resuspension can be explained by a positive relationship with population metabolic rate. Other mechanisms can strongly influence the outcome, such as bioturbation of deep sediment strata, biotic interactions with hydrodynamic stress and overlapping areas of influence must be further investigated. By relating the biota-mediated changes in resuspended sediment to metabolism, we can place our observations within the broader context of the metabolic theory of ecology and to formulate general expectations about changes in biota-mediated sediment resuspension in response to changes in population structure and climate change.

Formulating a general quantitative framework to measure bioturbation processes is made difficult by the high diversity in macrozoobenthic species in terms of distribution, characteristics and functional behaviours [24][25][26][27] . In addition to characteristics of benthic populations such as overall biomass and density, previous approaches to this topic have used functional trait classification for bioturbators that requires taxonomic and life-history expertise, which are not always available (e.g. 5,12,13,27,28 ). However, other empirical studies have pointed out that simple knowledge of an organism's biovolume may allow for relatively simple estimates of the reworking intensity of surface sediment across a range of functional diversity 29 . To develop a general approach to scale the effects of bioturbation on sediment erodability across a broad range of functional groups, we tested the effect of a fundamental ecological descriptor, the metabolic rate.
Theoretical framework. The positive scaling of bioturbation rates with metabolism can be assumed via the positive scaling of respiration, feeding and moving rates. An organism's metabolic rate I (W) is 30,31 expected to scale with individual size M (g), with an allometric exponent of ~0. 75.
where a is a normalization constant accounting for variability in metabolic rates independent from the effect of size (e.g. thermoregulatory strategy 32 ). The positive size scaling of metabolic rates expressed from Equation 1 is considered to be one of the most fundamental patterns in ecology [33][34][35][36][37][38] , however, the actual values of the measured scaling exponents can exhibit large environmental 39,40 , taxonomic (e.g. [41][42][43] ) and phenotypic (e.g. 44 ) variability so that the "universality" of the ~0.75 scaling exponent is subject to ongoing debate. Supporting a metabolic approach, Gilbert et al 29 . showed that individual biovolume is a good predictor for sediment reworking intensity. Significant and positive power laws have been commonly observed between the size of macrozoobenthic organisms, their standard metabolic rate [45][46][47] and other activities affecting stability, such as sediment ingestion/egestion (e.g. 48 for A. marina 49 , for benthic detritivores 50 , for C. edule), burial and transport (e.g. bioturbation potential 12 , burying depth of A. marina, 48 ; burying depth of L. balthica 51 ). Ideally, if macrozoobenthos rework sediment proportionally to their individual energetic requirements, the amount of sediment loosened by a bioturbator and made available for resuspension (R, g) should change with the bioturbator size proportionally to the metabolic rate I:R I (2) The relationship between body size and density is an essential link between individual and population traits, such as population spatial density 34,37,[52][53][54] . The overall energetic requirements of for a population of N (N of Ind. m −2 ) bioturbators of M size, I TOT (W m −2 ), may be estimated as the product of the individual metabolic rates and the individual density 55 : TOT It follows that, at the variation of I TOT , the overall amount of resuspended sediment R TOT (g m −2 ) may be approximated from the linear model

TOT T OT
where c is a normalization constant accounting for the amount of sediment that is resuspended due to physical forces only (i.e. in absence of bioturbation) and the coefficient d accounts for the relation between variation in I TOT and variation in R TOT . The coefficient c is expected to be strongly dependent on the physical conditions under which the sediment resuspension happens 56,57 . The coefficient d may vary according to the modalities that different functional groups of bioturbators have in reworking sediment 58 .

Experimental test.
Using artificially created single-species or similar functional group assemblages allows us to isolate particular aspect of bioturbation, and thereby provide much needed mechanistic insight 13 . In this study, we test the hypothesis that the amount of resuspended sediment due to bioturbation action is proportional to the overall activity of the bioturbator population, expressed as a linear combination of individual metabolic rates and density of individuals. Thus, we used an experimental set-up that excluded variation in physical conditions (e.g sediment granulometry, cohesiveness, compaction, shear stress 57 ,) or metabolic and behavioural changes in response to environmental cues (e.g. sediment composition, 9 ; acidification, 39 ; temperature, 59 ; food availability 60 ; shear stress 61 ,). This simplification was made to test if bioturbation fits within the ecological framework of size-dependent energetic theories (e.g. 33 ). We tested our hypotheses using an empirical dataset. This dataset uses the mass of resuspended sediment (MRS, g m −2 ) as a measure of bioturbation effects on R TOT . The MRS in the water is coupled with the mass of bottom sediment by a dynamic balance between deposition and erosion. Increasing bottom shear stress has the effect to increase the sediment erosion and decrease the sediment deposition, thus increasing the MRS. By loosening the sediment with their activities, benthic bioturbators increase the amount of sediment available for resuspension at a certain current velocity, shifting the balance point between erosion and deposition to a higher MRS 5 . Analogously to previous studies (e.g. 3,5 ) we did not measure sediment deposition and we only consider the effect of bioturbations on the equilibrium MRS (deposition rate = erosion rate) reached at a given level of bed shear stress from water motion.
The dataset encompasses a range of functional diversity (from shallow to deep bioturbators), individual densities (13 to 382 Ind. m −2 ) and body sizes (10 to 1136 mg Ash Free Dry Weight, AFDW). Individual metabolic rates were estimated according to the empirical model of Brey 62 . Three functionally different types of bioturbation activity were accounted for in the analysis 58 : (i) shallow-burrowing bivalves that make very shallow perturbations in the sediment by crawling on the surface, shaking valves and producing pelleted pseudo-faeces, represented by C. edule; (ii) Intermediate Burrowing Bivalves (IBBs) that disrupt the sediment surface by inhaling the sediment through their siphons and depositing pseudo-faeces, represented by Abra alba, Scrobicularia plana, L. balthica and Ruditapes philippinarum; and (iii) deep-burrowing Polychaeta that swallow surface sediment through a feeding funnel and expel it in the form of pseudo-faeces, forming characteristic feeding pits and pseudo-faeces casts, represented by A. marina. By using three contrasting functional groups with each group containing individuals that greatly vary in body size, we were able to measure the changes in sediment resuspension due to different modes of mechanical destabilisation across a large range of body sizes and densities.

Results
The estimated individual metabolic rate of experimental organisms varied from 0.07 mW [+/−0.02 95  cumulative effects on errors on the animal measurements, on the conversion between length or wet weight to AFDW and of the conversion between AFDW and metabolic rate. On average, 32.25 g m −2 [+/−1.73 95 CI] of sediment was resuspended in the control runs without bioturbators at a bed shear stress of 0.18 Pa. When the same shear stress was applied after 48 h of active bioturbation in the flumes, we observed a major difference in sediment resuspension: R TOT increased up to almost 5 times in the bioturbated runs than in the non-bioturbated controls. The highest R TOT (145.29 g m −2 [+/−100.92 95 CI]) was observed during activity by the highest density of the largest A. marina. Only in one case (the lowest density of the smallest A. marina) did we observe a lower amount of resuspended sediment than in the defanauted control ( Fig. 1).
Spearman's rank correlation between overall population metabolism and R TOT is 0.82 (p-value < 0.001), consistently exceeding the corresponding correlations with individuals' density only (0.37, p-value < 0.05). A large part of the observed variation in R TOT (R 2 = 0.67) could be described as a linear function of the overall population metabolism: The estimated intercept (35.11 g m −2 [+/−6.60 95 CI]) is consistent with the resuspended sediment observed in the defanauted control. Stepwise-model simplification showed that differences across functional groups were not statistically significant (p-value > 0.1, Table 2).The trend is still significant (p-value < 0.05) and positive (scaling coefficient 0.17 [+/−0.15]) when the two observations with higher I TOT are removed. (Appendix C, Table C3). We did not detect linear correlation between the amount of sediment resuspended per unit of metabolic power R BIO = (R TOT − R CONTROL )/I TOT (g mW −1 ) and both bioturbator density and individual size (

Discussion
In this paper, we derived relationships across species to describe the potential effect of bioturbators on the amount of sediment made available for resuspension. We hypothesised that the increase per area of resuspended sediment due to mechanical destabilisation by bioturbation is proportional to the overall activity by the bioturbator population, so that it can be quantified as a positive function of individual metabolism and density of individuals.
We were able to test these hypotheses by measuring the effect of different functional groups of bioturbators on sediment erodability. Confirming our expectations, the positive scaling with population metabolic rate was able to explain a large portion of the variance (67%) in loosened and resuspended sediment due to bioturbation. The scaling relationship was not significantly different among the investigated groups, supporting the idea of common energetic constraints acting across a large range of functional diversity (Fig. 1, Table 2). Equation 5 provides an empiric description of bioturbation effects on sediment resuspension. However, the expected link may not be mechanistically transferable from the individual to the population level purely on the basis of body size-density relationships. As an example, individuals may have overlapping areas of influence, that vary with organisms' density 3,5 and size 63 . We observed a constant amount of resuspended sediment per unit of metabolic power (R BIO ) across density and size gradients ( Table 3), indicating that negative interference was not apparent in the data. However, it must be considered that such interference may occur at the increase of both the individuals' size and density. As an example, at L. balthica densities (500, 1000, 1500 Ind m −2 ) higher than those we tested in our experiment (maximal tested density 382 Ind m −2 ), the per capita amount of bioturbated sediment decreases with increasing density 3 . It follows that a limiting function should be used to apply Equation 5 to very high individuals densities (e.g. >500 Ind m −2 for L. balthica 3,5 ,). As an example, Willows et al. 3 used an asymptotic function to model the relationship between L. balthica density and resuspended sediment, and van Prooijen et al. 5 used an exponential function to model the probability of overlap in influence area. However, such high densities of bioturbators are not very likely to be realized in nature, e.g. L. balthica has an individual density lower than the maximal one we tested, 382 Ind m −2 , in the 97% of the records collected between 2005 and 2011 in the Westerschelde and Oosterschelde, 64 . Thus, the use of a limiting function for individuals' density may eventually be neglected for broad field applications of Equation 5.
Our measures, similarly to those of previous studies 3,5 focus on sediment resuspension and neglect effects on deposition that may be relevant at high individuals density or at the decrease of the bed shear stress. Because bioturbators can increase the sediment surface roughness with their physical shape (i.e. protruding shells of C. edule) or with the structures they produce (i.e. feeding funnels of Intermediate Burrowing Bivalves (IBBs), pseudo-faeces casts of A. marina), high densities of bioturbators could dampen the near-bottom water flux and reduce the bed shear stress, thereby trapping the sediment in saltation before resuspension. The reworking of the sediment by the animals could also change the characteristics of the sediment: larger flocs can be broken, or excreted sediments can be pelletized and compacted. The finest part of the sediment, exposed to the buoyant action of water from sediment mixing, may be eroded at first, increasing the average granulometry and decreasing the cohesiveness of the bottom sediment 65 . It is possible that these processes may change the sediment response to the water action, inducing changes in the equilibrium MRS. Finally, it is also possible that individuals could   decrease their metabolic/activity rates with increasing density (e.g. 66,67 ), generating negative covariance between I and N. The aforementioned sources of negative correlations between bioturbators size, density and metabolic rate may generate non linearites in the relationship between R TOT and I TOT . At the scale of our analysis, these non-linearities are intrinsically included in Equation 5 and may result in a smaller scaling coefficient c. Detailed measurements are required to a deeper understanding of these interactions. It must be considered that our model has also a substantial prediction error and might be affected from the two observations with higher overall metabolic rates (largest sizes and densities of A. marina and C. edule, Fig. 1). This indicates that the experimental design could be improved by increasing replications and distributing experimental effort more evenly across the range of overall population metabolism. However, excluding the two observations with higher leverage the scaling trend is still positive and significant (p-value < 0.05 Appenidx C, Table C3). The logistic efforts necessary for empirical testing (in our case, each single run took ca. 1 week of preparation followed by 1 week of experiments) it is the main reasons that our dataset is limited. While the dataset we collected may be regarded as not being "optimal", it is one of the most complete experimental datasets (to our knowledge) on biota-mediated sediment resuspension that has been measured according to metabolic and density gradients across different species.
Comparison with previous approaches. Positive values for size and density scaling have previously been proposed as descriptors of bioturbation potential. For example, the bioturbation potential index (BP c 12,24,27 ,) for marine ecosystems approximates the effects of benthic bioturbators as a linear combination of abundance, average size and behavioural traits. Applications of BP c showed that benthic community structure can be actually used to predict the process of bioturbation in real ecosystems 13 . However, this approach has some limitations. First, it requires extensive knowledge of species life histories and high taxonomic expertise. In particular, the BP c index relies on detailed functional classification of organism traits associated with sediment mixing and motility 27 . The paucity of such information for the majority of marine species is a source of concern, as this information is needed to project potential changes in BP c under future policy or environmental scenarios 27 . Second, combining individual size, density and functional traits is potentially redundant in the sense that different functional groups and motility classes may already contain relationships combining body size and size/density ratios 68 . Third, the BP c index, being based on a linear combination of abundance, biomass and behavioural traits, represents changes in bioturbation potential as a function of changes in species composition (which is difficult to predict in any detail) and in the biomass/abundance ratio on which typical body size is calculated 27 . Our experiment showed that different functional behaviours and motility classes did not have significant influence the amount of sediment resupended due to bioturbation. The largest part of observed variance in resuspended sediment could be explained simply in terms of bioturbators' population metabolic rate. Hence, a coarser and less taxonomically demanding classification of bioturbators than used for the BP c index would still able to predict the effect of bioturbation processes on sediment resuspension. This is supported by other studies showing that indicators based on community size structure, rather than on species-specific characteristics, can be used for describing the ecosystem status for some functional aspects 69 . Scope for extrapolation. Sediment resuspension can originate from very diverse interactions and activities 4 , which may also depend on external conditions. While our measures allow for general quantifications of scaling effects across a range of bioturbator diversity, other issues must be investigated for applications to field contexts. Seawater temperature is a key regulator of metabolic rates in ectotherms, such as macrofaunal invertebrates 70 . Individual metabolic rates are expected to increase with higher temperature within the ranges of physiological thermal tolerance 32,59,[70][71][72][73][74] . From this perspective, investigations into the effects of temperature on bioturbation rates are relevant for estimating the potential effect of seasonal, latitudinal and global climate change on biota-mediated sediment resuspension. Increases in bioturbation activity with increasing temperature have been observed (e.g. [75][76][77], supporting the hypothesis that there is a metabolic dependence of large-scale bioturbation effects. Although our experiments were performed at a fixed temperature, the amount of resuspended sediment was related to the population metabolism than to density, biomass or biovolume (e.g. 29 ). Having such a metabolism-based relationship does enable speculative extrapolations on how predicted temperature increases may influence bioturbation effects.
The temperature dependence of metabolic rates can be described according to the Boltzmann-Arrhenius exponential model e −E/kT , where k is Boltzmann's constant, E is the mean activation energy of metabolic reactions and T is the Kelvin temperature 32,59,[70][71][72][73][74] . Assuming that all other mechanisms involved in sediment resuspension do not change with temperature variations, expectations about the effect of increasing temperature on R TOT via metabolic rates can be formulated including a Boltzmann-Arrhenius temperature scaling 59  increase in summer temperature of 3 °C (as it is expected to happen in the southern part of the North Sea before the end of this century 78 ,) should imply an increase in metabolic rates of 30% with respect to our reference temperature. Assuming that increases in individual metabolism will not be energetically balanced by a reduction in population density 40 , the positive scaling of metabolism with temperature would imply a consistent increase in bioturbation and biota-mediated sediment resuspension due to global warming at the end of this century.
Our estimations of changes in metabolic/bioturbation activity according to temperature are in the range of published empirical observations of variation in bioturbation activities [75][76][77]79 with respect to temperature change. For practical applications, population metabolic rates and their effect on sediment resuspension can be estimated with good approximation from survey or predicted data on benthic community composition by using empirical models (e.g. 62 ) or from theoretical expectations on size scaling of community metabolic rates (e.g. 33 ). Together with fluctuations in population structure and size density ratio 80,81 , resource availability and behavioural adaptations 82 , scaling relationships between temperature, metabolism and R TOT may contribute to explain and model the temporal 83 and spatial 84 variation observed in field bioturbation activity and contribution to sediment resuspension. Other parameters influencing aquatic invertebrates metabolism and energy allocation (e.g. depth, sex, age 44,62 ) may be included to adapt the metabolic rate according to the environmental conditions and physiological status of the bioturbating population to which Equation 5 is applied.
Beyond their effect on metabolic rates, the influence of local environmental conditions (e.g. sediment granulometry, current velocity, anthropogenic disturbance) should be assessed because this can affect bioturbator distribution (e.g. [85][86][87], as well as their interaction with the physical factors involved in sediment resuspension (e.g. shear stress 61 , temperature and acidification 39 , sediment composition 9 , and wave exposure 88 ). These factors can strongly influence the outcome of bio-mediated sediment resuspension. It is straightforward that different combinations of environment physical properties as sediment properties and bed shear stress will result in a different coefficient for physical erosion c in Equation 4 57 . Also the coefficient d relating population metabolism and sediment resuspension may vary as a function of environmental conditions. For example, it is possible that at values of bed shear stress higher than those tested in our experiment (i.e. allowing the water current to erode the deeper sediment strata and to completely smooth out all the roughness of the sediment surface generated by high densities of bioturbators) may result in steeper relationships being observed. However, such an increase in water energy may lead to physical factors and unpredictable erosion patterns (scouring) overcoming the importance of biological factors in determining sediment resuspension. Also, (i) intense shear stress is generally associated with non-cohesive sediment, where bioturbators have no or limited effect on sediment resuspension 9 and (ii) high densities and large sizes of bioturbators are generally associated with low shear stress in nature 64,86 , which makes the combination of strong shear stress, cohesive sediment, high bioturbator density and large bioturbator size unrealistic.
Finally, our measurements were focused on single species and sizes classes in order to emphasize scaling relationships. However, the effects of individual species on sediment resuspension in a mixed benthic community may be complex, depending on how interspecific interactions affect the activity of the involved species; these must also be accounted for in order to extrapolate mesocosm observations to field contexts 11,89 . For example, positive correlations have been observed between benthic diatom abundance and both C. edule 90 and A. marina 91 biomass. Benthic diatoms are well known sediment stabilisers, able to glue together sediment grains by producing extracellular polymeric substance 92 and to increase sediment resistance to erosion 93,94 . By disrupting and grazing the diatom film, benthic bioturbators may have a much higher relative impact on mudflat morphology than what we measured in our flumes because they are able to trigger the resuspension of sediment that is otherwise stabilised by diatoms 95 .

Conclusions
Empirical descriptions of the behaviour of organisms are needed for integrated modelling of bio-mediated physical processes 5,96 . With this study, we developed a general approach to scale the effects of bioturbation on sediment erodability across a broad range of functional groups. By following the general predictions from body size effects on the energy expenditure of individual activities 33,62,97 , and scaling up the energetic budget from the individual to population level 55 , we showed that the effect of bioturbators on cohesive sediment resuspension can be described simply in terms of bioturbators' population metabolic rate. While our quantitative estimation of increasing R TOT with community metabolism must be treated with caution due to the relatively limited extent of our dataset and because we did not directly test the effect of temperature or other influential factors, it can still be indicative of trends in ecosystem functioning. Being based on such a highly fundamental ecological descriptor as metabolic rate, our observations can be placed within the framework of general ecological allometric theories and in particular to the Metabolic Theory of Ecology 33 , allowing to formulate general expectations about present and future trend in biotic contribution to sediment resuspension based on expected variations in community composition (e.g. 80,81,98 ).

Material and Methods
Target organisms. The ecosystem engineers used for this experiment were the bivalves Cerastoderma edule, Limecola balthica, Abra alba, Scrobicularia plana and Ruditapes philippinarum, and the Polychaeta Arenicola marina. These organisms share a common habitat (mainly muddy intertidal flats), but they live in the sediment at different depths: from very shallow (C. edule, shells usually emerge from the sediment surface) to intermediate depths of 3 to 10 cm (L. balthica, A. alba, S. plana and R. philippinarum; grouped together as Intermediate Burrowing Bivalves or IBBs) to relatively deep depths below 10 cm (A. marina) 25,26 . Accordingly, their feeding modes vary from obligate suspension feeding (C. edule) to a mixture of suspension and deposit feeding (IBBs) to obligate deposit feeding (A. marina) 25 suspension feeders have been observed to have a lower size density ratio with respect to obligate deposit feeders, implying a different resource availability for the two guilds 68 .
The selected species are representative of three qualitatively different types of bioturbation activity 58 . C. edule reworking of sediment is mostly related to bio-deposition, vertical and horizontal movements and valve adduction 99 . L. balthica, A. alba, S. plana and R. philippinarum all burrow within ca. 5 cm from the sediment surface, and they can disrupt the sediment surface by inhaling the sediment with their siphons to graze on benthic diatoms 63 . A. marina swallow surface sediment through a feeding funnel and expel it in the form of pseudo-faeces, forming characteristic feeding pits and pseudo-faeces casts 6 . Given their similar function with respect to their sediment reworking modality, which sets them clearly apart from both A. marina and C. edule, we grouped these species into a single homogeneous group, i.e. IBBs. We also used this pooled approach to generate enough variation in body size and density for the intermediate burrowing, as this could not be realized at the species level in contrast to A. marina and C. edule. Not enough data were available to apply a robust statistical test to assess for inter-specific differences between size and density scaling within the IBB group. However, we did not detect any significant inter-specific differences when comparing the amount of suspended sediment observed from experiments using similar densities and sizes of the four species included in this group (see Appendix C, Tables C1 and  C2). In addition, our direct observations and those of many other authors collected over the years (e.g. 27,58,63,100,101 ) indicate that these organisms share common lifestyles, and modes of feeding and mobility. Finally, three of these species (A. alba, L. balthica and S. plana) belong to the same suborder (Tellinacea), two of them (A. alba and S. plana) to the same family (Semelidae), and they possess many physiological and structural similarities 102 .
The tested combinations of bioturbators's body sizes and densities were selected in a way to cover the natural range of each analysed species (e.g. 6,25,26,48 ) and according to the availability of experimental organisms (Table 1). Bioturbators' mass (mg AFDW) was estimated from the specimens length (mm, bivalves) or wet weight (mg, A. marina) according to the relationships provided from the NIOZ -Yerseke Monitor Taskforce. Bioturbators' individual metabolic rates were estimated according to the empirical model for acquatic macroinvertebrates respiration of Brey 62 using a trait classification for sessile intertidal satiate Anellida and Bivalvia Heterodonta operating at 18 °C and assuming an average energy density of 21.5 J mg −1103 . See Appendix A for more details about specimens' measurements and calculation of metabolic rates.
Considering the large total number of flume runs needed (32 different combinations of size, density and functional groups × 2 replicates × 3 runs for combination = 192 runs, Table 1) and the time-consuming character of each flume experiment, the animals were collected between May 2011 and May 2012 from the intertidal flats of the Oosterschelde and Westerschelde. At the time of collection, average daily water temperature was between 14 and 17 °C. After collection, the animals were always allowed to acclimate for 1 week in containers filled with sediment and aerated filtered marine water that was kept at 18 °C (water temperature in the Westerschelde during full summer). Considering the relatively limited difference in temperature between field and mesocosms, one week of acclimation (rather than the two weeks usually adopted in macrozoobenthos studies) should be sufficient to reduce the risk of temperature shocks that could severely affect bioturbator metabolic rates. However, it is still possible that small deviations in bioturbator basal metabolic rates due on partial acclimation may have induced some minor bias in our estimates. Experiments were performed directly after this week of acclimation. Each flume run always used homogeneously sized individuals of a single species. For each species, different densities of individuals were tested in separate runs. Experimental equipment. The annular flumes we used are a variation of the design described by Widdows et al. 61 (Appendix B, Figures B1-B3). The annular channel has a surface of 0.157 m −2 . In the majority of the cases, we used flumes with an overall height of 40 cm. A modified version with an overall height of 80 cm was used to have a higher sediment column and allow the largest sized A. marina to settle properly. To avoid abiotic variability in resuspended sediment due to different sediment characteristics 94 in the experiments, homogeneous, wet, muddy sediment (median grain size = 100 μm, silt content 12%) was put in a flume, mixed to a smooth mass and allowed to consolidate until creating a layer of ca. 10 cm height in the shorter flumes and of ca. 50 cm in the taller ones. Excess water in the sediment was drained through a pebbled bed placed at the bottom of the flume. After 48 hours, the flumes were filled with 31.4 L of filtered seawater (height of the water column 20 cm). To prevent damaging the sediment surface, a sheet of bubble plastic was placed on top of it before gently spraying water into the flume. A water current of 30 cm sec −1 was applied, corresponding to a bed shear stress of 0.18 Pa, which should be sufficient to resuspend the bioturbated sediment 61 . To apply the current, we used a smooth, adjustable rotating disk, which was driven by a microprocessor-controlled engine. An acoustic Doppler velocimetry probe was used to calibrate water velocity as a function of engine rotation speed. Water turbidity, as a proxy of resuspended sediment, was measured using an optical backscatter sensor (OBS 3+, Campbell scientific) facing the water perpendicularly to the current direction at 10 cm from the sediment surface, and converted into suspended sediment concentration (g L −1 ) based on calibration by gravitometric analysis (Appendix B, Table B1). To express sediment resuspension in spatial units, we converted the suspended sediment concentration to Mass of Resuspended Sediment (MRS, g m −2 ). Previous studies 5,61 have shown that in cohesive sediment, mostly supply-limited erosion occurs, i.e. after the water motion has started, the MRS reaches equilibrium due to limitation of erodible material 56 . In our experiments, the equilibrium MRS was usually reached after ca. 5 minutes of applying current. Experimental protocol. Every experiment (2 replicates) included a preliminary run, a control run without added animals and an experimental run with benthic animals. All runs lasted 20 min and were repeated at intervals of 48 h. The aim of the first, preliminary run was to further smoothen and homogenise the sediment surface of the flume bottom. As a consequence of the limited erosion that occurred during this control run, a uniform, less than 0.5 mm-thick layer of fine sediment was deposited on the sediment surface of each flume within a few hours from the end of the run. Pilot experiments conducted in flumes without fauna, involving several sequential daily runs, showed some small differences across flumes but no increase in sediment resuspension compared to the control run on the first day. This implies that all the sediment available for resuspension in the absence of bioturbation had been eroded during the first run and was suspended again during the subsequent runs. This observation allows us (i) to use the second run without fauna as an independent measure for the amount of sediment resuspended due only to shear stress and (ii) to subsequently use this measure as an internal control to quantify the bioturbation effect, while minimizing the small differences across flumes.
Immediately after the control run, animals were introduced into the flumes and evenly distributed over the sediment surface. The aim of the third, experimental run was to measure the change in sediment resuspension with respect to the abiotic control resulting from the action of the bioturbators during the 48 hours in which they remained in the flumes. The choice of a longer time interval (48 h) compared with the typical interval between erosion stress peaks (typically 12 or 24 h in a tidal system) was necessary to give to the animals a chance to properly settle in the new environment and recover from manipulation stress. The vast majority of them were buried within a few minutes after being placed in the flume. A. marina generally did not move from the initial settlement point and produced a single feeding pit with a pseudo-faeces cast for each individual. Some individuals of C. edule and IBBs crawled on and below the sediment surface, leaving evident tracks. Data analysis. We tested the hypothesis that the effect of bioturbators on sediment resuspension (approximated as equilibrium Mass of Resuspended Sediment, g m −2 at a current shear stress of 0.18 Pa) is proportional to the overall population metabolic activity (I TOT , mW m −2 , calculated as the product of the individual metabolic rate I, mW, and the population density N, N of Ind. m −2 ). To assess for differences in intercepts and scaling coefficients across the three groups of bioturbators, multivariate linear regression models were fitted with full interaction terms, including the dummy predictive variable "Functional Group". Selection of predictive variables and interaction terms was assessed by a bi-directional elimination stepwise procedure.
The correlation between individuals' size, density and amount of resuspended sediment per unit of metabolic power, R BIO = (R TOT − R CONTROL )/I TOT , g mW −1 , was investigated to asses for interferences between individuals at increasing individuals' density and size. Also in this case, selection of predictive variables and interaction terms was assessed by a bi-directional elimination stepwise procedure.
All analyses were performed within the free software environment R 3.3.2 104 .