. Density-dependent growth of bivalves dominating the intertidal zone of Banc

Abstract


Introduction
Density dependence is a key process contributing to long-term persistence and stability of populations (Murdoch, 1994).There is increasing evidence that temporal and spatial variability in environmental conditions, such as resource availability, energetic demands, predation risk or disease risk, affect the magnitude of density dependence (e.g.Wang et al. 2006, Finstad et al. 2009, Lok et al. 2013, Ford et al. 2016).For example, in a migratory bird species, density dependence in survival varied seasonally (Lok et al. 2013), whereas in fish, spatial variation in density-dependent growth and condition was driven by spatial variation in shelter availability (Finstad et al. 2009, Ford et al. 2016).
In addition to environmental conditions, species-specific traits can also influence the patterns and magnitude of density dependence.For example, group-living species, such as shoaling coral reef fish, may increase performance at higher densities as a result of better resource patch detection and reduced predation risk (Yeager et al. 2014).Consequently, understanding how environmental conditions and species-specific traits affect the magnitude of density dependence is crucial for making accurate predictions of population dynamics.
In soft-bottom-dwelling marine benthos, the extent to which populations are regulated by density-dependent processes may depend on feeding guild (Levinton 1972).In fact, Levinton proposed that infaunal deposit feeders, relying on rather constant and partially selfrenewing food supplies in the surrounding sediment, should be food-specialists occurring in densities at which they are limited by food availability.In contrast, infaunal suspension feeders, relying on the notoriously seasonally and locally variable phytoplankton in the overlying water, should be food-generalists, occurring in variable densities not closely regulated by food availability.The clarity of these predictions spawned a rich portfolio of descriptive and manipulative tests in soft-bottom marine benthic invertebrates, mostly conducted in temperate coastal systems (e.g.Ólafsson 1986, Peterson & Beal 1989, Kamermans et al. 1992).
As predicted by Levinton (1972), infaunal deposit feeders generally decrease shell growth rates at experimentally increased densities (e.g.Branch & Branch 1980, Ólafsson 1986, Kamermans et al. 1992).However, contrasting results were found for infaunal suspension feeders.Whereas some studies found no evidence for density dependence in growth of infaunal suspension feeders (e.g.Ólafsson 1986Ólafsson , Kamermans et al. 1992)), others did (e.g.Peterson & Beal 1989, Bijleveld et al. 2015).This indicates that the magnitude of density-dependent growth in soft-bottom marine benthos is not determined by feeding mode alone, but that environmental conditions also plays an important role.Yet, we are unaware of any experimental tests that have explored how species-specific feeding traits and environmental conditions (e.g., habitat characteristics, food availability) interact to affect the magnitude of density dependence in growth of soft-sediment marine benthic invertebrate populations.
When organisms colonize patches in proportion to the amount of resources available in each patch (i.e. according to an ideal free distribution, sensu Fretwell & Lucas 1970), this can obscure estimates of density dependence in demographic parameters (Wilson & Osenberg 2002, Shima & Osenberg 2003).This phenomenon is referred to as "cryptic density dependence" (Shima & Osenberg 2003).To identify density dependence, it is therefore crucial to experimentally manipulate densities.If organisms are ideal free distributed, and all available resources are used and converted into growth, an experimental doubling of natural densities is expected to result in an approximate halving of natural individual growth rates, independent of patch-quality and (associated) naturally occurring densities.However, if available resources are not fully used (i.e. if the strength of density dependence is lower), then the effect of doubling densities on individual growth rates will be smaller.Note that this scenario only holds in the absence of interference competition, a reasonable assumption for most soft-sediment benthic communities (reviewed by Peterson 1979).
In this study we investigated the role of species-specific feeding traits, habitat characteristics and season in shaping density dependence in shell growth within soft-bottom bivalve communities in a tropical intertidal ecosystem, Banc d'Arguin, Mauritania.By experimentally doubling natural densities of bivalve communities in two contrasting seasons, winter and summer, and by exploiting the spatial variability in habitat properties among sites within our study system (ranging from bare sandy to seagrass-covered muddy sediments), we tested for seasonal and/or local habitat effects on the magnitude of density dependence in shell growth of the three most abundant local bivalve species. .In addition to two filter-feeding bivalve species (i.e.Senilia senilis (Arcidae) and Pelecyora isocardia (Veneridae)), our study includes the first assessment of density dependence in shell growth of a bivalve species that mainly lives off carbon products provided by sulphide-oxidizing chemoautotrophic bacteria living inside its gills (i.e.Loripes orbiculatus (Lucinidae)) (van der Geest et al. 2014).In this symbiotic association, the bivalve host enhances chemosynthesis by its gill-symbionts by facilitating the supply of sulphide, carbon dioxide and oxygen to its gills.In exchange, the bacterial gill-symbionts fix carbon, fuelling their own energetic and biosynthetic needs, in addition to those of their host (Stewart et al. 2005).Assuming that the bacterial chemosynthesis of the food of 'chemosymbiotic' bivalves is dependent on the transport of resources (i.e.sulphide, oxygen and carbon dioxide) from the local environment (i.e.surrounding pore water) to the gill-symbionts by the host, we hypothesize that when experimentally increasing densities of chemosymbiotic bivalves, these resources may become locally depleted.Therefore, we predict density-dependent growth in chemosymbiotic L. orbiculatus.Moreover, assuming that our focal species obey an ideal free distribution, we predict that the effect of our density treatment on L. orbiculatus will be independent of local habitat characteristics (i.e.sediment grain-size and seagrass biomass).
In contrast to chemosymbiotic bivalves, filter feeders are less dependent on food produced in their immediate vicinity, but instead depend on water flow for delivery of their food produced elsewhere (Fréchette & Bourget 1985).As a result, competition for food may be low over a wide range of numerical densities.Yet, by reducing hydrodynamics (Fonseca et al. 1982), seagrass may reduce suspended food refreshment rates (Reusch & Williams 1999, Gonzalez-Ortiz et al. 2014) and physical disturbance of filter-feeding activity (Irlandi & Peterson 1991, Irlandi 1996).Moreover, lower predation risk within seagrass beds, due to seagrass providing refuges, may allow bivalves to be less vigilant and spend more time filterfeeding (Irlandi & Peterson 1991).Together, these factors may contribute to local food depletion by filter feeders inhabiting seagrass beds, when their natural densities are experimentally enhanced.Therefore, we predict that the negative effect of our density treatment on shell growth in our two focal filter-feeding species will increase with seagrass biomass.Assuming that the intensity of competition is greater during periods of resource shortage (Wiens 1977), and that resource availability may vary seasonally for both filterfeeding and chemosymbiotic bivalves, we predict an interaction between season and our density treatment on shell growth in all focal species.By evaluating the potential impact of feeding mode (filter-feeding versus chemosymbiotic) and environmental conditions (i.e.habitat characteristics and season) on the magnitude of density dependence in shell growth within seagrass-associated tropical bivalve communities, we aim to contribute to a better understanding of the factors that determine tropical bivalve population dynamics.This knowledge may help to identify threats posed by environmental change (e.g.rapid loss of seagrasses worldwide (Waycott et al. 2009)) and to guide conservation strategies.

Study site
The study area around the fishing village of Iwik (Fig. 1) is an accessible part of the intertidal area of , off the coast of Mauritania.This intertidal ecosystem is characterized by tidal flats of which ~80% are covered by seagrasses (mainly Zostera noltei Hornemann; Wolff & Smit 1990) that retain consistent aboveground biomass throughout the year (Vermaat et al. 1993, El-Hacen et al. 2018).The Banc d'Arguin is an important wintering site for migratory shorebirds, hosting almost two million individuals in winter (Oudman et al. 2017).Previous studies have indicated that the food web of this intertidal area is mainly supported by local benthic primary production (i.e.seagrass and microphytobenthos), with low contributions of phytoplankton, macrophytes and epiphytes to the food-web (Wolff et al. 1993b, Carlier et al. 2015).Being adjacent to the Sahara desert, Banc d'Arguin does not receive freshwater inflow from rivers and precipitation is limited to occasional thunderstorms that occur at irregular intervals, sometimes several years apart (Wolff & Smit 1990, van der Geest et al. 2014).
In the Iwik region, seawater temperature varies between ~20°C in winter (i.e.January) and ~30°C in late summer (i.e.September) and salinity ranges 40-44‰ (Wolff & Smit 1990, van der Geest et al. 2014).The tide is semi-diurnal and the tidal range is 1 to 2 m; maximum current speeds are about 1 m s -1 (Wolff & Smit 1990).

Experimental design
To avoid artefacts that may be imposed by commonly used enclosure experiments (for details see Text S1 in the Supplement), we used an experimental set-up to manipulate bivalve densities while keeping growing conditions as natural as possible, as will be described below.
The study area was divided into seven sub-regions, spread over an area of about 36 km 2 (Fig. 1).In each sub-region, we haphazardly assigned an annulus (hereafter named station) with an outer radius of 200 m and an inner radius of 100 m.Within each station, 10 sampling sites were randomly selected.Our sampling procedure thus yielded 70 sampling sites.
In autumn, between 13 October and 5 November 2007, and spring, between 13 April and 4 May 2008, benthic samples were taken at each sampling site during low tide.The samples were obtained by pushing a 30 cm diameter PVC ring (15 cm high) into the sediment and collecting all sediment within the ring to a depth of 10 cm, using a small shovel.To speed up the process of sorting benthic samples, we sieved the samples over a 2.8-mm mesh.By doing so, we may have missed bivalves < 2.8 mm.Yet, the density of bivalves < 2.8 mm was relatively small compared to the density of larger bivalves (van der Geest 2013).Moreover, being more fragile, these smaller bivalves would have been unlikely to survive the experimental handling and transportation.
The retained material on the sieve was transported to the field station near Iwik.Here, all living bivalves were separated from the matrix of seagrass remains and shell fragments.
From all clams of S. senilis, P. isocardia and L. orbiculatus, we haphazardly took 1 to 18 individuals per species of which shell heights (H1) were measured to the nearest 0.1 mm.This range from 1 to 18 individuals was the result of varying bivalve sizes and densities of our focal species in our benthic samples.Subsequently, these bivalves were individually marked with uniquely labelled glue-on shellfish tags (type 'FPN 4 mm circle tag' or type 'FPN 8x4m oval tag', Hallprint, Australia).Until redeployment, all bivalves were covered with a layer of seawater-moist Zostera noltei leaves to protect them from desiccation and stored at ~25°C.
Within 24 hours after collection, all bivalves were relocated in their natural environment as described below.
As sampling inevitably leads to disruption of sediment structure, all bivalves (both tagged and untagged, often comprising several species) from a benthic sample were relocated altogether at an undisturbed site 5 m away from the original sampling site.The relocation site was marked with two PVC sticks that were placed 50 cm apart.In the middle between the two sticks, a 30 cm diameter PVC ring was placed on top of the sediment.The bivalves were spread out over the enclosed area and gently pushed into the sediment to a depth of ~1 cm, after which the PVC ring was removed.As the minimum average size of a bivalve patch in our study area had been estimated to be ~10 m (Oudman et al. in press), we expected the composition of the bivalve community to be similar across a 5 m distance.Thus, by relocating all bivalves collected in a benthic sample at 5 m distance from their original sampling site over a surface area equal to the sampled surface area at the original sampling site, we roughly doubled the density of the total bivalve community at the relocation site.
For the focal bivalve species of which individuals were collected for the density treatment, we collected 1 to 8 additional clams in the vicinity of the sampling sites to use as control treatments.These clams were also measured and individually tagged as described above.In contrast to the treatment where bivalve densities were approximately doubled, these tagged clams were each placed 1 m apart and within 3 m of the site where bivalve densities were doubled.With mean densities of S. senilis, P. isocardia and L. orbiculatus being 242, 219 and 16 individuals per m 2 in bare sediments and 59, 316 and 339 individuals per m 2 in seagrass-covered sediments in our study area (Honkoop et al. 2008), bivalve densities per m 2 were only marginally changed by the addition of just a single individual.Therefore, we considered these sites as control treatments with natural (unmanipulated) densities.These control sites were again marked by two short PVC sticks placed 50 cm apart, the middle point between two short PVC sticks marking the spot where the tagged clam was relocated.In total, 1431 clams were tagged (S. senilis, N = 529; P. isocardia, N = 357; L. orbiculatus, N = 545, Table 1).
Approximately half a year later, in spring (13 April -4 May 2008), and in autumn ( 20October -14 November 2008), we sampled our experimental relocation sites in a similar way as described above, with the only difference that control sites, to which only one tagged clam was allocated, were sampled with a smaller (15 cm diameter) sediment core to reduce sampling effort.As such, we investigated the effect of our density treatment on shell growth in both winter (i.e.October 2007-April 2008) and summer (April 2008-November 2008).
Each benthic sample was sieved over a 2.8-mm mesh and the material retained was put in a plastic bag and transported to the scientific field station near Iwik.Here, all living clams were sorted per sample and when a tagged clam was recovered, its shell height (H2) was measured again (precision 0.1 mm).At each sampling, tagged clams were either found alive, found dead as empty shells, or missing.Missing clams would have been a consequence of either sampling error, emigration, removal by scavengers, post-mortem transport, or depredation.In some occasions we recaptured single tags in our samples that were not attached to a clam anymore.
Because we could not determine whether these tags belonged to clams that were alive, dead or missing from our study plot, we labelled their fate as 'unknown'.
As the tagging procedure involves physical handling and temporary removal of clams from their natural environment, it may cause retarded growth due to stress, in particular in L. orbiculatus, the species with somewhat fragile shells (van der Geest et al. 2011).Yet, as this would affect the growth of clams in both the control and density treatments, it will not invalidate any conclusions about the factors that determine density-dependent shell growth in our focal bivalve species.

Site-specific habitat characteristics
To investigate if habitat characteristics (i.e.seagrass biomass and grain-size of the sediment) affect the magnitude of density-dependent growth, half-way during the experimental interval (i.e. between 13 April and 2 May 2008) at each site a seagrass core (7 cm, internal diameter) and a sediment core (internal diameter, 2 cm) was taken to a depth of 10 cm.The content of the seagrass core was sieved over a 500-µm mesh.The material retained on the sieve was stored in a plastic bag, frozen at -18 °C and transported to The Netherlands, where for each sample all living seagrass parts (i.e.leaves, rhizomes and roots) were sorted.The ash-free dry mass (AFDM) of all living seagrass parts was determined via the loss-on-ignition method.
Samples were dried at 60°C for a minimum of 72 h, weighed and then incinerated at 550°C for 4 h after which the remaining ashes were weighed again.The difference between the first and the second weights gives the AFDM of the living seagrass parts in the sample (AFDM seagrass, in g m -2 ).
Content of the sediment core was stored in a plastic bag, frozen at -18°C and in The Netherlands samples were freeze-dried and grain-size distribution of each sample was determined using a particle-size analyser (Beckman Coulter Model LS 230).From the grainsize distribution the median was calculated (median grain-size, MGS).Given that MGS depends on wave exposure and factors that attenuate hydrodynamics (Fonseca et al. 1982, Paterson & Black 1999), we used MGS as a proxy for local water flow conditions, with larger MGS reflecting higher water flow (i.e.refreshment rates) at that site.

Shell growth rate
For the three focal species, the mean initial size of recaptured bivalves varied between our density treatments (for details see Table 1).To remove the effects of initial size on the magnitude of individual growth increments, we fitted Von Bertalanffy's growth function (VBGF) to our data, a commonly used equation when modeling indeterminate shell growth.
In this function, instantaneous growth rate dH/dt declines with an increase in shell height Ht in the following way: where  ∞ is a constant representing the mean maximum shell height and k is the growth coefficient representing the intrinsic rate (in day -1 ) at which asymptotic shell height is approached.To estimate k from tag-recapture data, the traditional VBGF must be modified using the derivation of Fabens (1965) increment model: where k is the estimated growth coefficient and H1 and H2 are defined as the shell heights at time of marking (t1) and recapture (t2), respectively, ∆t as the time interval in days (i.e. t2 − t1) and  ∞ is the mean maximum shell height.Rewriting this equation gives: We used the obtained growth coefficient k (day -1 ) for each recaptured clam as a proxy for shell growth rate in our analyses.

Statistical analysis
As sediment median grain-size (MGS) and seagrass biomass were negatively correlated (Pearson's r = -0.49),we performed a principle component analysis and used the first principle component (PC1) as a continuous variable to describe the habitat characteristics at a specific sampling site (Fig. 2).PC1 explained 74% of the proportion of the variance in sediment MGS and seagrass biomass.
For each focal species separately, we used linear mixed effects models to investigate whether density treatment (control or initial bivalve density doubled), season (winter or summer) and habitat (PC1) and their two-way interactions explained variation in growth coefficient k.Due to the nested structure of the data, station and site (nested within station) were included as random effects.We compared all possible combinations of these explanatory variables.Model selection is based on Akaike information criterion adjusted for small sample size (AICc) (Burnham & Anderson 2002).Parameter estimates and approximate 95% confidence intervals of the most parsimonious model are reported, with the most parsimonious model being the model with the fewest parameters within 2 ΔAICc of the top model (Burnham & Anderson 2002).Per focal species, the R 2 LR goodness-of-fit measure was calculated from the most parsimonious linear mixed effects model output by using a loglikelihood ratio test as described by Magee (1990).
We tested for heterogeneity in the residuals following the procedure described by Zuur et al. (2009), by comparing models that described the variance as different functions of the explanatory variables.We determined the appropriate value for  ∞ for each species by gradually increasing the value of  ∞ (starting from the H1 of the largest clam in the dataset) until the most parsimonious variance function no longer includes H1 (for details of this method, see Text S2 in the Supplement).This turned out to be at 76.3, 17.1 and 11.4 mm for S. senilis, P. isocardia and L. orbiculatus, respectively.Using these values for  ∞ , there was still evidence that variance was season-dependent in all three focal species for which we corrected using a 'varIdent' structure (Zuur et al. 2009), which allows different error variances for different factor levels.This variance structure was retained when investigating the statistical support for any of the fixed effects.As there is some individual variation around  ∞ , we performed a sensitivity analysis with respect to the maximum value of  ∞ (for details see Text S3 in the Supplement).
All analyses were performed in program R (R Core Team 2018, version 3.4.4).For linear mixed effects models, the R-package 'nlme' (Pinheiro et al. 2015) was used.Models with all possible combinations of main effects and their two-way interactions were compared simultaneously based on AICc, using the R-package 'MuMIn' (Bartoń 2016).
Live recaptured clams (N = 338) were distributed over 68 of the 70 experimental sites.
The habitat characteristics (as described by PC1) of the sites where tagged clams were recaptured differed significantly among the three focal species (one-way ANOVA, F2,119 = 5.84, P = 0.004).A post-hoc Tukey test showed that the habitat characteristics at the sites occupied by S. senilis differed significantly from those occupied by P. isocardia (P = 0.03) and L. orbiculatus (P = 0.007), with recaptured S. senilis being significantly more restricted to bare sandy sediment sites with relatively low PC1 values, compared to clams of P. isocardia and L. orbiculatus which were recaptured in more muddy seagrass-covered sites with relatively high PC1 values (Fig. 3).There was no significant difference between habitat characteristics at the sites occupied by P. isocardia and those occupied by L. orbiculatus (P = 0.99).
For S. senilis, the most parsimonious model to account for variation in growth rate included an effect of season only (Table 2), which explained 39% of the variation in shell growth (R 2 LR = 0.39).Growth coefficient k was about twice as high in summer than in winter (Table 3, Fig. 4a).Although the model including both season and density treatment as fixed effects was best supported (Table 2), it required an extra parameter without reducing AICc by 2 points (Table 2).This implies that there is only limited support for an effect of density (Burnham & Anderson 2002), as also indicated by the 95% confidence interval that included zero (βdensity = 3.87*10 -5 , 95% CI = -0.80*10 - -8.54*10 -5 ).
The most parsimonious models to account for variation in growth rate in P. isocardia and L. orbiculatus included effects of season and density treatment (Table 2) (R 2 LR = 0.43 and R 2 LR = 0.32, respectively).For both P. isocardia and L. orbiculatus, k was higher in summer than in winter and was lower in the density treatment than in the control (Table 3, Fig. 4b,c).
While the best-supported model for P. isocardia included an effect of habitat, this extra parameter reduced AIC by only 0.02 points, and the 95% CI included zero (βhabitat = -2.57*10- 4 , 95% CI = -5.30*10 - -1.63*10 -5 ), implying only limited support for a habitat effect.There was no support for the statistical interactions between our density treatment and habitat characteristics or season in any of the focal species (Table 2).These results were insensitive to the value of maximum shell height  ∞ across a wide range of values for  ∞ (P.isocardia, 17.1-23 mm; L. orbiculatus 11.4-14 mm, S. senilis 76.5-81.7 mm) (for details see Text S3 in the Supplement).

Discussion
By quantifying the combined effects of density treatment, season and habitat on shell growth of bivalves from two different feeding guilds, this experimental study yields important insights in the factors that determine density dependence in shell growth of tropical bivalves.
We show that shell growth was negatively density dependent in filter-feeding P. isocardia and chemosymbiotic L. orbiculatus, the two species that mainly inhabit seagrass sediments, but not in filter-feeding S. senilis, a species that dominates bare sandy sediments (Fig. 3).We found no season or habitat-specific effect of density on shell growth.These results suggest that the bivalve community of seagrass-covered sediments was closer to carrying capacity than that of adjacent bare sediments, regardless of species-specific feeding mode or season.

Intra-versus interspecific competition
As we manipulated the density of the total bivalve community, we estimated the joint effects of intra-and interspecific competition on shell growth in our three focal bivalve species.In view of the different feeding modes (filter-feeding vs. feeding mainly on chemosynthetically produced food), and the spatial segregation of the three numerically dominant bivalve species studied here (S. senilis being restricted to more bare sandy sediments and P. isocardia and L. orbiculatus mainly inhabiting muddy seagrass sediments; Fig. 3), intraspecific competition may have been more prevalent than interspecific competition.Yet, L. orbiculatus has recently been shown to grow faster when P. isocardia is depleted, which was attributed to L.
orbiculatus being a mixotroph, thus potentially competing for suspended resources with filterfeeding P. isocardia (van Gils et al. 2012).Moreover, exploitative interspecific competition for oxygen cannot be ruled out (Ferguson et al. 2013), especially given the detected oxygen deficit due to benthic respiration at night in the seagrass-covered sediments of Banc d'Arguin (Clavier et al. 2011).Irrespective of the relative importance of intra-versus inter-specific competition, density dependence in shell growth was stronger in bivalve communities of seagrass-covered sediments compared to those in bare sediments.

Density-dependent growth in filter feeders
Our density treatment strongly affected shell growth in P. isocardia, while there was only limited support for an effect of our density treatment on shell growth in S. senilis.Although we cannot exclude the possibility that species-specific morphological or life history traits play a role in explaining the differential effect of density on filter-feeding S. senilis and P.
isocardia, the difference in the main habitats occupied by the two species likely contributes to this effect.S. senilis mainly lives in bare sandy sediments whereas P. isocardia occupies the more muddy seagrass-covered sediments (Honkoop et al. 2008;Fig. 3 of this study).
Filter feeders depend on the flow of water to supply them with food (Fréchette & Bourget 1985).When flow velocities of the water column and subsequent food refreshment rates are low, filter feeders can locally deplete their food (Reusch & Williams 1999, Gonzalez-Ortiz et al. 2014).Hence, the observed negative correlation between sediment MGS and seagrass biomass (Fig. 2), suggests that inside seagrass canopies water flows (i.e.food refreshment rates) are reduced (Paterson & Black 1999).This would result in enhanced competition for suspended food particles among the filter feeders of the seagrass-covered sediments of Banc d'Arguin.
Alternatively, one could argue that the lack of density-dependent growth in bivalves of bare sediments (i.e. S. senilis), is simply caused by the disappearance of the density treatment (doubling densities) by bivalves moving away from the experimental plots more freely in bare sediments than in seagrass-covered sediments.However, as the percentage of missing tagged S. senilis clams at resampling was higher at the control sites compared to sites where clam densities were doubled (45% vs. 43% in spring, and 57% vs. 45% in autumn; Table 1), this was unlikely the case.

Density-dependent growth in chemosymbiotic bivalves
To our knowledge, this study is the first to reveal density-dependent growth in chemosymbiotic bivalves.Chemosymbiotic bivalves generally dominate the macrobenthic infauna of seagrass sediments (Honkoop et al. 2008, van der Heide et al. 2012), where they can reach densities of up to 4000 individuals per m 2 (van der Geest et al. 2011).Here, these bivalves and their sulphide-oxidizing gill-symbionts can be considered part of a nested symbiosis with seagrasses, which may be essential to the health and ecological success of seagrasses (van der Heide et al. 2012); while the bivalve-bacteria consortium profits from sulphide that is indirectly provided by seagrasses, and from oxygen released by seagrass roots, the removal of toxic sulphide by the bivalve-bacteria consortium stimulates growth of the seagrasses (van der Heide et al. 2012).In light of this tripartite mutualism, the observed density-dependent growth in chemosymbiotic L. orbiculatus, suggests a limit to their capability to detoxify sulphide from seagrass sediments, which would have major implications for seagrass community functioning and persistence.For example, negative density-dependent growth in chemosymbiotic bivalves could have contributed or even initiated the recently described breakdown of the mutualism between seagrass and chemosymbiotic L. orbiculatus, which accelerated landscape-scale intertidal seagrass collapse in our Banc d'Arguin study system (de Fouw et al. 2016).Clearly, whether and how densitydependent growth in chemosymbiotic bivalves affects the resilience of this seagrass-bivalve mutualism provides scope for future studies.

Effects of environmental conditions on density-dependent growth
We found pronounced seasonal variation in shell growth in all three bivalve species (Tables 2    & 3, Fig. 4).This is consistent with previous studies at Banc d'Arguin that investigated seasonality in shell growth of S. senilis (Lavaud et al. 2013)    * 100% using the parameter estimates in Table 3).This may indicate that, compared with chemosymbiotic L. orbiculatus, seasonality in (phototrophic) food availability was more intense for filter feeders, especially those inhabiting seagrass sediments (i.e.P. isocardia).
The lack of support for a statistical interaction between season and density on growth in P. isocardia and L. orbiculatus did not confirm our prediction that the density treatment would have a stronger effect during periods of reduced growth conditions (i.e., in winter).We suggest that this counterintuitive result may be due to the seasonality in predation pressure at our experimental sites.Being responsible for about 80% of all mollusc consumption by vertebrate predators (van Gils et al. 2012), the migratory red knot Calidris canutus canutus is the main molluscivore predator at Banc d'Aguin, where their diets mainly consist of P.
isocardia and L. orbiculatus (Onrust et al. 2013).These birds are able to deplete their favourite P. isocardia and L. orbiculatus prey over the winter, before they return to their Arctic breeding grounds in late spring (Ahmedou Salem et al. 2014).Since red knots preferably forage at sites with high densities of P. isocardia and L. orbiculatus (van Gils et al.

2015)
, and since we have not excluded predators from our experimental plots, the relatively high predation pressure in winter may have reduced the effect of our density treatment on shell growth in P. isocardia and L. orbiculatus, which may have compensated for the reduced growth conditions in this season, eventually resulting in the lack of an interaction between season and our density treatment.Indeed, the percentage of missing tagged P. isocardia clams at sites where local bivalve densities were experimentally doubled was higher in spring (74% clams missing) compared to autumn (61% clams missing), which may reflect higher predation rates on P. isocardia in winter (Table 1).
Whereas the lack of support for a statistical interaction between density treatment and habitat characteristics confirmed our prediction for chemosymbiotic L. orbiculatus, we had expected a positive relationship between seagrass biomass and density dependence in growth, hence a statistical interaction between density treatment and habitat characteristics, for the filter-feeding S. senilis and P. isocardia, as refreshment rates of suspended food generally decrease as seagrass biomass increases (e.g.Reusch & Williams 1999, Gonzalez-Ortiz et al. 2014).That there was no support for such a statistical interaction in the two filter-feeding species could be explained by the fact that S. senilis and P. isocardia occupied only a limited range of habitats (bare sandy sediment and seagrass-covered mud respectively, Fig. 3), making it less likely to statistically detect any within-species habitat effects.

Density-dependent regulation
Density-dependent regulation is essential for the long-term persistence of populations (Murdoch, 1994); it can reduce susceptibility to environmental fluctuations (Anderson et al. 2008), resulting in lower risk of extirpation or extinction.As fecundity and survival in bivalves generally increase with size (e.g.Paine 1976, Peterson 1986), the observed density dependence in shell growth of L. orbiculatus and P. isocardia is likely to result in negative density dependence in reproduction and survival rates driving regulation of their populations through negative feedbacks.If true, we would expect rather constant biomass densities in L.
orbiculatus and P. isocardia, and more variable densities in S. senilis over time.Indeed, a comparison of biomass densities at the tidal flats of Banc d'Arguin in 1986, 1988and 2007 showed that L. orbiculatus and P. isocardia biomass were very constant in these three periods, varying between 1.1-2.6 g ash-free dry mass (AFDM) m -2 in L. orbiculatus and 0.1-0.9g AFDM m -2 in P. isocardia, while biomass densities of S. senilis were highly variable, varying from 8.1 g AFDM m -2 in 1986 to 0.8 g AFDM m -2 in 1988 up to 20.3 g AFDM m -2 in 2007 (Wolff et al. 1993a, Jansen et al. 2008, Wolff & Michaelis 2008).

General implications
This study suggests that tropical bivalve communities of seagrass-covered sediments are closer to carrying capacity and presumably more 'regulated' by density-dependent processes than those of adjacent bare sediments.Seagrasses are ecosystem engineers (Jones et al. 1994) in the sense that they create, modify and maintain their own habitat by causing changes in biotic and abiotic conditions that modulate the availability of resources to themselves and other species (e.g.Bos et al. 2007, Folmer et al. 2012).As a result of ecosystem engineering, seagrasses create more stable and predictable environmental conditions for those species that depend on them (Bertness & Leonard 1997).We suggest that in tropical realms like Banc d'Arguin, where seagrasses generally retain consistent aboveground biomass throughout the year (Duarte 1989, Vermaat et al. 1993, El-Hacen et al. 2018)

Fate
of clams for our three focal bivalve species over two time periods as a function of density 653 treatment.Winter = October 2007-April 2008, Summer = April 2008-November 2008, 654 Ntagged = number of clams tagged, Nrecap = number of tagged clams that were recaptured, A = 655 percentage recaptured alive, D = percentage recaptured dead, M = percentage missing, U = 656 percentage of clams of which the fate was unknown, H1 tagged = mean initial shell height (mm) 657 of tagged clams, H1 recap = mean initial shell height (mm) of tagged clams that were 658 recaptured, H2 recap = mean shell height (mm) of tagged clams that were recaptured.

Fig. 4 .
Fig. 4.The effect of density treatment on growth coefficient k (day -1 ) in both winter (i.e.

Table 2 .
Model selection results for growth coefficient k (day -1 ) as a function of density treatment (d), season (s) and habitat (h; modeled as a continuous PC1 variable) and all possible two-way interactions per focal bivalve species (i.e.Senilia senilis, Pelecyora isocardia and Loripes orbiculatus).The most parsimonious model is shown in bold.K denotes the number of parameters.Only models with a model weight of > 0.02 are shown.