Introduction

Coastal environments are exposed to several anthropogenic disturbances, which can affect organisms at several biological scales, ranging from biogeochemistry and physiology up to the community level. Taxa that are sensitive to the disturbance can be killed or displaced, while the survival and proliferation of resistant taxa might be facilitated. As a result, there can be changes in populations and assemblages and a shift in the ecological processes shaping the communities with, often, a decrease of biodiversity (Underwood and Peterson 1988).

Marine beaches and other intertidal and subtidal habitats are the centres of outdoor leisure activities, which can have a strong impact on the ecology of these habitats. Tourist activities such as sun-bathing, collection of animals as baits or food and scuba-diving can disturb the ecosystem through a variety of ways, including their repeated trampling on the substratum (Keough and Quinn 1998, Davenport and Davenport 2006). Effects of trampling have been demonstrated to be detrimental for terrestrial plants and for organisms inhabiting sand-dunes (Liddle 1991; Sun and Walsh 1998 and references herein). Similarly, relatively recent research has demonstrated that trampling can be harmful also for seaweeds and seagrasses as well as for their associated fauna in the marine environment (Hawkins and Roberts 1993; Keough and Quinn 1998; Schiel and Taylor 1999; Brown and Taylor 1999; Eckrich and Holmquist 2000; Milazzo et al. 2004).

The effect of human trampling has been measured for populations and assemblages in marine rocky shores, coral reefs and recreational beaches, but less so for mudflats, which is partly explainable by the fact that tourist activities are lower on this type of habitats. On mudflats, the effect of trampling has been mainly considered as a side-effect due to the harvesting of baits or commercially important bivalves, which includes the disturbance due to animal removal and to the method of capture (bait-pumps or digging) (Peterson 1977; Wynberg and Branch 1997; Contessa and Bird 2004). However, there can be other relevant type of anthropogenic disturbance, in the form of trampling, such as leisure mud-walks that are now becoming a popular tourist attraction. Benthic ecologists can also be a further, unintended source of disturbance when they repeatedly walk on the mudflats to collect samples for monitoring programs or for extensive manipulative experiments. To our knowledge, the sole published work to have studied directly the effect of human footsteps on macrofauna was done on a footpath used by pilgrims visiting the Holy Island in the Lindisfarne National Natural Reserve (UK) (Chandrasekara and Frid 1996).

On hard substrata, trampling can affect marine organisms in a variety of ways: directly, by removing and crushing organisms or indirectly by displacing other species that interact through competition, predation or habitat provision (Brosnan and Crumrine 1994). For example, trampling often reduces cover of coral colonies and erected macroalgae directly, which has consequences on other species and assemblages because of habitat loss and niche reduction (Hawkins and Roberts 1993; Keough and Quinn 1998; Schiel and Taylor 1999; Brown and Taylor 1999; Eckrich and Holmquist 2000; Milazzo et al. 2004). Similarly, in seagrass meadows, physical action removes plants and rhizomes directly and, in turn, alters the habitat for the associated species (Eckrich and Holmquist 2000). Correspondingly, on mudflats, the mechanical disturbance of trampling can bury the animals living on the surface and bring deep burrowers to the surface, where they will die if not allowed to return quickly to depth. Footsteps can also disrupt the benthic biofilm and destroy animal burrows. This might vary the strength of biological interactions and have consequences on other organisms (Peterson 1977; Wynberg and Branch 1997; Contessa and Bird 2004 and references therein). On mudflats, trampling can also alter topographic complexity, which can indirectly affect recruitment and spatial distribution of microalgae (Wynberg and Branch 1994) and macrofauna (Rossi and Chapman 2003; Cruz-Motta et al. 2003). Furthermore, compaction of the sediment might alter the exchange of nutrients and oxygen between the sediment and the overlying water, change sediment accumulation rate and, again, modify population dynamics and distribution of animals in the mudflat (Contessa and Bird 2004).

This study aims at quantifying the effect of human trampling on biological and environmental variables in an intertidal mudflat, where marine ecologists had used a walking path during an extensive and prolonged field campaign in 2005. In particular, we focused this study on (1) the abundance and species composition of macrofauna, considering the population dynamics of the most abundant species of bivalves, on (2) benthic microalgal biomass and (3) sediment properties, including grain-size composition and concentration of the dissolved inorganic nitrogen (NH4 and NO3) in the pore water as an indicator for changes in the biogeochemistry of the sediment. Since trampling occurred on a single location and it was not possible to sample the area before the impact, we used an asymmetrical design, with data collected only after the disturbance (Glasby 1997). This approach falls within the so-called beyond BACI designs (Before-After-Control-Impact), which compares changes in single response variables at the putatively impacted location(s) to multiple controls, thereby taking into account natural patterns of variability and variability-induced by the source of disturbance under investigation (Underwood 1992, 1994). Recent developments have provided useful tools to extend the analytical framework of beyond-BACI designs into a multivariate context (Terlizzi et al. 2005). In this study we defined an impact as a difference between disturbed and control locations in the structure of macrofauna assemblages or in any other relevant ecological variable that may have been impacted by trampling.

Material and methods

Study area and sampling design

The study was done on the intertidal mudflat of Paulina Polder (51°21′ 23″ N, 3°42′ 49″ E; Westerschelde, The Netherlands) in autumn 2005.

The Paulina Polder is a typically Westerschelde mudflat—salt marsh ecosystem, characterised by an extended mudflat with a rich benthic community dominated by polychaetes (Heteromastus filiformis, Arenicola marina, Pygospio elegans) and molluscs (Macoma balthica, Cerastoderma edule, Hydrobia ulvae). The mudflat covers an area of approximately 1.0 km2 and has a mean tidal range of 3.9 m, with a semidiurnal regime. Its surface is mostly level because there are few minor drainage channels and major bedforms. The sediment is in general muddy with an average mud content (% < 63 μm) of 50%.

The extent of the study area was limited by a stone groyne on the western side and a creek on the eastern side. The impact area was a 25 × 25 m location in the centre of the mudflat (Fig. 1), where manipulative experiments had been run during 2005. Data relative to undisturbed sediment as sampled in November 2005 showed a macrofauna assemblage dominated by C. edule (mean ± SE: 2.5 ± 1.26 individuals/78 cm2), H. ulvae (mean ± SE: 8.25 ± 4.44 individuals/78 cm2), M. balthica (mean ± SE: 2.33 ± 0.76 individuals/78 cm2) and H. diversicolor (mean ± SE: 2.5 ± 0.87 individuals/78 cm2). There were also some juveniles of C. edule (1.25 ± 0.25 individuals/78 cm2) and of M balthica (3.75 ± 1.44 individuals/78 cm2) (Van Colen K, Rossi et al. unpublished data). At this location, walking paths were created among the experimental plots to allow movement and facilitate sampling. Since March 2005, several researchers (from 3 to 7 people) visited the location, first daily for a week and, then, weekly or fortnightly until the end of September 2005. Overall, five people of an average weight of 70 kg visited the location twice every month from March to September 2005 and walked on the footpath during low tide for a period between 3 and 5 h. The intensity of trampling was such that the entire walking path was covered by footsteps without any visible undisturbed areas. The average depth of the footsteps was between 3 and 5 cm and the disturbance was still visible at the next visit. This location will be referred as the disturbed, putatively impacted location (hereafter indicated as I). The two control locations (25 × 25 m, hereafter indicated as C1 and C2 or, globally Cs), were chosen at an equal shore elevation approximately 100 m to each side of the impact area to consider all the range of possible conditions in the area and to avoid any other disturbance caused by the creek and the stone groyne (Fig. 1).

Fig. 1
figure 1

Map of the study mudflat and position of the three sites in each of the three locations: the impacted location (I) and the two controls, C1 and C2

Full size image

Three sites of 4 × 4 m were randomly selected within the trampled area of I, in C1 and C2 and sampled with three replicates each. Logistical constrains prevented higher replication at the level of control locations, sites and units of observation. In addition, we sampled several response variables given the lack of studies in the literature providing indications on which response-variables could have been affected by the disturbance (thus indicating “an impact”) at the scales of our study.

The sampling started 4 days after the last visit of the researchers and it was repeated after 18 and 40 days to test for patterns of differences between I-v-Cs through time. Unfortunately, during the first campaign (“T1” 4th October 2005) it was impossible to completely sample C2 due to time constraints and the data relative to time 1 were not considered in the analyses. The other two campaigns were performed on 18th October and 8th November 2005 and are referred as T2 and T3 in the paper.

Data collection and analyses

At each site and time, three 50 × 50 cm plots were randomly chosen. In each plot we counted the number of footprints and of visible burrows of the lugworm Arenicola marina, we estimated the percentage of surface depressions and took photographs of the surface. In each plot, we took five measurements of the solar reflected upwelling radiance from the undisturbed sediment surface with a Ramses ARC diode portable spectrometer (Trios GmbH, Germany). Upwelling radiance from white polystyrene plates was then immediately measured, and the ratio of the two radiances was used to calculate the reflectance spectrum. The normalised difference vegetation index, NDVI based on a two wavelength algorithm ((750–675)/(750 + 675)) was calculated from reflectance values and used to estimate the surface chlorophyll a concentration (Kromkamp et al. 2006), which is a proxy for microalgal biomass. In each plot, one core (11 cm inner diameter, 30 cm depth) was taken to estimate the abundance of macrofauna. Two smaller cores (3 cm i.d.) were used to sample water content and grain-size at both 1 and 10 cm depth. The sediment from the two cores was, then, pooled to obtain one measure at each plot and depth. The analyses showed no differences between the top centimetre and the deep layers and only the results from the top centimetre are shown here for briefness. Finally, three additional cores (2 cm i.d.) were taken from the surface and pooled to estimate NH4 and NO3 in the pore water. All cores used represent standard sizes for macrofauna and sediment samples in most studies of intertidal sediments.

In the laboratory, samples for macrofauna were sieved at 1 mm, sorted alive, classified to species level and counted. Sieving at 1 mm and sorting alive excluded the collection of small polychaetes, which can, instead, be retained by a sieve of 0.5 mm. However, this procedure allowed us to process samples rapidly and the only common species not retained were Heteromastus filiformis and Pygospio elegans (Rossi et al., unpublished). The specimens of the clam Macoma balthica and the cockle Cerastoderma edule were subdivided into size-classes. The shell of the cockle was measured to the nearest millimetre and animals were separated into two size-classes (I: < 12 mm and II: > 12 mm) based on the bi-modal size-class distribution in the area (Rossi F, unpublished). The distribution was probably due to the fact that cockles were at the end of their recruiting and growing season. M. balthica was directly separated into size-classes according to its population dynamics, since individuals smaller than 10 mm are those that recruited during the year, while large individuals are older than 1 year. Very small animals (<5 mm) represent new recruits (Rossi et al. 2004). Thus, we considered three size-classes of M. balthica: I: < 5 mm (new recruits); II: between 5 and 10 mm (juveniles) and III: > 10 mm (adults).

Water content was measured by weighting sediment before and after desiccation by freeze-drying. The dried samples were resuspended in a buffer medium and shaken before grain size analyses using a Malvern Mastersizer 2000 laser diffraction analyser (Sperazza et al. 2004). Grain-size distribution were analysed in term of median grain-size in phi value (D50) and percentages of very fine sand (<125 μm) and silt and clay (<63 μm).

Porewater for DIN was extracted by centrifugation (2,000 rpm, 15 min) and the supernatant filtered (Millipore 45 μm). Samples were stored at −20°C and analysed within 3 weeks using automated colorimetric techniques.

Statistical analyses

Asymmetrical analyses of variance (Underwood 1992; Glasby 1997) were used to compare at I-v-Cs the following response variables: Arenicola marina burrows, porewater DIN and the abundance of the most abundant taxa. The model consisted of three factors: Location (L, one disturbed and two controls), Site (S(L) three levels, random, nested in location) and Time (T, 2 levels, random and crossed with L), with three replicates. For the analyses of NDVI, we included an additional source of variation, Plots (P, 5 levels nested in T × S(L) interaction).

For analyses, the Location was partitioned into two portions: I-v-Cs and the variability among Cs. The Time × Location (T × L) interaction was similarly divided into T × I-v-Cs and T × Cs portions. The terms S(L) and T × S(L) were similarly divided into S(Cs) and S(I), and T × S(Cs) and T × S(I). Analogously, when applicable, the term P (T × S(L)) was divided into the portion relative to the control locations, P(T × S (Cs)) and P(T × S(I)). Finally, the residual variation was divided into the residual variability for observations within Cs and within I (i.e. Res Cs and Res I).

Denominators for F ratios were identified following the logic of beyond-BACI design (Underwood 1992) and the modification given in Terlizzi et al. (2005). The homogeneity of variances was tested using Cochran’s C test and data were transformed when variances were heterogeneous. Transformed data are indicated in the results.

Permutational multivariate analyses of variance (PERMANOVA, Anderson 2001, 2005) based on Bray-Curtis dissimilarities on square-root transformed data was used to estimate variation in the multivariate assemblages between I-v-Cs. The analyses tested the same hypotheses described for univariate asymmetrical ANOVA but in a multivariate context (Terlizzi et al. 2005). Each term in the analyses was tested using 4999 random permutations of the appropriate units (Anderson and ter Braak 2003). Where the number of permutation units was not enough to get a reasonable test by permutation, a P value was obtained using a Monte-Carlo random sample for asymptotic permutation distribution (Anderson and Robinson 2003). The analyses were done using the computer programs DISTLM.exe and PERMANOVA.exe (Anderson 2004).

Multivariate patterns were visualised by non-metric multi-dimensional scaling (nMDS) ordinations of site’s centroids on the basis of a Bray-Curtis dissimilarity matrix on square-root transformed data. Centroids and distances among them in Bray-Curtis space were obtained using the computer program PCO.exe (Anderson 2003); nMDS plots of distance matrices among centroids were then generated with PRIMER 6 software (Clarke and Gorley 2001). Species responsible for the multivariate patterns were identified with SIMPER routine available in PRIMER.

Results

Abiotic variables

During the first campaign, when only C1 and I were sampled, the impacted location clearly showed footprints (average number ± SE: 4.6 ± 0.3; 3.3 ± 0.6; 3.6 ± 1.5 in site 1, 2 and 3 respectively), while the surface appearance was very irregular during the following campaigns and footprints less visible (Fig. 2a). The sediment surface showed regular ripples in the location C1, with an interval between crests of 8 cm. The surface of the location C2 was, instead, smooth or characterised by an irregular surface with small depressions of a few centimetre (Fig. 2b). On average, the percentage of depressions was higher at the impact than at the control locations during the second sampling date, but it was considerably lower at the final sampling (Fig. 3). This interaction was, however, not detected by the asymmetrical analyses of variance, which showed large variability among sites in all locations, variable with the time of sampling (T × S(L) interaction in Table 1).

Fig. 2
figure 2

Pictures of the sediment surface of a the impacted location, 4 days after the last visit of researchers (T1), after 21 and 40 days (T2 and T3, respectively); b control locations (C1 and C2)

Full size image
Fig. 3
figure 3

Mean (±SE, n = 9) percentage of depressions at each location. Data are averaged over three sites and three replicates per site. White and grey column indicate C1 and C2, respectively, while black columns refer to I

Full size image
Table 1 Summary of asymmetrical analyses of variance testing for differences between impact (I) and the average of controls (Cs) on abiotic variables
Full size table

There were no significant changes among locations or in their interaction with time in the grain-size composition of the top centimetre of sediment as summarised by the median grain-size measurements, D50 (Table 1; Fig. 4a). Rather, the control location C1 was characterised by a high fraction of very fine sand and low silt and clay compared to C2 and I (Fig. 4b, c). The analyses of variance done on these fractions were similar to the median grain-size and are not reported. The water content in the top centimetre also showed no changes in the impact location (Fig. 4d).

Fig. 4
figure 4

Mean (±SE, n = 9) sediment grain size (a–c), water content (d) and DIN in the top centimeter pore water (e–f). Data are averaged over three sites and three replicates per site. White and grey column indicate C1 and C2, respectively, while black columns refer to I

Full size image

The pore water concentration of DIN (NH4 and NO3) also showed no differences between I-v-Cs or its interactions with time (Table 1; NO3 not shown), although the concentration of NH4, but not NO3 slightly increased in the impacted location at the second time of sampling (Fig. 4e, f). Eventually, the values of NH4 were negatively correlated with the percentages of depression in the impact (r = −0.38, n = 18), but not in the control locations (r = 0.08, n = 36).

Macrofauna

The burrows of Arenicola marina were more abundant in the C1 than in both the C2 and I (Fig. 5a). There was however larger variability among sites in the impact location than in both the controls (Table 2).

Fig. 5
figure 5

Mean (±SE, n = 9) number of Arenicola marina burrows (a), and the abundance of the abundant taxa as indicated in the title of each graph (c–h). Data are averaged over three sites and three replicates per site. White and grey column indicate C1 and C2, respectively, while black columns refer to I

Full size image
Table 2 Asymmetrical analyses of variance testing for differences between impact (I) and the average of controls (Cs) on burrows of A. marina and abundances of H. diversicolo r and H. ulvae
Full size table

The abundances of mobile fauna (Hediste diversicolor and Hydrobia ulvae) at I were in the range of natural variability that characterised this mudflat (Fig. 5b, c) and there were no significant difference between I-v-Cs in their mean numbers or in the variability at the smallest scale and their interactions with time (Table 2).

There was evidence that the trampling changed the distribution of the bivalves and that the effect depended on their size-class. At I there were more newly-recruited individuals of Macoma balthica (<5 mm) than in both controls, which resulted in a significant difference between I-v-Cs (Table 3a, Fig. 5d). The abundance of the juveniles (5–10 mm) and 1-year old M. balthica ( > 10 mm) decreased in I at both sampling dates (Fig. 5e, f). However, the large variability among sites did not allow the analyses of variance to detect clear-cut patterns at P = 0.05 for the abundance of juveniles, whereas it showed significant differences in the mean numbers of adults between I-v-Cs (Table 3a). Cerastoderma edule class I (<12 mm) did not show any response to the trampling, but large variability among sites in all the locations (Table 3b, Fig. 5g). Cerastoderma edule class II (>12 mm) showed, instead, a sharp decrease in the impact compared to the control locations (I-v-Cs significant in Table 3b, Fig. 5h).

Table 3 Asymmetrical analyses of variance testing for differences between impact (I) and the average of controls (Cs) on Size-classes of (a) Macoma balthica and (b) Cerastoderma edule
Full size table

Multivariate analyses provided evidence for statistical differences in the composition of the assemblages at I-v-Cs. There was also large variability among control locations and an effect of the time of sampling (Table 4). This was also portrayed by the nMDS, which separated the centroids of impacted sites (Black symbols in Fig. 6) from those of the two control locations (white and grey symbols in Fig. 6). SIMPER analyses identified 6 taxa contributing the most to the differences between I and Cs (Table 5). This included not only the three size-classes of M. balthica and the adult C. edule, but also Hydrobia ulvae and H. diversicolor.

Table 4 PERMANOVA based on the Bray-Curtis dissimilarities of the multivariate data (square-root transformed, 15 taxa)
Full size table
Fig. 6
figure 6

Non-metric multi-dimensional scaling ordination (nMDS) based on the Bray-Curtis dissimilarity for the centroids of each site at each location and time. Squares: Time 2; Triangles: Time 3

Full size image
Table 5 SIMPER analyses to evaluate the contribution of the species to the differences observed between Cs and I
Full size table

Benthic microalgae

The values of microalgal biomass as estimated by NDVI, showed a decrease in the impact location compared to the control locations and a general increase in all the locations between T2 and T3 (Fig. 7). However, the asymmetrical analyses of variance did not detect any changes between I-v-Cs or its interaction with time, but differences between the impact and the controls at the lowest spatial scales (i.e. among plots). The Res I and the Res Cs were in fact not equivalent for this data, suggesting that further tests on the mean effects at the highest spatial scales and any interaction with time should be interpreted with caution. Furthermore, the values of microalgal biomass averaged per plot were highly and negatively correlated with the percentages of depression in the trampled plot (r = −0.75, n = 18), but not in the control locations (r = −0.07, n = 36). In addition, there was a positive correlation of microalgal biomass with the concentration of NH4 in the pore water (I: r = 0.53, n = 18; Cs: r = 0.27, n = 36) (Table 6).

Fig. 7
figure 7

Mean (±SE, n = 45) microphytobenthic biomass as estimated from the normalised difference vegetation index (NDVI). Data are averaged over three sites, three plots and five replicates. White and grey column indicate C1 and C2, respectively, while black columns refer to I

Full size image
Table 6 Asymmetrical analyses of variance testing for differences between impact (I) and the average of controls (Cs) on NDVI
Full size table

Discussion

Trampling disrupts the sediment surface and increases the complexity of the habitat which, in turn, can affect sediment dynamics and grain-size composition as well as several biogeochemical processes occurring in the sediment. The presence of overlying pools of water, for example, will lower the supply of atmospheric CO2 for photosynthesis by benthic microalgae during low tide, due to greatly decreased rate of diffusion (photosynthesis during high tide is low or zero in this estuary due to high turbidity of the water column). Reduced primary production can have important consequences on benthic productivity since microphytobenthos represent the main source of food for many infaunal species. In addition, when sediment is more compact, as in the depressions, the exchange of oxygen between the sediment and the water column as well as the penetration of water in the sediment can decrease considerably (Contessa and Bird 2004), which can regulate nitrogen cycling in the sediment.

Footprints were initially visible after the disturbance, but then faded and remained as wet depressions (Fig. 2). However, the percentage of substratum depressed in the impacted area did not differ consistently from the control locations. Rather, it showed variability among plots. Variability at small scales was also observed for grain-size and water content, NH4 in the pore water and for benthic microalgae, but, again, they did not show any clear response to the disturbance. However, with three replicates per site, and considering that biogeochemical processes can vary over very small scales, sample size was probably not enough to detect clear-cut patterns. More samples would be needed to increase the power of the analyses for these response variables. Therefore in this case, our failure to accept the alternative hypotheses should not be taken as an acceptance of the null hypotheses of no differences between the impacted and the control locations. Rather, we consider the results of the analyses uncertain with more resource to be invested in future investigations. There are, however, no published data on how trampling affects biogeochemistry of the sediment and our finding can be of great help for determining the sample size, the alternative hypotheses and the effect size needed to calculate the statistical power of future experiments. This could be especially true for benthic microalgae and NH4. It is worthy to notice that microalgal biomass was inversely proportional to the distribution of depressions in the impacted location, but not in the control locations. This suggested that footsteps might have changed the distribution of microalgae at some small spatial scales, and call for manipulative experiments to test these models. Furthermore, the concentration of ammonium in the pore water not only co-varied with the percentage of depressions, but also with the biomass of benthic microalgae. Microalgae are known to play an important role in the ammonium uptake and may regulate the denitrification/nitrification rate in the surface sediment (Risgaard-Petersen 2003). Although the descriptive nature of our study does not allow identifying causal processes, the correlations might imply an interaction between benthic microalgae and footsteps in regulating nitrogen cycling.

There was a very clear effect of the trampling on the macrofauna composition. Available accounts on the effect of trampling associated to the digging for bait collection report a general decrease in the numbers of tube-dwelling, sub-surface deposit-feeders and deep-burrowing species (Wynberg and Branch 1994, 1997 and reference herein). Chandrasekara and Frid (1996), on the contrary found that, except for two sub-surface deposit-feeders, overall animals did not respond and abundance of some animals even benefited from the trampling. Here, trampling did not affect mobile species such as Hydrobia ulvae and Hediste diversicolor, but it changed distribution and population dynamics of some bivalves (Macoma balthica and Cerastoderma edule). Trampling not only impacted negatively on adults, but it also increased recruitment of juvenile M. balthica. If we had considered the total abundance of those species, as in other studies (Chandrasekara and Frid 1996) we would have erroneously concluded that there was no impact or even a positive effect on macrofauna. Thus, it might be of interest for future studies on trampling to consider population dynamics of species, rather than only their abundance.

Trampling can damage directly the species, or produce a mixed response when other biological and environmental variables are changed (indirect effects, Keough and Quinn 1998). For instance, surface animals with not-well developed structure for burial can be killed by direct crashing or burial, which stop them migrating to the surface fast enough to avoid asphyxia (Chandrasekara and Frid 1997). The cockles, which normally inhabit the top 2–3 cm of the sediment, were probably unable to escape even small depths of burial during trampling. The sensitivity of cockles to burial was also shown during laboratory tests that measured an increase in mortality with increasing depth of burial and lethal effects at a depth of 10 cm (Jackson and James 1979; Chandrasekara and Frid 1997 and references therein). This, however, does not explain why we did not detect any changes for the smallest cockles, which live at the same depth into the sediment or even shallower. In this area, many benthic species, including both M. balthica and C. edule, have reproductive periods between April and October (Rossi et al. 2004), and during trampling larvae and juveniles would have been present in the water column to replace those displaced by the trampling. This might determine fast recovery of small cockles similarly to mobile animals (e.g. Hydrobia ulvae and Hediste diversicolor in our study), which recolonise trampled areas and overcome the effect of burial very rapidly (Chandrasekara and Frid 1997).

The collapse of burrows is considered another direct mechanism that kills animals living deep in the sediment and feeding on the surface, such as ghost shrimps (Peterson 1977; Wynberg and Branch 1994, 1997; Contessa and Bird 2004). Likewise, the clam Macoma balthica lives to a depth of 10 cm and extrude the siphon to feed at the surface-water interface (Rossi et al. 2004; Kamermans 1994). Thus, adult M. balthica were probably killed because footsteps destroyed their connection with the surface. In addition, the clam is an important prey for many aquatic birds and the vertical mixing of sediment could bring infauna to the sediment surface, exposing animals to predation (Chandrasekara and Frid 1996). The increased variability in Arenicola burrows might be also due to the disruption of part of their burrows by footsteps, but their heterogeneous distribution over the mudflat rendered difficult the detection of clear-cut patterns.

The response of juvenile M. balthica cannot be, however, explained if only direct impact (damage by burial) was considered, even including the possibility of fast recovery for juveniles. Adult cockles can out-compete larvae and spats of other bivalves and reduce their settling and recruitment by disturbing the sediment while feeding and moving or filtrating planktonic larvae (Andre and Rosenberg 1991; Kamermans et al. 1992). This was also observed in the present mudflat during another experiment, where there was a decrease in the numbers of M. balthica spats after adult cockles were transplanted into the sediment (Rossi et al. unpublished data). Additionally, it has been largely recognised that topographic heterogeneity can affect the recruitment rate. Although the percentage of depressions was not increased by trampling, these results were not conclusive (see above) and this explanation cannot be excluded. Small depressions, like those of footprints might have provided suitable micro-habitats for settlement and recruitment (Olafsson et al. 1994; Herman et al. 1999).

Although we did not measure the time of recovery from trampling, evidence suggests that after the direct damage, indirect effects can increase the opportunities for the populations of bivalves to recover, at least during the growing seasons. However, in this study only one species benefited from the disturbance and in the long-term the balance between Macoma balthica and Cerastoderma edule might shift towards a dominance of M. balthica in the macrofauna assemblage. As a consequence, there could be a change in ecosystem functioning because these species can play important, yet distinct roles for the mudflat. Deposit-feeding of Macoma balthica can modify the biogeochemistry of the sediments, while suspension feeding Cerastoderma edule can affect more the exchange of oxygen and nutrients at the boundary layer, but has fewer effects on the sediment (Widdows et al. 1998, 2000). Therefore, the effects of trampling should be considered seriously for intertidal mudflats. Long-term, well-designed monitoring plans and manipulative experiments will quantify the impact of this source of disturbance and its effects on biodiversity-ecosystem functioning relationship.