Royal Society Publishing


Despite the complexity of natural systems, heterogeneity caused by the fragmentation of habitats has seldom been considered when investigating ecosystem processes. Empirical approaches that have included the influence of heterogeneity tend to be biased towards terrestrial habitats; yet marine systems offer opportunities by virtue of their relative ease of manipulation, rapid response times and the well-understood effects of macrofauna on sediment processes. Here, the influence of heterogeneity on microphytobenthic production in synthetic estuarine assemblages is examined. Heterogeneity was created by enriching patches of sediment with detrital algae (Enteromorpha intestinalis) to provide a source of allochthonous organic matter. A gradient of species density for four numerically dominant intertidal macrofauna (Hediste diversicolor, Hydrobia ulvae, Corophium volutator, Macoma balthica) was constructed, and microphytobenthic biomass at the sediment surface was measured. Statistical analysis using generalized least squares regression indicated that heterogeneity within our system was a significant driving factor that interacted with macrofaunal density and species identity. Microphytobenthic biomass was highest in enriched patches, suggesting that nutrients were obtained locally from the sediment–water interface and not from the water column. Our findings demonstrate that organic enrichment can cause the development of heterogeneity which influences infaunal bioturbation and consequent nutrient generation, a driver of microphytobenthic production.


1. Introduction

It is rare in nature to find an entirely uniform habitat or for the distribution of organisms to be completely regular. Most organisms exhibit a patchy distribution reflecting the heterogeneous nature of the environment (Tilman et al. 1994; Williams et al. 2006). Therefore, it is surprising that the natural heterogeneity of ecosystems has rarely featured in the experimental analysis of ecosystem processes (Cheng et al. 2007; Holzschuh et al. 2007). Heterogeneity has functionally important consequences for productivity and other ecosystem services provided by an ecosystem, particularly if the transmission of material and resources between patches is slow or restricted (Strayer 2005). Heterogeneity is also known to be important in the maintenance of species diversity (Sommer 2000), habitat (Levinton & Kelaher 2004) and material and energy flow (Franklin 2005), such as nutrient cycling (Bengtson et al. 2006). It is clear that both local processes (Levinton & Kelaher 2004) and the landscape matrix, which they form, are important in determining habitat quality (Williams et al. 2006). If we are to fully understand the role of species in mediating ecosystem processes, particularly at larger scales, it is essential to integrate heterogeneity effects when considering overall habitat performance. Investigation of the spatial distribution of specific populations is common (Noren & Lindegarth 2005; Bengtson et al. 2006; Grenyer et al. 2006; Jones et al. 2006; Condeso & Meentemeyer 2007), and some studies now include links to related functional attributes. For example, Jesus et al. (2005) provided a detailed analysis of microphytobenthos (MPB) distribution on the surface of an estuarine mudflat and linked it to the photosynthetic functionality at a centimetre scale. However, the inclusion of spatial distribution patterns has not yet been incorporated as a treatment in studies of biodiversity and ecosystem functioning.

Coastal zones and estuarine ecosystems have proven to be valuable sites for the investigation of relationships between biodiversity and ecosystem function (BEF; for review, see Covich et al. 2004). Different attributes of the marine environment have been incorporated into experimental systems to test empirical relationships (e.g. flow, Biles et al. 2003; regional attributes, Emmerson et al. 2001; grazers, Duffy 2002, Hagerthey et al. 2002) using an approach (Raffaelli et al. 2003) that is analogous to those used in other systems (Schmid et al. 2002). In many instances, the rate or flux of nutrients has been used as a measure of ecosystem function (e.g. Ieno et al. 2006) and, for such point processes, spatial heterogeneity becomes important when considering nutrient cycling at larger scales.

In intertidal areas, one natural and reproducible element of heterogeneity is the patchiness of macroalgae (Hagerthey et al. 2003) and the associated physicochemical variability of the sediment bed (Raffaelli 2000). Buried algae decays rapidly providing resources for infaunal organisms (Rossi & Underwood 2002) but may also lead to sediment anoxia, thus the overall effect on organisms may be positive or negative. This may lead to opposing organizational forces (localized detrital input versus mobility of consumers) in deposit-feeding marine communities that exert structural control at the landscape scale (Levinton & Kelaher 2004). The major primary producers in mudflat systems are the MPB (Paterson & Hagerthey 2001), and it is known that their distribution can be patchy, varying over spatial scales of less than 1 cm (Jones et al. 2006), in response to environmental variables and macrofaunal composition (Christie et al. 2000; Hagerthey et al. 2002). The biomass of MPB can be assessed by non-destructive pulse-amplitude modulated (PAM) fluorescence techniques (Honeywill et al. 2002; Consalvey et al. 2004a; Jesus et al. 2006b), which allow repeated measurements over restricted spatial and temporal scales (Jesus et al. 2005).

Here we manipulated the spatial heterogeneity within mesocosm systems by the selective addition of detrital algal material to a defined region of sediment. The influence of this induced heterogeneity on ecosystem function was assessed using MPB biomass distribution as a proxy for photosynthetic capacity (Honeywill et al. 2002; Consalvey et al. 2004b; Jesus et al. 2006a). The factorial experiment was designed to examine the influence of species identity, species density and algal enrichment (as a mechanism for inducing heterogeneity) on microphytobenthic primary production. We hypothesized that (i) macrofaunal distribution (identity and biomass) would influence production capacity, but (ii) this would be influenced by the patchiness created in the experimental system.

2. Material and methods

Sediment was collected from the Ythan estuary, Aberdeenshire, Scotland, UK (57°20.085′ N, 02°0.206′ W) and sieved (500 μm) in seawater to remove unwanted macrofauna, and left to settle for 24 hours to retain the fine fraction (less than 63 μm). Excess water was removed, the sediment slurry homogenized and distributed between mesocosms (opaque aquaria, 21×15×14 cm). Sediment was added to each mesocosm to a depth of 3 cm. Enrichment was achieved by the addition of dried and ground Enteromorpha intestinalis collected from the Ythan Estuary. Perspex sheets were used to divide the mesocosms into equal halves and 1g of algae was added to enrich selected patches (equivalent to 126 g m−2, within levels found naturally; Raffaelli 2000). Mesocosms were initially filled with 2.5 l seawater (UV-sterilized, 10 μm pre-filtered, salinity≈33), left for 24 hours and refilled with seawater to eliminate nutrient pulses associated with assembly (Ieno et al. 2006). Mesocosms were placed in a controlled temperature room (11±2°C), aerated and the photoperiod was set to a 12 hour light–dark cycle (26 mm Ø white fluorescent tubes, model GE F36W/35; 36W, 3500 K). The experiment ran for 10 days.

To determine the effects of macrofaunal species identity, macrofaunal species biomass and algal enrichment on MPB biomass, 396 mesocosms were established, divided randomly and equally between two experimental runs. Two patches were established in each mesocosm. The deposit-feeder Hediste diversicolor (Polychaeta), the surficial grazer Hydrobia ulvae (Gastropoda), the regenerator Corophium volutator (Crustacea) and the suspension/deposit-feeding bivalve Macoma balthica were added on day 0. Macrofauna were confined to their initial patches for 24 hours using Perspex dividers. Combinations of macrofaunal biomass (0, 25, 50 and 100% of natural density in both the left and right patches, i.e. 16 possible combinations) were established for all possible interface combinations of patch arrangements (E|E, E|NE, NE|E and NE|NE where ‘|’ represents the interface and E=enriched and NE=non-enriched; the measured patch is on the left of | and a neighbouring patch is on the right of |) for each of the four macrofaunal species (figure 1). Each configuration was replicated three times (n=396). For M. balthica and H. diversicolor, whole individuals were counted and four individuals per patch were taken as analogous to the natural density on the Ythan estuary. For H. ulvae and C. volutator, the natural wet weight biomass was determined (2 and 1 g per patch, respectively) and appropriate proportional wet weights added to the mesocosms. In addition, replicate (n=3) control mesocosms containing no macrofauna were established for each interface configuration (n=12) to determine the effect of the presence of macrofauna, irrespective of identity.

Figure 1

Overview of experimental design. The species density gradients across the patch interface were established at the start of the experiment, using the relative levels of 0, 25, 50 and 100% natural density at the study site. These combinations were used for each of the four interface treatments (enriched is shaded and non-enriched is not shaded), and every species density–interface combination was used for each of the four species (Cv, Corophium volutator; Hu, Hydrobia ulvae; Mb, Macoma balthica; Hd, Hediste diversicolor).

An emergent property of the experimental design allowed the analysis of the influence of the difference between the initial and final biomass of the macrofauna set in adjacent patches. The initial density difference was expressed numerically as the difference in biomass between the measured and the adjacent patches, such that: +4=all macrofauna (at maximum biomass) were in the measurement patch; 0=equal distribution in each patch; and −4=all macrofauna (at maximum biomass) were in the adjacent (non-measured) patch.

Measurements of MPB biomass were taken on day 6 based on PAM fluorescence (Consalvey et al. 2004a) using a Hansatech FMS2 meter. A 6-day interval was appropriate because this was the combination for optimum MPB biomass and species activity (Defew et al. 2002). Mesocosms were dark adapted for 15 min to optimize MPB biomass estimates from the Fo15 output, which is a proxy for Chl a content (Honeywill et al. 2002). To reduce variability, two measurements of Fo15 were taken from each patch and averaged, and three replicate mesocosms were measured for each treatment (n=3).

A generalized least squares (GLS; Pinheiro & Bates 2001) statistical mixed modelling approach was used to assess the two experimental hypotheses. A GLS framework is preferential over linear regression using transformed data because it retains the structure of the data while accounting for unequal variance in the variance–covariate terms. As a first step, a linear regression model was fitted. Model validation was applied to verify that underlying statistical assumptions were not violated; normality of residuals was assessed by plotting theoretical quantiles versus standardized residuals (QQ plots), homogeneity of variance was evaluated by plotting residuals versus fitted values, and influential data points were identified using Cook's distance method (Quinn & Keough 2002). The validation procedure showed that there was no evidence of nonlinearity but there was evidence of unequal variance among the explanatory variables. The GLS model was refined by manual backwards stepwise selection using maximum likelihood (ML) to remove insignificant terms, and the final model was presented using restricted maximum likelihood (REML; West et al. 2007). The highest potential level of interaction that was assessed was the three-way interaction terms. The statistical outputs of these models are based on the comparisons of the first level within each term with all other levels; no other within-level comparison is made. To assess the importance of individual independent variables, a likelihood ratio test was used to compare the full minimal adequate model with models in which the independent variable, and all the interaction terms it was involved with, were omitted. As a complementary indicator of the importance of these individual variables, in each case we calculated the decrease in the adjusted R2 value for the model without that variable as compared with the full model. Analyses were performed using the ‘R’ statistical and programming environment (R Development Core Team 2005) and the ‘nlme’ package (Linear and nonlinear mixed effects models; Pinheiro et al. 2006).

3. Results

The minimal adequate model was a linear regression with a GLS extension incorporating four two-way interaction terms and four single terms (adjusted R-squared=0.49). Single factors were species identity, interface, species density and initial density difference. Two-way interactions were species identity×interface, species identity×species density, interface×species density and species density×initial density difference. There were no significant three-way interactions. The variance–covariate terms were species identity, interface and species density.

(a) Independent terms

Species identity had the greatest influence on MPB biomass (L-ratio=755.18, d.f.=15, p<0.0001, decrease in adjusted R-squared Embedded Image), followed by interface type (L-ratio=425.78, d.f.=15, p<0.0001, Embedded Image), species density (L-ratio=218.34, d.f.=33, p<0.0001, Embedded Image) and initial density difference (L-ratio=9.78, d.f.=2, p<0.005, Embedded Image). As the source of nutrients fuelling MPB growth can either originate from bottom-up (sediment) or top-down (water column) processes, we compared MPB biomass in non-macrofaunal control mesocosms. These analyses showed that the focus patches (left) had a significant effect on MPB biomass, while neighbouring patches (right) had no significant effect (two-way ANOVA: left patch, F=5.93, d.f.=1, p=<0.05; right patch, F=0.26, d.f.=1, p=0.627), indicating that bottom-up processes were determining MPB biomass.

(b) Two-way interaction terms

The significant two-way interaction terms, in order of importance, were species identity×interface (L-ratio=39.18, d.f.=33, p<0.0001; figure 2), species identity×species density (L-ratio=38.13, d.f.=39, p<0.0001; figure 3), species density×interface (L-ratio=24.15, d.f.=39, p<0.0001; figure 4) and species density×initial density difference (L-ratio=4.42, d.f.=41, p=0.036; figure 5). Hediste diversicolor had the weakest effect in terms of MPB response (highest MPB biomass) followed by M. balthica, H. ulvae and C. volutator (lowest MPB biomass). There was a consistent pattern in that enriched patches (E) maintained higher levels of MPB biomass than non-enriched (NE) patches (figure 2). However, while the fully enriched condition (E|E) maintained the highest biomass of MPB for H. diversicolor and H. ulvae, this was not the case for M. balthica or C. volutator, where the heterogeneous condition (E|NE) maintained the highest level of biomass. Within the interaction species identity×interface, M. balthica×E|NE (p=0.019, coefficient (95% CI)=70.26 (11.54–129.00)), C. volutator×E|NE (p=0.017, coefficient (95% CI)=56.11 (10.29–101.94)) and C. volutator×NE|E (p=0.027, coefficient (95% CI)=43.80 (4.93–82.66)) were significant compared with H. diversicolor×E|E. The nature of the interaction was to increase the influence of the E|NE condition (figure 2), so that for these two species, the interface condition positively influenced MPB biomass.

Figure 2

Graphical representation of the effect of the two-way interaction term species identity×interface. Vertical lines represent species identity: Hd, Hediste diversicolor; Mb, Macoma balthica; Hu, Hydrobia ulvae; Cv, Corophium volutator. Horizontal bars represent predicted values from the optimal regression model for each heterogeneity treatment, focal ‘patches’ are represented by the expression on the left of ‘|’ while neighbouring patches are on the right. The four horizontal lines are the average for control mesocosms (containing no macrofauna) at interface treatment E|E (solid line), E|NE (dashed line), NE|E (dot-dashed line) and NE|NE (dotted line). As the GLS framework allows for different spreads in the data, individual data points are omitted to prevent misinterpretation.

Figure 3

Graphical representation of the effect of the two-way interaction term species identity×species density. Lines represent species identity: H. diversicolor (solid line); M. balthica (dashed line); H. ulvae (dotted line); and C. volutator (dot-dashed line). Species density is expressed as a percentage of the natural densities at the study site. As the GLS framework allows for different spreads in the data, individual data points are omitted to prevent misinterpretation.

Figure 4

Graphical representation of the effect of the two-way interaction term interface×species density. Lines represent heterogeneity: E|E (solid line); E|NE (dashed line); NE|NE (dotted line); and NE|E (dot-dashed line), where E is an enriched patch, NE is a non-enriched patch and ‘|’ is the interface between each patch. Analysis is based on the left patch and coded for neighbouring patch on the right. Species density is a percentage of the natural densities found at the study site. As the GLS framework allows for different spreads in the data, individual data points are omitted to prevent misinterpretation.

Figure 5

Graphical representation of the effect of the two-way interaction term species density × initial density difference. Lines represent the initial density difference: −4 (solid line); 0 (dashed line); 4 (dotted line), where initial density difference ranges from a maximum density in the right-hand patch and no macrofauna in the left-hand patch (−4) to a maximum density in the left hand patch and no macrofauna in the right-hand patch (4). Species density is a percentage of the natural densities found at the study site. As the GLS framework allows for different spreads in the data, individual data points are omitted to prevent misinterpretation.

For each species, there was an overall reduction in MPB biomass with increasing density (figure 3). At low-density levels, H. diversicolor had least effect (highest MPB biomass), M. balthica and H. ulvae had similar effects and C. volutator had the greatest effect on MPB biomass (lowest MPB biomass). As species density increased, the rate of decline in MPB biomass was similar for all species with the exception of M. balthica, for which it was less pronounced (p=0.0033, coefficient (95% CI)=11.48 (3.82–19.15)).

The interaction species density×interface showed an overall reduction in MPB biomass as the density of each species increased in all treatments except for NE|NE (figure 4). The rate of change in MPB biomass was similar between the three declining treatments. At low densities, MPB biomass varied with interface treatment, with the highest biomass associated with E|NE followed by E|E, NE|E and NE|NE.

The interaction species density×initial density difference was also significant, but very weak (L-ratio=4.42, d.f.=41, p=0.036). Model visualization (figure 5) indicates that the level of MPB biomass declined as species density increased. The rate of decline was greatest in mesocosms with the maximum biomass in the focal patch and zero biomass in the neighbouring patch, followed by treatments with initial densities evenly distributed between patches, and mesocosms with the maximum biomass in the neighbouring patch and zero biomass in the focal patch.

4. Discussion

The mesocosm experiments were designed to examine the influence of spatial heterogeneity on MPB production. It was established that the addition of powdered algae, E. intestinalis, provided a suitable mechanism to induce system heterogeneity. The highest MPB biomass was recorded from enriched (E|E) mesocosms and the lowest in mesocosms that had not been enriched (NE|NE). The macrofaunal species selected for this work were known to be consumers of diatoms; therefore, there was an a priori reason to expect an effect on MPB biomass, both through the direct effects of grazing and through the indirect influence of nutrient release through sediment bioturbation. However, the macrofauna used have different bioturbatory characteristics, and these are variable depending on environmental conditions (Biles et al. 2003). Effects of species density on MPB biomass will therefore reflect the behaviour of individual species.

To date, BEF effort has included studies on species identity, diversity, biomass and functionality but without reference to the inherent natural variability of habitat. While the impact of spatial heterogeneity on ecosystem function has been considered (Lovett et al. 2005), empirical data are largely lacking. This contribution represents an initial empirical step to consider the role of spatial heterogeneity. It should be noted, however, that the classical mesocosm approach to measure nutrients as a proxy for ecosystem functioning (Raffaelli et al. 2003) when considering spatial heterogeneity is of limited value in marine benthic systems since the measured effects are integrated at the water column level, and a local contribution cannot be ascertained. The localization capability of PAM fluorescence allows the variation of MPB biomass between patches to be measured conveniently at a range of spatial scales.

The current experimental approach was firmly based on previous studies (Emmerson et al. 2001; Raffaelli et al. 2003; Ieno et al. 2006), but we acknowledge that the construction of synthetic assemblages in a mesocosm system is open to criticism (Carpenter 1996). It is important to reiterate here that, despite the apparent limitations of mesocosm systems, they allow the theory to be tested and global-scale environmental problems to become amenable to experimental endeavour (Benton et al. in press). Such systems are not an accurate representation of natural systems, rather they allow the development of theory that can be later tested under more natural conditions as the theory and practice of BEF research develop (Loreau et al. 2001; Naeem & Wright 2003; Balvanera et al. 2006; Raffaelli 2006).

The statistical model indicated that species identity, type of interface (heterogeneity) and species density were the strongest determinants of ecosystem response. The influence of species identity and density is unsurprising and consistent with numerous studies (for review, see Covich et al. 2004). Of particular significance, however, is that the macrofaunal species used in this study represent varied functional attributes and have clear trophic connections with the response variable, yet heterogeneity (interface type) was a driver for two of the three strongest interaction terms in the model. It is clear that spatial heterogeneity is of absolute importance and that point measurements of function may lead to qualitatively different and scale-dependent interpretations that are not relevant when considering processes at an ecosystem scale.

Decaying macroalgae within an estuarine system can increase the organic content and, under the right conditions, increase nutrient levels within the sediment in the immediate vicinity (Raffaelli 2000). MPB can also use these resources during photosynthesis to enhance production and levels of biomass. Localized enrichment has also been shown to influence macrofaunal behaviour (Levinton & Kelaher 2004). Previous work found that H. diversicolor and H. ulvae move towards enrichment, whereas C. volutator moves away (Lawrie et al. 2000; Raffaelli 2000), and M. balthica shows very little movement (De Goeij & Luttikhuizen 1998). Here, H. diversicolor had a positive effect on MPB biomass compared with the other species, although this did decrease with increasing biomass. This positive effect was possibly due to its relatively large size and bioirrigatory capacity (Magni & Montani 2006) increasing nutrient turnover, as well as stimulating microbial activity (Hansen & Kristensen 1997). In contrast, H. ulvae, although highly active, had limited impact on sediment nutrient turnover (consistent with Ieno et al. 2006; Orvain 2006) while the behaviour of M. balthica, which tends to siphon-feed in still water (Kamermans 1994; De Goeij & Luttikhuizen 1998), is unlikely to impact on MPB biomass. Although C. volutator has been shown to be highly active and mediate the release of comparatively large quantities of nutrient (NH4–N; Emmerson et al. 2001), the low MPB biomass levels found in C. volutator treatments appear to be influenced by a secondary effect caused by the behaviour of this species. Sediment resuspension during burrow maintenance causes the water column to become turbid, attenuating light and reducing photosynthesis at the sediment surface (De Deckere et al. 2000).

MPB can use nutrients generated in the enriched sediments at the sediment–water interface or from nutrients previously released into the water column. Nutrients released into the water column become available for the whole mesocosm and any response is likely to be effected over the entire system. It follows, therefore, that if the MPB obtains nutrients locally from the sediment–water interface, any observed responses in our experimental system would only occur in algal-enriched sediment. Overall, the highest MPB biomass was in enriched patches and the lowest MPB biomass was in non-enriched patches, irrespective of the neighbouring patch type. It is clear, therefore, that the source of nutrients for MPB is derived locally from the sediment–water interface rather than the water column itself and that sediment heterogeneity is an important determinant of MPB production.

Heterogeneity was induced by the addition of allochthonous carbon that may have both direct and indirect effects on the functional response of the system. The principal direct effect was expected through the release of resources (nutrients) that enhance MPB biomass at the sediment surface. In addition, the presence of organic material will influence the behaviour and migration of the macrofauna (Raffaelli et al. 1991; Rossi 2006) with a consequent feedback on MPB biomass (Hagerthey et al. 2002). This feedback is difficult to predict, as the effect may be positive (bioturbation releasing nutrients) or negative (grazing of MPB). Our results suggest that the important independent variables for MPB, in order of greatest effect, are macrofaunal species identity, the nature of the interface between two patches, macrofaunal density and the gradient in macrofaunal biomass between two patches. Although the interactions between these factors were more complex, the influence of system heterogeneity is clearly a significant factor for MPB performance, particularly in the case of C. volutator and M. balthica. When these species are present, the statistical model indicated that the functionality was higher than expected, suggesting that any negative effect of the species (direct grazing) was more than compensated for by the positive effects of bioturbation, such as increased nutrient turnover. This point is not trivial, as it has important ecological consequences since growth may be enhanced sufficiently to compensate for grazing pressure and result in increased standing stock (production). This suggests that the landscape matrix is more important than local ecosystem structure in determining MPB production (Williams et al. 2006) and may, in the longer term, have consequences for macrofaunal fitness and reproductive capacity. The model does not allow for more specific determination of interaction terms (suitable post hoc analyses are not possible), but it does highlight the overall importance of the interface. Elucidating the mechanistic effect requires further work but is likely to be a combination of species movement expressed through bioturbation, grazing and nutrient recycling.

5. Conclusions

Spatial heterogeneity plays an important role in determining MPB production, interacting with both macrofaunal species identity and density, even at the restricted level of patches within our experimental mesocosms. In nature, these effects are likely to be widespread. Attention must now be given to the development of novel methodologies capable of incorporating these interactions, to further elucidate the nature of the relationship between habitat heterogeneity and ecosystem function and the mechanisms underlying them, as well as the consequences for the conservation of biodiversity and ecosystem services in changing environments.


This work was supported by NERC studentship NER/S/J/2003/12648, tied to NERC research grant NER/A/S/2003/00577. D.M.P. and K.E.D. acknowledge the support by the MarBEF Network of Excellence ‘Marine Biodiversity and Ecosystem Functioning’ which is funded by the Sustainable Development, Global Change and Ecosystems Programme of the European Community's Sixth Framework Programme (contract no. GOCE-CT-2003-505446). This publication is contribution number MPS-07044 of MarBEF. The authors thank Craig Dearing, Anne Holford, Patricia Lastra, Owen McPherson and Leigh Murray for assistance in the field, and Highland Statistics Ltd for statistical advice. We are grateful for the comments of two anonymous referees, which greatly improved the manuscript.


    • Received July 9, 2007.
    • Accepted July 16, 2007.


View Abstract