Much of what we know about the role of biodiversity in mediating ecosystem processes and function stems from manipulative experiments, which have largely been performed in isolated, homogeneous environments that do not incorporate habitat structure or allow natural community dynamics to develop. Here, we use a range of habitat configurations in a model marine benthic system to investigate the effects of species composition, resource heterogeneity and patch connectivity on ecosystem properties at both the patch (bioturbation intensity) and multi-patch (nutrient concentration) scale. We show that allowing fauna to move and preferentially select patches alters local species composition and density distributions, which has negative effects on ecosystem processes (bioturbation intensity) at the patch scale, but overall positive effects on ecosystem functioning (nutrient concentration) at the multi-patch scale. Our findings provide important evidence that community dynamics alter in response to localized resource heterogeneity and that these small-scale variations in habitat structure influence species contributions to ecosystem properties at larger scales. We conclude that habitat complexity forms an important buffer against disturbance and that contemporary estimates of the level of biodiversity required for maintaining future multi-functional systems may need to be revised.
Natural ecosystems are open systems composed of interconnected gradients, patches and networks, in which matter, energy, information and organisms are continually exchanged [1–4], yet human activities and natural disturbances continually modify these inter-connections by increasing the number of patches, decreasing individual patch sizes and altering their relative isolation . Such changes in habitat configuration and complexity have profound effects on biodiversity, altering the relative abundance (evenness) and number of species over both ecological and evolutionary timescales [5–7], with the potential to mediate present and future levels of ecosystem functioning [8–10].
Despite the widely acknowledged importance of dispersal and movement between patches for the maintenance of diversity [11,12], and the recognition that habitat heterogeneity and the structure of biotic communities (abundance, evenness, distribution and diversity) are inextricably linked [13–15], the explicit incorporation of spatial components as structuring features in biodiversity–ecosystem function experiments are rare. Where components of habitat structure have been considered, species that were previously found to be of minor importance in mediating ecosystem properties in homogeneous experimental systems (e.g. [16,17]) have subsequently been shown to be of much greater importance in equivalent heterogeneous systems . Moreover, the way species interact with their environment has the potential to alter the nature of trophic interactions  and the fundamental form of the relationship between biodiversity and ecosystem properties [20,21], rekindling debate over the relevance and capacity of manipulative experiments to inform management, conservation and policy [22–26].
Theory (e.g. [12,20]) and empirical (e.g. [7,27]) evidence suggest that the dispersal of individuals between habitat patches can enhance biodiversity effects on ecosystem functioning and community stability because environmental heterogeneity promotes niche diversity [28,29]. Dispersal between resource patches, however, may also lessen variability in grazer abundance and reduce diversity effects on primary production . Nevertheless, it is clear that spatial processes have the potential to generate nonlinear effects . Thus, the functional consequences of altered biodiversity at scales larger than single habitat patches may be very different from predictions based on contemporary views formed from experimentation or small-scale observations [30–32].
Here, we report the results of an experiment in which a multi-patch mesocosm approach was adopted to determine the effects, alone and in combination, of benthic macrofaunal composition (three single-species treatments and one multi-species treatment) and habitat structure (spatial heterogeneity, homogeneous versus heterogeneous resource distribution; patch connectivity, 0, 50 and 100 per cent connectivity between patches; and patch location, seven locations per mesocosm) on ecosystem properties at two (patch and multi-patch) spatial scales. At the patch scale we quantify the intensity of infaunal sediment mixing (bioturbation) and the net movement of invertebrates between patches, while at the multi-patch scale we quantify the effects of species composition and habitat structure on nutrient concentration. Bioturbation varies with environmental context as well as benthic community structure [33,34], and affects carbon processing, nutrient cycling and oxygen dynamics , thus constituting the mechanistic link between the environment, biodiversity and ecosystem functioning in the marine benthos [3,36]. Here, we test the hypothesis that ecosystem processes and functioning are positively affected by species composition, but that the level of ecosystem functioning is mediated by habitat structure.
2. Material and methods
Two types of sediment (dune sand and estuarine mud), the green algae Ulva intestinalis (Chlorophyceae) and three macrofaunal invertebrates (the polychaete Hediste diversicolor, the amphipod Corophium volutator and the gastropod Hydrobia ulvae) were collected from the Ythan Estuary, Scotland (57°20.085′ N, 02°0.206′ W). Mud was sieved (500 µm mesh) in a sea water bath to remove macrofauna, allowed to settle for 24 h to retain the fine fraction (less than 63 µm) and homogenized. Sand was dry-sieved (500 µm mesh) to remove debris and plant material.
Multi-patch mesocosms (n = 150) consisted of transparent Perspex tanks (internal measurements, length 66 × width 13 × height 32 cm) containing four patches that were either unconnected (0 per cent connectivity) or interconnected by dispersal corridors of constant length (9 cm) but varying width (50 per cent connectivity, width = 6.5 cm; or 100 per cent connectivity, width = 13 cm), giving a total of seven locations per mesocosm (electronic supplementary material, figure S1). These dimensions are consistent with the observed levels of patchiness at the study site (S. M. Lawrie 1996, unpublished PhD thesis). In the 0 per cent connectivity treatments, the patches and corridors were separated by a transparent Perspex barrier (width 13 × height 32, wall thickness = 0.5 cm) to prevent invertebrate relocation. In the 50 per cent connectivity treatments, corridor widths were reduced to half the width (6.5 cm) of the mesocosm using two U-shaped Perspex barriers (each length 10 × width 3.25 × height 32 cm; electronic supplementary material, figure S1). Each barrier incorporated a mesh- (355 µm) covered hole (diameter = 4 cm) located at 10 cm water depth to allow transfer of water and nutrients, but not invertebrates, between patches. In the 100 per cent connectivity treatments, there were no such barriers. The patches contained sieved mud (approx. 1.2 l per patch, depth = 10 cm) and the corridors contained a 1 : 1 mixture of mud and sand (approx. 1.2 l per corridor, depth = 10 cm) to create a less desirable, but not inhospitable, habitat for the faunal species. To ensure that all mesocosms contained an identical volume of sediment, the U-shaped barrier lumen in the 50 per cent connectivity treatments also contained the sand–mud mixture. Each mesocosm was filled with sea water (18.5 l, UV-sterilized, 10 µm pre-filtered, salinity ≈ 33), held in a constant-temperature room (11 ± 2.0°C) with a 12 h L : 12 h D regime (2 × 26 mm diameter white fluorescent tubes, model GE F36W/35; 36 W, 3500°K) and continually aerated (one airline per patch). Sediment and sea water were added to each mesocosm and sea water was replaced after 24 h to remove excess nutrients associated with the assembly.
To incorporate habitat heterogeneity, patches were either non-enriched (NE) or enriched (E) with dried, ground U. intestinalis [18,19]. We assembled two heterogeneity configurations (HA and HB). In HA, locations 1 and 5 were enriched with algae (i.e. a patch sequence of E|NE|E|NE, where ‘|’ represents a corridor), while in HB locations 1 and 3 were enriched with algae (i.e. E|E|NE|NE). HA thus represented a more heterogeneous distribution of resources relative to HB.
Replicate (n = 5) faunal communities (n = 120) were assembled in monoculture and in three species mixtures. Mesocosms without macrofauna were also included (n = 30). In each mesocosm, the total macrofaunal biomass was fixed at 4 g per mesocosm (initially distributed as 1 g per patch and 0 g per corridor), similar to levels found at the study site (e.g. ). In mesocosms with 50 and 100 per cent connectivity, macrofauna were confined to their initial patches for the first 24 h using Perspex dividers. To visualize particle movement and quantify bioturbation for each patch, 2 g dry weight of luminophore tracers (dyed natural sediment that fluoresces under ultraviolet light ; 125–250 µm diameter; orange colour; Partrac Ltd, UK) were added to each patch 24 h after introducing the macrofauna. Owing to the large number of combinations involved (n = 150), mesocosms were randomized across five experimental runs (n = 30 mesocosms per run) of 10 days duration (sufficient time to allow faunal-mediated nutrient generation, while avoiding vertical homogenization of luminophore tracers through bioturbation; e.g. ).
At the end of each experimental run, a sediment core (diameter = 3.6 cm, 10 cm deep) was taken from the centre of each patch and vertically sliced in 0.5 cm intervals to a depth of 2 cm and in 1 cm intervals from 2 cm to a depth of 10 cm (i.e. 600 cores at 12 slices per core = 7200 slices). Each sediment slice was homogenized, evenly spread on a Petri-dish (diameter = 90 mm) and dried for 24 h at 55°C and 10 per cent humidity in an environmental chamber (VC 4100, Vötsch Industrietechnik). Images of each slice were taken (Canon EOS 400D; image size 3888 × 2592 pixels, effective resolution 11.97 µm per pixel) in an imaging box illuminated by an ultraviolet light source (2 × 8 W tubes, Sylvania Blacklight). To maximize the distinction between the luminophore particles and the sediment, a lighting filter (medium yellow no. 10, Lee Filters, UK) was mounted in front of the camera lens. Following Solan et al. , images were saved in RGB colour mode with JPEG compression and analysed using a custom-made semi-automated macro that runs within ImageJ (v. 1.40), a Java-based public domain program (available at http://rsb.info.nih.gov/ij/index.html). Bioturbation intensity (biodiffusion coefficient, Db) was determined for each macrofaunal species using the one-dimensional diffusion model by Crank  and model fit (r) was estimated following Francois et al. .
Net movement of macrofauna between patches was estimated as the change in biomass (Δg) within each patch between the start and end of the experiment. Pre-filtered (0.45 µm, Nalgene) water samples were taken from approximately 5 cm water depth at location 4 (centre; see electronic supplementary material, figure S1) of each mesocosm. NH4-N and PO4-P concentrations were determined using standard protocols with a Tecator flow injection auto-analyser (FIA Star 5010 series) using an artificial sea water carrier solution.
(a) Statistical analysis
(i) Patch-scale models
Statistical models were developed for each of the dependent variables, net movement and bioturbation. The models incorporated species composition (five levels, including a no-macrofauna treatment), location (net movement, four habitat patches and three corridors; bioturbation, four habitat patches), heterogeneity (two levels) and patch connectivity (net movement model, two levels; bioturbation model, three levels) as nominal explanatory variables. For the net movement model, the 0 per cent connectivity treatment (constrained movement) and the no-macrofauna treatment were removed from the analysis.
(ii) Multi-patch-scale models
The statistical models developed for the concentration of each nutrient (dependent variable) incorporated species composition (five levels, including a no-macrofauna treatment), heterogeneity (two levels) and patch connectivity (three levels) as nominal explanatory variables.
(iii) Model selection
We adopted the same selection approach for both the patch-scale and multi-patch-scale models as outlined by Diggle et al.  and routinely adopted in ecology (e.g. [17,18,44]). As a first step, a linear regression model was applied to check that underlying statistical assumptions were not violated. As the validation procedure indicated normality but heterogeneity of variances, we used linear regression with a generalized least-squares (GLS) estimation procedure [45,46] that incorporates variance–covariate terms to model the variance structure. For responses at the patch scale, we accounted for the dependence among patch-scale observations within the same mesocosm by directly modelling the variance–covariance structure of the response using the within-group component (i.e. location). To determine the minimal adequate model, the optimal structure in terms of random components was determined using restricted maximum-likelihood (REML) estimation. The optimal fixed-effects structure was then determined using maximum-likelihood (ML) estimation. To determine the optimal random structure, we compared the full linear regression models to the equivalent GLS models incorporating specific variance structures using Akaike information criteria (AIC) and by inspection of model residual patterns. The optimal fixed structure was obtained by applying a backward selection using the likelihood ratio test obtained by ML estimation. Following West et al. , the numerical output of the minimum adequate model was derived using REML estimation and only the highest-order terms were interpreted . The importance of each explanatory variable was assessed by comparing a reduced model (with all terms involving the factor of interest removed) with the full model, using the likelihood ratio test. All analyses were performed using the nlme package (v. 3.1 ) in the R statistical and programming environment .
At both the patch and the multi-patch scale, there were clear effects of both species composition and habitat structure (heterogeneity, connectivity and location) on ecosystem properties. For brevity, the details of the initial and minimal adequate models are presented in the electronic supplementary material.
(a) Effects at the patch scale
(i) Net movement
Location had the greatest influence on the net movement of invertebrates between patches (L-ratio = 676.85, d.f. = 96, p < 0.0001), followed by patch connectivity (L-ratio = 173.39, d.f. = 56, p < 0.0001), species composition (L-ratio = 167.55, d.f. = 84, p < 0.0001) and habitat heterogeneity (L-ratio=118.63, d.f. = 56, p < 0.0001), but effects were strongly dependent on complex interactions between species composition and all components of habitat structure (heterogeneity × species composition × connectivity × location, L-ratio = 32.05, d.f. = 18, p < 0.05; figure 1). Overall, movement between patches was more pronounced in more connected systems (panels (a) versus (b) and (c) versus (d), figure 1) and/or where resources were more heterogeneously distributed (compare panels (a) and (b) with (c) and (d), figure 1). While the variability of species movement between patches was most strongly associated with the level of heterogeneity (HA exhibited greater variability than HB, compare panels in figure 1), the direction of movement was more consistent across all treatments. Irrespective of patch connectivity or habitat heterogeneity, and consistent with findings elsewhere (e.g. ), C. volutator tended to move out of enriched patches, while H. ulvae and H. diversicolor generally moved into the enriched patches. The biomass distribution of C. volutator between patches was more even in the more homogeneous and less connected systems. Movement of H. ulvae out of patches was most pronounced from non-enriched patches with two interfaces (location 3 or 5) where only one neighbouring patch was enriched with algae. The movement of H. diversicolor between patches was also closely related to patch enrichment, with greatest movements occurring out of non-enriched patches, especially in mesocosms with greater heterogeneity. Net movement in the species mixture reflected the behaviour of individual species, incorporating the effects of interspecific interactions, such that differences between treatments were less distinct.
(ii) Bioturbation intensity
The transport of luminophore particles approximated a biodiffusive profile with depth in all habitat patches (rmean ± s.e. = 0.019 ± 0.003, n = 149; electronic supplementary material, figure S2). Bioturbation (Db) was dependent on a two-way interaction between species composition and the level of connectivity between patches (species composition × connectivity, L-ratio = 21.69, d.f. = 8, p < 0.01; figure 2a), as well as an independent effect of location (L-ratio = 14.86, d.f. = 8, p < 0.01, figure 2b). Species composition was of greatest importance (L-ratio = 253.63, d.f. = 12, p < 0.0001), followed by connectivity (L-ratio = 30.92, d.f. = 10, p < 0.01) and location (L-ratio = 14.86, d.f. = 8, p < 0.01). When species were in monoculture and isolated in patches (0% connectivity), bioturbation intensity (coefficient ± s.e.; electronic supplementary material, table S1) was highest for C. volutator (24.50 ± 2.45, t = 9.99, p < 0.0001), followed by H. diversicolor (0.35 ± 0.08, t = 4.38, p < 0.0001) and H. ulvae (0.09 ± 0.06, t = 2.63, p < 0.01), but in all cases, increasing openness of the system decreased bioturbation. When species were in mixture, the intensity of bioturbation was high relative to monocultures of H. ulvae and H. diversicolor, and also decreased with increasing openness of the system. Irrespective of species composition, Db was lowest in the 50 and 100 per cent connectivity treatments (figure 2a). The effect of location on Db was less pronounced (figure 2b), but was significantly greater in locations 3 and 7 relative to locations 1 and 5 (p < 0.05 in all cases, except location 3 versus 5, p = 0.09).
(b) Effects at the multi-patch scale
Nutrient concentrations were dependent on the independent effects of species composition (NH4-N, L-ratio = 86.71, d.f. = 4, p < 0.0001, figure 3a; PO4-P, L-ratio = 44.63, d.f. = 4, p < 0.0001, figure 4a), patch connectivity (NH4-N, L-ratio = 80.14, d.f. = 2, p < 0.0001, figure 3b; PO4-P, L-ratio = 48.42, d.f. = 2, p < 0.0001, figure 4b) and, in the case of PO4-P, marginal effects of heterogeneity (L-ratio = 3.79, d.f. = 1, p = 0.052, figure 4c). The highest concentrations of nutrients (coefficients ± s.e.) were in treatments containing H. diversicolor (NH4-N, 1.97 ± 0.27, t = 7.17, p < 0.0001; PO4-P, 0.040 ± 0.008, t = 4.90, p < 0.0001), followed by the species mixture (NH4-N, 1.58 ± 0.19, t = 8.25, p < 0.0001; PO4-P, 0.029 ± 0.005, t = 5.49, p < 0.0001), C. volutator (NH4-N, 1.43 ± 0.20, t = 7.12, p < 0.0001; PO4-P, 0.012 ± 0.004, t = 2.80, p < 0.01) and H. ulvae (NH4-N, 0.46 ± 0.16, t = 2.91, p < 0.01; PO4-P, 0.001 ± 0.005, t = 0.30, p = 0.768). In the case of the latter, PO4-P concentration was not significantly different from treatments containing no macrofauna.
The effects of patch connectivity on nutrient concentration were consistent across both nutrients. The estimated coefficients suggest that, relative to 0 per cent patch connectivity, nutrient concentrations (coefficients ± s.e.) were highest with 50 per cent patch connectivity (NH4-N, 0.38 ± 0.16, t = 2.44, p < 0.05; PO4-P, 0.015 ± 0.003, t = 4.44, p < 0.0001) and intermediate with 100 per cent patch connectivity (NH4-N, 0.31 ± 0.16, t = 1.96, p = 0.052; PO4-P, 0.013 ± 0.003, t = 3.91, p < 0.001). For both nutrients, concentrations did not differ between 50 and 100 per cent patch connectivity treatments (NH4-N, −0.07 ± 0.16, t = −0.48, p = 0.633; PO4-P, −0.002 ± 0.003, t = −0.537, p = 0.592).
Reduction in heterogeneity (i.e. HA → HB, figure 4c) had marginal negative effects on PO4-P concentration (coefficient ± s.e. = −0.006 ± 0.003, t = −1.96, p = 0.052).
In natural ecosystems, a suite of complex associations between habitat structure and physical, chemical and biological processes operate at varying spatial and temporal scales to affect the functioning of ecosystems [3,4,30,50]. Resource availability and distribution, in particular, are fundamental in influencing community structure and species behaviour [4,51], and may also strongly modify ecosystem properties [18,19]. In contrast to previous studies that have demonstrated clear effects of biodiversity on ecosystem function under largely homogeneous conditions (e.g. [52–55]), our analyses demonstrate that species effects on ecosystem process (net movement and bioturbation intensity) and function (nutrient concentration) are non-additive, species-specific and mediated by habitat structure. Indeed, the inclusion of one or more components of habitat structure in all of our models confirms the importance of environmental complexity in mediating species interactions (e.g. [56,57]) and shaping biodiversity–function relations  (although see  for a counterexample).
At the patch scale, species identity effects on net movement and bioturbation intensity were dependent on complex interactions with habitat structure. Net movement between patches was most pronounced in more connected and more heterogeneous systems, strongly reflecting species-specific habitat preferences. Algal enrichment changes the physical and chemical characteristics of sediments and, although enrichment increases food availability for some benthic invertebrates, it can also cause enhanced sediment anoxia as a result of microbial processes (e.g. ). The observed single-species faunal responses in our study are therefore either a consequence of species-specific sensitivities to changes in the chemical characteristics of the sediment or to food availability . Corophium volutator can be very sensitive to oxygen depletion caused by algal enrichment  and therefore moves away from enriched patches [18,19]. In contrast, H. ulvae and H. diversicolor are more tolerant of anoxic sediment conditions and are capable of exploiting patchy resources. Both of these species tend to move towards enriched patches, deplete the available algal resource and then relocate to other patches that contain elevated levels of resource (sensu ). In the species mixtures used here, negative interspecific interactions may have also affected movement between patches in response to overlapping habitat use . Irrespective of the mechanism, the ability of fauna to move between patches and through corridors creates a highly dynamic system, in which species composition and density distribution (but not species richness) are continually changing, emphasizing the importance of species evenness  in maintaining ecosystem functioning . These movement dynamics underpin the observed patterns in local bioturbation intensity and probably disproportionately contribute to nutrient concentrations at the multi-patch scale because the resource-dependent relocation of infauna forms an additional source of bioturbation to that linked with other activities, such as feeding and burrow maintenance. Faunal movement, coupled with the effects of density dependence in preferential patches (e.g. [65,66]), resulted in the relocation of some fauna (especially C. volutator) to a less desirable habitat (the corridors), thereby reducing patch bioturbation intensity. Indeed, corridors may form a viable habitat if the corridor quality is perceived to be higher than that of the other surrounding habitat . The patchy distribution of resources among the mesocosms therefore resulted in individual species favouring different patches, leading to spatial turnover in intra- (or inter-) specific density distribution that negatively affected the patch-scale ecosystem process. In a recent meta-community study, France & Duffy  demonstrated that allowing mobile grazers to move and select preferential patches reduced the spatial variability in grazer abundance, which, in turn, altered the level of ecosystem productivity. Hence, the movement of species between patches may alter local and regional processes by changing the composition and density distribution of species assemblages [27,30,68]. Here, overall nutrient concentration was positively affected by enhanced patch connectivity, presumably because in more open systems, movement between patches probably resulted in more extensive burrow systems.
Habitat fragmentation may not necessarily have strong effects on community structure if the scale of the patch–corridor system relates to the size and mobility of the organisms . Here, the dimensions of the alternative habitat patches, and the distance between those patches, matched those observed in natural systems (S. M. Lawrie 1996, unpublished PhD thesis), but we did not detect reduced levels of infaunal movement with decreasing corridor width in addition to the effects of differential habitat quality. Meta-community experiments using below-ground soil decomposer organisms also demonstrate weak effects of habitat fragmentation (reviewed in ), but both marine sediment and terrestrial soil systems are known to host a highly heterogeneous distribution of decomposer resources (e.g. [35,71]), such that decomposer species may be less susceptible to short-term changes in habitat structure. Indeed, Matthiessen et al.  proposed that on small spatial and temporal scales, diversity loss (or the loss of a dominant consumer species) will have a stronger effect on local ecosystem processes when species dispersal plays a minor role. Here, a patch-scale ecosystem process (bioturbation intensity) was greater in the more closed experimental system, while ecosystem functioning at the multi-patch scale was positively affected by increasing openness of the system.
Collectively, our results indicate that local community dynamics associated with habitat structure may be of greater importance in moderating ecosystem processes and/or function than have previously been appreciated, and these emergent effects are likely to be dependent on scale and ecosystem process/ecosystem function. While these observations emphasize the importance of the biodiversity–habitat structure couple for maintaining multi-functional ecosystems , they also highlight the context dependency of this coupling [3,8], and reinforce the need for experiments aimed at understanding the long-term functional consequences of altered biodiversity and community dynamics in less connected habitats . If we are to fully appreciate the ecological consequences of biodiversity loss under future global change, studies need to embrace these dynamics and account for uncertainties generated by the modifying effects of multiple interacting biotic and abiotic stressors.
Supported by the Natural Environment Research Council (awarded to J.A.G., grant NER/S/A/2005/13734). We thank Leigh Murray, Anne Holford, Martha Higuera, Martha Pascual and Carlos Sanz Lazaro for assistance in the field and help with setting up the experiments. We are also grateful to A. Jamieson for the technical drawing.
- Received November 5, 2010.
- Accepted December 13, 2010.
- This Journal is © 2011 The Royal Society