Burrowing Behavior of a Deposit Feeding Bivalve Predicts Change in Intertidal Ecosystem State

Behavior has a predictive power that is often underutilized as a tool for signaling ecological change. The burrowing behavior of the deposit feeding bivalve Macoma balthica reflects a typical food-safety trade-off. The choice to live close to the sediment surface comes at a risk of predation and is a decision made when predation danger, food intake rates or future fitness prospects are low. In parts of the Dutch Wadden Sea, Macoma populations declined by 90% in the late 1990s, in parallel with large-scale mechanical cockle-dredging activities. During this decline, the burrowing depth of Macoma became shallow and was correlated with the population decline in the following year, indicating that it forecasted population change. Recently, there has been a series of large recruitment events in Macoma. According to the food-safety trade-off, we expected that Macoma should now live deeper, and have a higher body condition in association with this change in depth of living. Indeed, we observed that Macoma now lives deeper and that living depth in a given year forecasted population growth to the next year, especially in individuals larger than 14 mm. As living depth and body condition were strongly correlated in individuals larger than 14 mm, larger Macoma could be living deeper to protect their reproductive assets. Our results confirmed that burrowing depth signals impending population change and, together with body condition, can provide an early warning signal of ecological change. We suggest that population recovery is being driven by improved intertidal habitat quality in the Dutch Wadden Sea, rather than by the proposed climate-change related effects. This shift in ecosystem state is suggested to include the recovery of diatom habitat in the top layer of the sediment after cockle-dredging ended.


INTRODUCTION
Terrestrial and marine environments are currently changing at unprecedented rates and scales (Vitousek et al., 1997;Halpern et al., 2008). Understanding how human induced changes affect animal populations is a continuous challenge, especially when the drivers of change follow indirect pathways and/or are undetectable (Wong and Candolin, 2014). For many animals, it is well known that behavioral changes can signal habitat alteration (Brown and Kotler, 2004;Kotler et al., 2007). This is because a change in behavior can relatively quickly improve the prospects of surviving and reproducing in changing habitats (Wong and Candolin, 2014). However, few studies have examined the capacity for behavioral indicators to forecast changes in population growth at the landscape scale (Lima and Zollner, 1996;van Gils et al., 2009).
In the Dutch Wadden Sea, as in most coastal systems of the human inhabited world, intertidal ecosystems have been affected by human induced perturbations (e.g., fisheries and climate change; Lotze et al., 2006;Eriksson et al., 2010;Wolff, 2013). One of the most dramatic changes in recent years in the Dutch Wadden Sea was the decline (∼90%, from late 1990s to 2006) of the Baltic tellin, Tellinidae (Macoma balthica, Linnaeus, 1758;Piersma et al., 2001;van Gils et al., 2006;Kraan et al., 2007). Debate ensued as to why Macoma populations declined with competing ideas being: (1) negative effects of climate change on recruitment success, partially caused by climaterelated changes in the temporal presence of crustacean predators of bivalve spat (Philippart et al., 2003;Beukema et al., 2009, and (2) adverse effects of large-scale mechanical cockledredging (Cerastoderma edule) on sediment characteristics and benthic communities Kraan et al., 2007), possibly affecting the feeding opportunities of Macoma (van . As the debate went on, Macoma might have already signaled reasons for its decline. Specifically, Macoma lived shallow during the population crash period, but deep prior to the population decline . Furthermore, population growth and living depth were correlated, which indicated that burrowing behavior could be used as a predictor of population change . Depending on its internal physiological state, an animal must decide whether to face the risk of predation and obtain a higher food intake, or evade predation risk and live on a lower food intake (Houston et al., 1993;McNamara and Houston, 1996). This is the so-called "food-safety trade-off " (Lima and Dill, 1990;Lima, 1998a,b). For burrowing bivalves, living shallow in the sediment is a "risky behavior, " as there is a greater risk of being caught by predators Zaklan and Ydenberg, 1997). However, as shallow living is associated with a higher food intake, when body condition is low this is a risk worth taking (Zwarts, 1986;Lin and Hines, 1994;Zwarts et al., 1994). Conversely, living deeper is "safer" (Zwarts and Wanink, 1989) but comes at the cost of a lower food intake (Lin and Hines, 1994;de Goeij and Luttikhuizen, 1998). Furthermore, as the marginal value of energy for animals in a good body condition is low (Brown and Kotler, 2004), individuals with more energy stores should choose for safety and have a lower intake rate, a principle called the "asset protection principle" (Clark, 1993). Based on the asset protection principle, the collective behavior of many individuals should reflect the future prospects of the population like population growth .
In recent years, Macoma has started to recover across the Dutch Wadden Sea (Compton et al., 2013). This gives us an opportunity to verify and extend the ideas presented by van . Here we analyze 17-years of data, based on a spatially extensive and consistent sampling design over summer, to examine whether Macoma's depth of living is a trait signaling ecological change. We expected that Macoma should be living deeper compared with the period of population decline and that this should be correlated with yearly population growth estimates in all 17-years. An association between deeper living and a higher body condition would confirm that Macoma protects its future reproductive assets. We will explore whether living depth can be used to forecast ecological change.

MATERIALS AND METHODS
Sampling of Macoma was carried out across 14775 ha of intertidal flats in the western Dutch Wadden Sea (53 • 10 ′ -53 • 20 ′ N, 4 • 50 ′ -5 • 30 ′ E, Figure 1) in the summer months of June-August from 1997 to 2014. Over these 17-years, sampling was conducted using a grid-based sampling design (Figure 1, 500-m distance-grid). To estimate the population parameters of Macoma balthica, the same sampling stations that were revisited annually were used (see Table 1, van Gils et al., 2009;Compton et al., 2013). A grid-based sampling design provides a high statistical power for accurately determining species abundances (van der Meer, 1997; Bijleveld et al., 2012). Sampling was conducted over a period of 3 months and sampling stations were not revisited in the same order each year.
At each sample station, samples were taken using a benthic sediment core (total area sampled at each point: 0.0177 m 2 by foot and 0.0174 m 2 by boat). To estimate the living depth of Macoma, each sediment sample was partitioned into a top (4 cm) and bottom fraction (maximum depth of 30 cm) before sieving the sediments (1 mm sieve) and counting the Macoma in each fraction. The division of the core is based on the knowledge that the main shorebird predators of Macoma can only probe up to ∼4 cm in depth (Zwarts andWanink, 1989, 1991;. All samples were brought back to the laboratory for further analysis. In the laboratory, Macoma were recounted, measured, and weighed. After measuring the shell length of Macoma (∼0.1 mm accuracy), body condition of each individual was estimated by separating the body flesh from the shell, and placing the flesh into a crucible and drying it for 48 h. Once the flesh was dried, the crucibles with the dried flesh were placed into a desiccator to cool and then weighed to the nearest 0.1 mg. These crucibles were then placed into an incinerator where they were burnt at 560 • C for 5 h. After burning, the crucibles with the ashed flesh were cooled in a desiccator and then weighed. The flesh ash-free dry mass (AFDM flesh ) of an individual Macoma was determined by subtracting the flesh dry mass from the flesh ash mass. With the laboratory work for 2014 not yet finished, we used the fraction of Macoma living in the top part as assessed in the field, separately for two size classes, i.e., small (<10 mm) and large individuals (≥10 mm).

Analysis
To examine how the densities of adults and recruits changed over time, we described the stock-recruitment relationship. The shell length of the recruits was defined based on an examination of a subset of aged Macoma in this area. Macoma can be aged by counting growth rings on the shell (Lammens, 1967;   Cardoso et al., 2007). Based on the subset of aged Macoma, we observed that new recruits, i.e., bivalves with no growth ring, were consistently smaller than 10 mm in length. We thus defined recruits, in year t+1 , as individuals with a shell length smaller than 10 mm. Recruits were defined based on visual inspection in 2014. Stock densities in year t were estimated by excluding the recruits from the data (shells <10 mm), and densities of the recruits in year t+1 were estimated by excluding the stock (shells ≥ 10 mm) from the data. In the cases when either recruits or stock were excluded and nothing was observed, these samples were included in the density calculations as absences. Also note that the stock-recruitment relationship stayed the same even when maximum recruitment sizes were 6, 8, or 10 mm.
To describe the trends in the population growth of Macoma, we examined changes in the per capita population growth rate (r), which describes the magnitude and direction of population change for a given time period, as a fraction of the initial population. It provides an integral measure of growth and mortality processes in a population (Sibly et al., 2003). For all shell lengths of Macoma, population growth rates (r) were calculated by dividing the annual average density (D) in year t+1 by the annual average density (D) in year t and then logtransforming this ratio: To examine the depth of living in Macoma, the fraction in top was estimated by dividing the density of Macoma in the top 4 cm (D top ) by the total density of Macoma (in top and bottom) at each sampling station for a specific size range of Macoma at sites where this information was recorded: Although, it is known that Macoma larger than 10 mm and smaller than ∼18-20 mm do not differ in living depth over the summer period (Zwarts and Wanink, 1991 Macoma smaller than 10 mm were excluded, as they cannot live deep due to physical limits set by their siphon length (see Appendix 1 in Supplementary Material). As larger individuals (>18 mm) are known to have a higher incidence rate of parasite infection relative to smaller Macoma (Swennen and Ching, 1974;Pekkarinen, 1984;Edelaar et al., 2003), they often live closer to the sediment surface due to a low body condition associated with parasitism (see Appendix 1 in Supplementary Material and P. de Goeij pers. comm.). Thus, we excluded individuals above 18 mm from the analyses. For each of these size ranges, changes in living depth were described by calculating the annual average fraction of Macoma living in the top 4 cm of the sediment. The relative body condition in Macoma was estimated by calculating the body mass index (BMI) of each individual. The body mass index (BMI) was defined as the flesh ash-free dry mass (AFDM, mg)/shell length 3 (cm 3 ).
To deal with dependencies in the data, e.g., spatial autocorrelation, we ran all of the above calculations using a bootstrap sampling approach. Bootstrap sampling uses sampling with replacement to generate a sampling distribution of an estimate (Davison and Hinkley, 1997). Here, we randomly subsampled the data, with replacement, 1000 times. From this subsample we estimated the mean annual population growth, depth of living and body condition. The 1000-bootstrap estimates were used to estimate 95% confidence intervals in the representations of the trends over time. The confidence intervals were defined based on the percentile method (Davison and Hinkley, 1997). In estimating population growth, we subsampled 200 stations per bootstrap (∼50% of all stations visited). In the case of living depth, 10 stations were sampled for each bootstrap and for body condition 40 individuals were sampled per bootstrap.
To estimate whether mean annual population growth rate and mean annual depth of living were correlated in the three size groups, we ran linear regressions using (a) the raw data and (b) the bootstrapped estimates of mean population growth and mean depth of living. Regressions between the bootstrapped estimates provided a level of confidence in this fitted relationship by providing frequency distributions of the (1) slope magnitude and direction and (2) fit of the estimates, i.e., R 2 . The same procedure was followed for the association between mean depth of living and mean body condition in the years 1997-2013 for the three size groups. The percentage of times that the regressions were significant was calculated from the number of times the pvalue fell below 0.05 in the 1000 bootstraps. The residuals of these regressions were checked, but there were no patterns suggesting dependencies or missing covariates.
All analyses were carried out in R (version 3.0.3, R Core Team, 2014) within R Studio (version 0.98.932). The package Frontiers in Ecology and Evolution | www.frontiersin.org ggplot2 was used to create Appendix 1 in Supplementary Material (Wickham, 2009).

RESULTS
Densities of Macoma in the western Wadden Sea declined from ∼200 n/m 2 in 1998 to 25 n/m 2 in 2011 (Figure 2A). During this period, recruitment also remained low (Figure 2B). In the years 2012 to 2014, mean annual recruit densities increased up to 264 n/m 2 due to an influx of recruits. Notably, these recruitment events were not associated with higher adult densities, as adult densities were low from 2011 to 2013 (Figure 2B, <30 n/m 2 ).
The annual average estimates of population growth, and their bootstrapped confidence intervals, showed that from 1997-2005 population growth in year t was negative or constant (r = 0, Figure 3A, see Appendix 2 in Supplementary Material for size range 10-18 mm), but that since 2006 population growth turned mainly positive (r > 0). The bootstrap confidence intervals showed that the overall changes in Macoma's population growth remained consistent even when 50% of the sampling stations were randomly selected.
Over the 17-years of this study, it appears that irrespective of shell length the fraction of shells in top, for shells larger than 10 mm, followed a trend over time (Appendix 1 in Supplementary Material). We observed that Macoma's living depth changed in a unimodal way over time across all three defined age groups. Specifically, the annual fraction of Macoma living in the top was initially low from 1997 to 2000 (fraction: 0.4-0.5), then became high from 2001 to 2006 (fraction: 0.8-0.9) and then started to decline after 2006 (fraction: <0.75; Figure 3B; see Appendix 2 in Supplementary Material for size range 10-18 mm). The body mass index of Macoma was variable over time, but with a trend for the body mass index to be initially high and variable from 1997 to 2000 (BMI: ∼13), then variable but generally becoming lower from 2001 to 2007 (BMI: ∼11) and then increasing from 2007 (BMI: >13) ( Figure 3C, see Appendix 2 in Supplementary Material for size range 10-18 mm). Bootstrapped confidence intervals, estimated based on random subsamples of the data, showed that these trends in living depth and body mass index were not due to spatial dependencies in the data.
There was a negative correlation between population growth of Macoma and the fraction living in the top for all size groups not only for the raw data ( Figure 4A 14-18 mm, Table 2, see Appendix 3 in Supplementary Material for size range 10-18 mm), but also for the majority of correlations based on the bootstrapped data (Table 3 all size groups). However, the fit of this relationship was weak for the raw data for all three groups (R 2 < 0.2, Table 2), with low R 2 values for the correlations based on the bootstrapped data (Table 3, <5% of the correlations had an R 2 of 0.3). There was a weak correlation (p = 0.09) for the size group of 14-18 mm between population growth and fraction in the top, but not for the other two size groups ( Table 2). The larger size group had the highest percentage of significant correlations (∼13% , Table 3), whereas in the other two size groupings the correlations of the bootstrapped data were less often significant (∼5% of cases, Table 3).
The correlation between the mean annual fraction of Macoma living in the top and mean annual body mass index was negatively correlated for all three size groupings for the raw data, as well as the majority of correlations based on the bootstrapped data ( Figure 4B 14-18 mm, Table 2 all size groups). The fit of this relationship was best for the larger size group (14-18 mm, R 2 = 0.4, Table 2), and was weak for the smaller size group (10-14 mm, R 2 = 0.006, Table 2). Correlations of the bootstrapped data also showed that the R 2 of the larger size group was stronger (32% with R 2 < 0.3, Table 3) than for the smaller size group (0.6% with R 2 < 0.3, Table 3). The correlation between the mean annual fraction of Macoma living in the top and mean annual body mass index was significant for the larger size group, and almost significant for all sizes, and not significant for the smaller size group based on the raw data (10-14 mm, Table 2). The percentage of times that this correlation was significant, based on

Model 1 (M1) was the population growth (growth) as a function of the fraction living in the top (fraction in top). Model 2 (M2) was the fraction in top as a function of body mass index (bmi). The F-values (which always had 15
• of freedom) and p-values are given. The magnitude and direction of the slope and fit (R 2 ) are also shown. the bootstrapped data, was also highest for the largest size group (Table 3).

DISCUSSION
The return to positive population growth in Macoma was associated with recent large recruitment events from 2012 to 2014. Coincidentally, Philippart et al. (2014) reported an increased probability of encountering Macoma larvae in the water layer from 2011 to 2012 compared to the years of 2006 to 2010. The observation of relatively large recruitment events in a continually warming Wadden Sea (Beukema et al., 2014, Appendix 4 in Supplementary Material) contradicts the current idea that a warming Wadden Sea and the associated epibenthic predator increases are driving recruitment failure in Macoma (Philippart et al., 2003;Beukema et al., 2009. A factor known to be associated with higher recruitment in Macoma is cold winter temperatures (Beukema et al., 1998). However, average mean winter temperatures in the years of higher recruitment (2011 ∼4.8 • C, 2012 ∼5 • C, 2013 ∼5 • C) were no different to the years prior to this (1997-2010 of ∼5 • C, Appendix 4 in Supplementary Material). Higher adult densities were also not the main factor driving higher recruitment, as there was no apparent relationship between the densities of adults and recruits. The lack of a single and simple association between recruitment success and adult density has previously been observed in other studies Andresen, 2013). A full understanding of the factors driving recruitment success in Macoma clearly needs more work. Changes in the living depth of Macoma over the study period could not be explained by changes in the densities of Macoma's main shorebird predator, red knots Calidris canutus. Although, red knot numbers have declined (van Roomen et al., 2005;Laursen et al., 2010;Wolff et al., 2010), the densities of red knots per unit suitable foraging area remained constant at 10 knots/ha of suitable foraging area from 1996 to 2005 in our study area . Flatfish are also known to affect living depth of Macoma via the foraging on siphons (de Vlas, 1981;de Goeij et al., 2001;Meyer and Byers, 2004). However, as there has been a dramatic decline in flatfish since the 1980s (van der Veer et al., 2011), and siphon nipping occurs mainly in spring (de Goeij, 2001), whereas our sampling was conducted over summer, siphon nipping is unlikely to be the cause for the living depth changes in Macoma. Increased parasitism is also unsupported as another cause for the change in living depth. Space and energy consuming macroparasites were only found in ∼10% of Macoma in the Wadden Sea in 2004 (H. Andresen unpub. data), a year when Macoma was in steep population decline, and also in 2014 (D. Thieltges pers. comm.) when the densities of Macoma were highest. Finally, changes in living depth cannot be explained by spatial or temporal (monthly) autocorrelation, as we have randomized our data such that this signal would disappear if these factors were contributing to the changes in living depth each year.
In larger Macoma (>14 mm) we observed an increase in body condition over time. As body mass prior to spawning is correlated with numbers of eggs produced, and larger Macoma individuals produce more eggs than smaller individuals (Honkoop and van der Meer, 1998;Beukema et al., 2009), our results are consistent with the idea that the reproductive potential of larger Macoma has improved in recent years. Interestingly, it has also been observed that larger adult Macoma suffered lower mortalities from 2006 to 2010 in the Balgzand area (Drent et al., 2008) compared to relatively higher mortalities from 1998 to ∼2005 (Beukema et al., 2009).
The significant correlation between depth of living and body condition in large Macoma, and our knowledge of the foodsafety trade-off in burrowing bivalves, suggests that the observed changes in living depth over time are associated with a change in food intake by Macoma. Macoma is a facultative deposit feeder that suspension and deposit feeds (Kamermans, 1994;Rossi et al., 2004), probably in relation to its depth of living (de Goeij and Luttikhuizen, 1998). Although, the measure of average phytoplankton in the water column (chlorophyll-a) has remained constant and variable from 1977 to 2009 (Philippart et al., 2003;, diatom or epiphytobenthos abundances on the sediment surface have not been measured over time. However, it is known that continual dredging at the same spot by fisheries disturbs the sediment surface, and thus destroys diatom mats, and increases turbidity and sediment resuspension (Mercado-Allen and Goldberg, 2011). Thus, as Macoma co-occurs with the fished cockle species in our study area (Kraan et al., 2007), and we here observed that Macoma's behavioral change overlapped closely with the end of mechanical cockle-dredging (in , van Gils et al., 2006, we suspect that the benthic diatom community was negatively influenced not only in the fished areas but also the neighboring areas during the time of cockle-dredging. The living depths of Macoma in year t and population growth between year t and year t+1 were correlated over time, such that a high per capita population growth rate was associated with deeper living, a correlation that was strongest in large Macoma. These results thus build on those of van  and showed that Macoma can now make the decision to live deeper at the cost of a lower food intake to protect its future reproductive potential. We propose that the depth of living (and the body condition) of Macoma have signaled a state change (sensu Scheffer et al., 2001) in the intertidal ecosystem of the Dutch Wadden Sea since the termination of mechanical cockledredging.

CONCLUSION
We here show that living depth in Macoma signals subsequent population changes in the western Dutch Wadden Sea, and is correlated with adult condition (enhanced reproductive output) and settlement success of the surviving juveniles. Consequently, we conclude that Macoma's burrowing behavior is an important signaling trait that is a first response when habitat conditions are altered and populations have yet to recover. We hypothesize that after the termination of cockle-dredging the intertidal ecosystem of the western Dutch Wadden Sea shifted back to a state amenable to Macoma. We propose that the messages carried by state-dependent behaviors such as burrowing depth in bivalves should receive more scrutiny as they signal ecosystem changes, whereas densities do not. We here suggest that ecosystem managers can use Macoma's burrowing behavior as an instrument to forecast ecosystem change.

AUTHOR CONTRIBUTIONS
TC and WB ran the statistical analyses and wrote the paper with TP and JvG. AD, SH, JtH, and NM were responsible for the data collection and laboratory analysis, whereas AK was responsible for data management. All authors contributed to the content of the paper.

ACKNOWLEDGMENTS
Sampling on such an extensive scale would have been impossible without the crew of the RV Navicula. We thank numerous volunteers and colleagues that have worked on this project. We acknowledge Casper Kraan for his involvement in the collection and organization of early data. The project Waddensleutels (funding by Waddenfonds, awarded to Han Olff and TP) made it possible for WB to contribute. JvG was funded by an NWO-VIDI grant (ALW 864.09.002) and TP and AD by the project Metawad (funding by Waddenfonds, WF 209925). Data collection started as part of the NWO-PIONIER project awarded to TP, with the continued sampling funded by NIOZ, NAM and the NWO ZKO grant. We thank Henk van der Veer, Gábor Zsolt Szabó and Micha Rijkenberg for valuable comments on previous drafts, and two anonymous reviewers for further useful comments.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fevo. 2016.00019