Philosophical Transactions of the Royal Society B: Biological Sciences
You have accessArticles

Marine biodiversity–ecosystem functions under uncertain environmental futures

Mark T. Bulling

Mark T. Bulling

Oceanlab, University of Aberdeen, Main Street, Newburgh, Aberdeenshire AB41 6AA, UK

[email protected]

Google Scholar

Find this author on PubMed

,
Natalie Hicks

Natalie Hicks

Gatty Marine Laboratory, University of St Andrews, East Sands, St Andrews, Fife KY16 8LB, UK

Google Scholar

Find this author on PubMed

,
Leigh Murray

Leigh Murray

Oceanlab, University of Aberdeen, Main Street, Newburgh, Aberdeenshire AB41 6AA, UK

Google Scholar

Find this author on PubMed

,
David M. Paterson

David M. Paterson

Gatty Marine Laboratory, University of St Andrews, East Sands, St Andrews, Fife KY16 8LB, UK

Google Scholar

Find this author on PubMed

,
Dave Raffaelli

Dave Raffaelli

Environment Department, University of York, Heslington, York YO10 5DD, UK

Google Scholar

Find this author on PubMed

,
Piran C. L. White

Piran C. L. White

Environment Department, University of York, Heslington, York YO10 5DD, UK

Google Scholar

Find this author on PubMed

and
Martin Solan

Martin Solan

Oceanlab, University of Aberdeen, Main Street, Newburgh, Aberdeenshire AB41 6AA, UK

Google Scholar

Find this author on PubMed

    Abstract

    Anthropogenic activity is currently leading to dramatic transformations of ecosystems and losses of biodiversity. The recognition that these ecosystems provide services that are essential for human well-being has led to a major interest in the forms of the biodiversity–ecosystem functioning relationship. However, there is a lack of studies examining the impact of climate change on these relationships and it remains unclear how multiple climatic drivers may affect levels of ecosystem functioning. Here, we examine the roles of two important climate change variables, temperature and concentration of atmospheric carbon dioxide, on the relationship between invertebrate species richness and nutrient release in a model benthic estuarine system. We found a positive relationship between invertebrate species richness and the levels of release of NH4-N into the water column, but no effect of species richness on the release of PO4-P. Higher temperatures and greater concentrations of atmospheric carbon dioxide had a negative impact on nutrient release. Importantly, we found significant interactions between the climate variables, indicating that reliably predicting the effects of future climate change will not be straightforward as multiple drivers are unlikely to have purely additive effects, resulting in increased levels of uncertainty.

    1. Introduction

    Anthropogenic activities continue to cause significant global climate change (IPCC 2007). Partly as a consequence of climate change, and also as a direct result of human activities (e.g. harvesting, habitat destruction and pollution), there has been a significant and rapid decline in global biodiversity (Vitousek et al. 1997; Sala et al. 2000). The speed of this decline has been estimated to be over 1000 times greater than rates of extinction calculated from the fossil record (Lawton & May 1995; Pimm et al. 1995, Millennium Ecosystem Assessment 2005). The recognition that ecosystems provide humans with goods and services that are vital for our survival and well-being has led to substantial concern over how such ecosystem services may be affected by both climate change and the decline in biodiversity (Naeem et al. 1994; Tilman et al. 1997; Kinzig et al. 2002; Hooper et al. 2005, Millennium Ecosystem Assessment 2005; Balvanera et al. 2006; IPCC 2007).

    Two of the potentially most important effects of climate change on ocean ecosystems are increasing temperatures (Hughes 2000; IPCC 2007) and rising atmospheric CO2 concentrations (Crowley & Berner 2001; Caldeira & Wickett 2003; Freely et al. 2004; IPCC 2007). The oceans act as a significant temperature sink (Levitus et al. 2005), absorbing an estimated 80 per cent of the additional heat generated by global warming (IPCC 2007). This has led to an increase of 0.1°C in the average temperature of the upper 700 m of the oceans over the last 40 years (IPCC 2007). The rate of sea temperature rise is predicted to increase over the next decades, with the estimated rise in the mean surface temperature for the end of the twenty-first century lying between 1.1°C and 6.4°C (IPCC 2007).

    The current concentration of atmospheric CO2 is approximately 380 ppm, having increased from preindustrial levels of 280 ppm (IPCC 2007). By the year 2100, assuming a ‘business-as-usual’ scenario, atmospheric CO2 concentration is predicted to be between 450 and 1000 ppm (IPCC 2007). The world's oceans act as a CO2 sink, with dissolved CO2 forming carbonic acid which dissociates into protons and bicarbonate ions, causing a decrease in ocean pH, a process known as ‘ocean acidification’ (Caldeira & Wickett 2003; Freely et al. 2004). Surface ocean pH is currently 0.1 unit lower than the preindustrial value, and by the end of the century is likely to be a further 0.14–0.35 units lower, leading to changes greater than any experienced in the past 300 Myr (Caldeira & Wickett 2003; Freely et al. 2004; IPCC 2007).

    Impacts of climate change on ecosystems, however, are unlikely to occur as isolated single drivers. With multiple impacts, the question then arises as to whether they are likely to affect ecosystems in an additive way or whether there will be significant interactions between the drivers (Folt et al. 1999; Christensen et al. 2006; Fabry et al. 2008; Przeslawski et al. 2008; Byrne et al. 2009). Understanding these relationships has important implications for our ability to predict the interactions between climate change, biodiversity loss and ecosystem services. To date, research on the effects of increasing temperature and ocean acidification on marine ecosystems has focused on single drivers in isolation, and there has been a particular emphasis on specific effects such as coral calcification (for overviews see Fabry et al. 2008; Przeslawski et al. 2008). Thus, there is currently a disconnect between the literature examining the biodiversity–ecosystem functioning (BEF) relationship and the literature examining the effects of climate change on ocean ecosystems (Schmitz et al. 2003), as well as a general lack of studies examining interactive effects of climate change on ecosystem functioning. The latter is particularly important because the interactions between multiple stressors are difficult to predict as the resulting net effect may not be additive but, instead, may either be greater (synergistic effects) or less (antagonistic effects) than anticipated (Folt et al. 1999).

    As a step towards understanding the uncertainty associated with current forecasts of the ecological consequences of environmental change, here we test for interactive effects of macrofaunal species richness, temperature and atmospheric CO2 concentration on levels of nutrient release from the sediment to the water column in a model marine benthic system. An important mechanism underpinning nutrient release is the process of bioturbation, the mixing of porewater solutes and sediment particles by the movement and activities of the benthic organisms (Richter 1952). Bioturbation affects sediment permeability, breaks down chemical gradients in pore waters and subducts organic matter, thus influencing rates of remineralization and inorganic nutrient efflux (Gray 1974; Rhoads 1974; Kristensen & Blackburn 1987). Benthic habitats process and recycle a significant proportion of nutrients that are necessary for primary production in coastal waters and hence underpin the delivery of ecosystem services from these systems, with ammonium being of particular importance in nitrogen-limited marine waters (Logan et al. 1995). Quantifying the impact of temperature and carbon dioxide, and their interactions, on the relationship between benthic macrofaunal biodiversity and ecosystem functioning is therefore of considerable importance for predicting the effect of climate change on these systems (Thrush & Dayton 2002; Lohrer et al. 2004).

    2. Material and methods

    We used a mesocosm approach to investigate the effects of changes in temperature and CO2 on ecosystem functioning in a model three-species estuarine benthic system. Our experimental design incorporated all possible combinations of three invertebrate macrofauanal species exposed to different levels of atmospheric CO2 and temperature. Three levels of atmospheric CO2 were used, ranging from present day levels (380 ppm) to predicted future levels (600 and 1000 ppm) that represent low-middle and upper values in the range of the estimates of the Intergovernmental Panel on Climate Change (IPCC) for the year 2100 (IPCC 2007). These levels of atmospheric CO2 were arranged in an orthogonal design with three temperature levels (6°C, 12°C and 18°C) that are within the range experienced by these species throughout the annual cycle. All combinations of species and environmental variables were replicated three times. Ecosystem functioning was assessed by measuring levels of NH4-N and PO4-P in the water column at the end of the experiment (7 days). In order to assess the mechanistic effects of different climate scenarios we also measured levels of bioturbation activity using fluorescent tracer particles (luminophores; Mahaut & Graf 1987), following the methodology described in Solan et al. (2004).

    (a) Mesocosms

    Sediment was collected from the Ythan estuary, Aberdeenshire, Scotland, UK (57°20.085′ N, 02°0.206′ W), sieved (500 μm) in sea water to remove unwanted macrofauna, and left to settle for 24 h, allowing us to also retain the fine fraction (less than 63 μm). Excess water was removed, and the sediment slurry was homogenized and distributed between mesocosms (Perspex cores 33 cm high with an internal diameter of 10 cm) to a depth of 10 cm (equivalent to 785 cm3), thus minimizing any effects of sediment heterogeneity between mesocosms and ensuring the absence of macrofauna.

    Sea water (UV-sterilized, 10 μm pre-filtered, salinity ≈ 33) was initially added to the mesocosms to a depth of 20 cm, left for 24 h and then refilled with sea water to eliminate nutrient pulses associated with assembly (Ieno et al. 2006) and to ensure that changes in nutrients could be attributable to activity during the experimental period. Mesocosms were aerated throughout the experimental period (7 days).

    (b) Regulation of atmospheric CO2, temperature and light

    Mesocosms were placed in environmental chambers (VC 4100, Vötsch Industrietechnik), which maintained a constant temperature environment (±0.1°C). Levels of atmospheric CO2 within the chamber were maintained (±30 ppm) using a CO2 monitor in the chamber connected to a valve system on a standard CO2 cylinder (BOC Gases Ltd., UK) via a digital controller (Technics horticultural carbon dioxide controller). The regulation of atmospheric CO2 was calibrated and validated using an Infrared gas analyser (ADC LCA3). Light was supplied by two 400 W metal halide bulbs (Newlec) in each chamber, on a 12 h light—12 h dark cycle. Mesocosms were arranged randomly within a chamber to minimize the effect of any spatial heterogeneity in light levels.

    (c) Macrofauna

    Three species of macrofauna were used: the gallery-forming Hediste diversicolor (Polychaeta), the surficial modifier Hydrobia ulvae (Gastropoda) and the regenerator Corophium volutator (Crustacea). These species represent a broad range of taxonomic groups, functional effects and functional responses and they fulfil key structural and functional roles within benthic communities. They dominate the study site both in terms of abundance and biomass, but have contrasting mobility, modes of bioturbation and bioirrigation. They also have differing effects on ecosystem functioning, particularly with respect to setting redox depth and controlling carbon and nutrient cycling. Macrofauna were collected from the Ythan estuary and stored in aerated sea water tanks for 24 h before they were placed in the mesocosms (day 0).

    (d) Measures of ecosystem functioning

    Pre-filtered (Nalgene, 0.45 μm) water samples were taken on the final day of the experiment. NH4-N and PO4-P concentrations were determined using standard protocols with a modular flow injection auto-analyser (FIA Star 5010 series) using an artificial sea water carrier solution.

    (e) Measures of ecosystem process

    We distributed 2 g of 90–125 μm luminophore particles evenly across the sediment surface 24 h after the introduction of the macrofauna. At the end of the experiment, a sediment core (3.5 cm diameter and 10 cm deep) was taken from each mesocosm and sliced (0.5 cm thick for the uppermost 2 cm, and 1 cm thick from 2 to 10 cm depth). The number of luminophore particles within each slice was determined using standard image analysis techniques based on thresholding (Solan et al. 2004). These counts provided a vertical luminophore profile for each mesocosm which could be used to determine a bioturbation index of activity (Db; Diaz & Cutter 2001; Gilbert et al. 2003; Solan et al. 2004). Here, Db is an index of the rate of redistribution of the luminophores due to bioturbatory activity and is expressed as an estimate of the increase in the variance of depths to which the luminophores would have moved over the time scale of a year.

    (f) Species richness treatments

    Within each environmental combination there were eight species richness treatments comprising all species richness combinations: a control (no macrofauna), three monocultures (H. diversicolor, Hd; H. ulvae, Hu; C. volutator, Cv), three two-species treatments (Hd + Hu, Hd + Cv, Hu + Cv) and one three-species treatment (Hd + Hu + Cv). This design minimizes hidden treatment effects (sensu Huston 1997) and eliminates pseudo-replication. The repetition (n = 3) of each permutation allows the generality of any diversity effects to be evaluated. A complementary analysis could also be performed where species combination replaces species richness as an independent factor, allowing ‘species combination’ effects on ecosystem function or process to be detected. To ensure that any observed changes in ecosystem process were directly attributable to treatment manipulations and not to species density, total biomass was held constant at 2.0 g per mesocosm (equivalent to 255 g−2), a level consistent with that found at the study site (e.g. Biles et al. 2003). As nutrient cycling is primarily a microbial process that is mediated by macrofaunal bioturbation, replicate cores in the absence of macrofauna were required so that the contribution of other components of the sediment (microbial and meiofaunal communities) could be distinguished from those of the macrofauna. There were three replicates for each species—environment (temperature × CO2) treatment combination, giving a total of 216 mesocoms.

    (g) Analysis

    A generalized least squares (GLS; Pinheiro & Bates 2000; Zuur et al. 2007) statistical mixed modelling approach was used, treating NH4-N, PO4-P and Db as dependent variables, and levels of CO2, temperature and species richness (or species combination) as independent fixed factors. A GLS framework was preferred over linear regression using transformed data because it retains the structure of the data while accounting for unequal variance in the variance–covariate terms. In each case, as a first step, a linear regression model was fitted. Model validation showed no evidence of nonlinearity but there was evidence of unequal variance among the explanatory variables. The GLS framework was then adopted in order to model this heterogeneity of variance. The most appropriate random structure was found by examination of AIC scores in conjunction with plots of fitted values versus residuals for models with different variance–covariate terms relating to the independent variables, using restricted maximum likelihood (REML, West et al. 2007). The fixed component of 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 REML. Following Underwood (1998), the highest potential level of interaction that was assessed was the three-way interaction, and nested levels within higher order interactions were not examined. 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 including it were omitted. 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

    (a) NH4-N

    The minimal adequate model with species richness as an independent term was a linear regression model with a GLS extension (see the electronic supplementary material), and incorporated all three two-way interaction terms: species richness × temperature, CO2 × temperature and species richness × CO2 (table 1). Species richness was the most influential variable (L-ratio = 150.26, d.f. = 15, p < 0.0001), followed by temperature (L-ratio = 65.08, d.f. = 12, p < 0.0001), and CO2 (L-ratio = 45.71, d.f. = 12, p < 0.0001).

    Table 1.Summary of significant terms found in the linear regression models with a generalized least squares extension, treating levels of NH4-N, PO4-P and particle bioturbation (Db) as dependent variables, and species richness (or species identity), atmospheric CO2 concentration and temperature as fixed explanatory variables.

    dependent variable significant terms L-ratio d.f. p
    species richness models
    NH4-N species richness × temperature 23.96 6 <0.001
    CO2 × temperature 19.40 4 <0.001
    species richness × CO2 12.39 6 0.054
    PO4-P temperature 116.14 2 <0.0001
    CO2 23.25 2 <0.0001
    Db species richness × temperature 30.25 6 <0.0001
    species identity models
    NH4-N species identity × temperature 67.14 14 <0.0001
    species identity × CO2 32.70 14 <0.01
    CO2 × temperature 26.10 4 <0.0001
    PO4-P species identity × temperature 26.72 14 0.021
    CO2 30.64 2 <0.0001
    Db CO2 × species identity ×temperature 60.11 28 <0.001

    There was a general increase in NH4-N concentration with increasing species richness, and this was most rapid at the intermediate temperature of 12°C (figure 1). At present day concentrations of CO2, the highest levels of NH4-N were obtained at 18°C. However at the higher CO2 concentrations, this peak in NH4-N concentrations generally occurred at 12°C.

    Figure 1.

    Figure 1. Predicted NH4-N concentrations from the minimal adequate regression model for varying species richness levels at three atmospheric CO2 concentrations (380 ppm (solid line), 600 ppm (dashed line) and 1000 ppm (dotted line)) and the three temperatures (a) 6°C, (b) 12°C and (c) 18°C.

    Within the pattern of increasing NH4-N concentration with increasing species richness, there was a decrease in concentration with increasing level of CO2. This effect was greatest at the higher species richness levels. At 18°C there was a marked decrease between NH4-N levels under present day CO2 concentrations and those under the two possible future CO2 concentrations. This coincided with a general flattening of the relationship between species richness and NH4-N concentration at 18°C.

    The minimal adequate model treating the particular combinations of species as an independent variable (species identity) was also a linear regression model with a GLS extension (see the electronic supplementary material) incorporating all three two-way interaction terms: species identity × temperature, species identity × CO2 and CO2 × temperature (table 1). Species identity was the most influential variable (L-ratio = 204.16, d.f. = 35, p < 0.0001), followed by temperature (L-ratio = 108.68, d.f. = 20, p < 0.0001), and CO2 (L-ratio = 73.14, d.f. = 20, p < 0.0001).

    The overall trend of increasing NH4-N concentration with increasing species richness could be discerned in this model (figure 2). However, within species richness levels, there was a marked variation in NH4-N concentrations, particularly at the two lower temperatures. Higher concentrations of NH4-N were found when C. volutator, H. ulvae or both were present. However, this pattern was absent at 18°C where marked species identity effects were absent.

    Figure 2.

    Figure 2. Predicted NH4-N concentrations from the minimal adequate regression model for each combination of species (Hd, Hediste diversicolor; Hu, Hydrobia ulvae; Cv, Corophium volutator) at three atmospheric CO2 concentrations (380 ppm (solid line), 600 ppm (dashed line) and 1000 ppm (dotted line)) and the three temperatures (a) 6°C, (b) 12°C and (c) 18°C.

    (b) PO4-P

    Levels of PO4-P in relation to species richness were best explained by a linear model with a GLS extension (see the electronic supplementary material) with only single additive terms for temperature and CO2 (table 1).

    Estimated coefficients for this model (figure 3) indicated that there was no significant difference between the concentrations of PO4-P at 6°C and 12°C, but concentrations at these temperatures were significantly different to those at 18°C. Estimated coefficients for the effect of CO2 (figure 3) indicated significantly different concentrations of PO4-P at each level of CO2 concentration, with the greatest PO4-P concentration at 600 ppm and the lowest at 1000 ppm. However, the actual differences in PO4-P concentrations across these CO2 concentrations were very small (less than 0.025 mg l−1).

    Figure 3.

    Figure 3. Predicted PO4-P concentrations from the minimal adequate regression model at three atmospheric CO2 concentrations (380, 600 and 1000 ppm) and the three temperatures 6°C (solid line), 12°C (dashed line) and 18°C (dotted line). This minimal adequate model resulted from the full model which included species richness in the full model. Effects of temperature and CO2 were additive with no significant interaction term.

    The resulting model for the dependent variable PO4-P with species identity as an independent variable (see the electronic supplementary material) contained the two-way interaction term species identity × temperature, and the single term CO2 (table 1). The most influential independent variable was temperature (L-ratio = 161.15, d.f. = 13, p < 0.0001), followed by species identity (L-ratio = 52.57, d.f. = 21, p < 0.001) and CO2 (L-ratio = 30.64, d.f. = 2, p < 0.0001).

    There were marked species identity effects on the levels of PO4-P within the two-way interaction term species identity × temperature (figure 4a). Patterns in these effects were generally similar between temperatures of 12°C or 18°C, but these were different to patterns in concentration observed at 6°C. At the higher temperatures, C. volutator and combinations of C. volutator and H. ulvae tended to be associated with higher concentrations of PO4-P, while combinations of C. volutator and H. diversicolor were associated with low concentrations of PO4-P. However, at the lowest temperature, treatments containing C. volutator were associated with lower concentrations of PO4-P, while the highest PO4-P concentrations were linked to the H. ulvae monoculture treatment. Concentrations of PO4-P were similar at 6 and 12°C, although the data showed a high degree of variation, and these PO4-P concentrations were generally significantly above those observed at 18°C. Coefficient estimates for the effect of CO2 (figure 4b) reveal that the highest concentrations of PO4-P were observed where CO2 was 600 ppm, followed by 380 ppm and then 1000 ppm.

    Figure 4.

    Figure 4. (a) Predicted PO4-P concentrations from the minimal adequate regression model (with species identity rather than species richness as an independent variable) for each combination of species (Hd, Hediste diversicolor; Hu, Hydrobia ulvae; Cv, Corophium volutator) at the three temperatures 6°C (solid line), 12°C (dashed line) and 18°C (dotted line). Atmospheric CO2 concentration is held constant. (b) Estimates from the minimal adequate regression model for PO4-P concentrations of coefficients associated with the effects of the higher atmospheric CO2 concentrations relative to the effect at 380 ppm. Error bars represent standard errors.

    (c) Particle bioturbation (Db)

    The minimal adequate model for Db as a dependent variable and with species richness as an independent variable simply consisted of the two-way interaction term species richness × temperature (table 1). The most influential variable was species richness (L-ratio = 185.73, d.f. = 9, p < 0.0001) followed by temperature (L-ratio = 37.56, d.f. = 8, p < 0.0001).

    There was a general decline in Db with increasing species richness (figure 5). However, the Db values for monocultures were markedly higher at 6°C and 12°C than at 18°C, making the rates of decline in Db with species richness much greater at these temperatures.

    Figure 5.

    Figure 5. Predicted Db (cm2 yr−1) levels from the minimal adequate regression model for varying species richness levels at the three temperatures 6°C (solid line), 12°C (dashed line) and 18°C (dotted line).

    The minimal adequate model with species identity as an independent variable (see the electronic supplementary material), contained the three-way interaction term CO2 × species identity × temperature (table 1). Species identity was the most influential variable (L-ratio = 327.11, d.f. = 63, p < 0.0001), followed by temperature (L-ratio = 148.14, d.f. = 48, p < 0.0001) and CO2 (L-ratio = 77.94, d.f. = 48, p < 0.01).

    In the monoculture treatments with a CO2 concentration of 380 ppm and a temperature of 6°C, the greatest Db values were associated with H. diversicolor (figure 6), with values for the other two species being markedly lower. However, when temperature or CO2 levels increased these Db values decreased dramatically, and at the highest CO2 concentrations all Db values for H. diversicolor monocultures were on par with those for the other species. Even under environmental conditions where Db was high for H. diversicolor monoculture treatments, treatments involving H. diversicolor in combination with other species tended to result in much lower Db values.

    Figure 6.

    Figure 6. Predicted Db levels from the minimal adequate regression model for each combination of species (Hd, Hediste diversicolor; Hu, Hydrobia ulvae; Cv, Corophium volutator) at three atmospheric CO2 concentrations (380 ppm (solid line), 600 ppm (dashed line) and 1000 ppm (dotted line)) and the three temperatures (a) 6°C, (b) 12°C and (c) 18°C.

    4. Discussion

    When data from BEF experiments are integrated (Balvanera et al. 2006; Cardinale et al. 2006; Schmid et al. 2009), the relationships between biodiversity and ecosystem functioning are mostly positive. We found a positive relationship between species richness and NH4-N concentration incorporating strong species identity effects, particularly H. diversicolor, but no relationship between species richness and PO4-P concentration, although species identity effects were found. These results generally agree with previous studies on this particular benthic estuarine system (Emmerson et al. 2001; Dyson et al. 2007; Bulling et al. 2008) and the results for NH4-N are broadly in line with findings elsewhere (e.g. Hansen & Kristensen 1997; Mermillod-Blondin et al. 2005). However, Caliman et al. (2007) did find a positive relationship between PO4-P concentrations and species richness in freshwater mesocosm experiments.

    Significant species composition effects were found for both measures of nutrient release (NH4-N and PO4-P) and for bioturbation (Db), and significant species richness effects were found for NH4-N and Db, although not for PO4-P. A cursory consideration of the differences between the two models for PO4-P may be interpreted as contradictory conclusions on the relative importance of biodiversity. However, the ability to incorporate species composition in the modelling process (achieved here because of the complete experimental design) has revealed mechanistic complexities that would otherwise remain hidden, including which species disproportionately affects functioning.

    The general negative relationship between species richness and bioturbation (Db) was modified by temperature. However, the model incorporating species composition as an independent variable revealed that this relationship was primarily driven by H. diversicolor. Monocultures of H. diversicolor generally resulted in markedly higher levels of bioturbation but this effect decreased in the presence of other species and with increasing temperature and, generally, with increasing CO2 concentration. This dominance of H. diversicolor is consistent with previous results (Ieno et al. 2006), and is a consequence of H. diversicolor's ability to deeply burrow into the sediment. As Db is strongly influenced by the depth of the luminophore profile obtained, the decline in Db at higher levels of species richness is likely to be due to a change in behaviour of H. diversicolor, reducing the extent of its movement into deeper sediment and altering its behaviour (bioirrigation more than particle displacement) when other species were present and in relation to temperature levels and CO2 concentrations (Ouelette et al. 2004). Experimental work examining the effect of sea water acidification on the polychaete worm Nereis virens (Widdicombe & Needham 2007) found the reverse effect of increasing NH4-N concentrations with decreasing pH. Burrow characteristics were unchanging between pH treatment levels and it was concluded that the impact of sea water pH on nutrient flux was probably due to changes in the microbial community (Bertics & Ziebis 2009). The shorter duration of our experiment (7 days) relative to that for N. virens (five weeks), however, means that changes in microbial community composition are unlikely to fully explain the patterns observed here. A more parsimonious explanation is a combination of context-dependent switching behaviour (Riisgård & Kamermans 2001) and biogeochemical effects related to lowered pH levels (Wood et al. 2009). Context-dependent changes in behaviour have been shown in experiments investigating the effect of temperature on the polychaetes Nereis virens (Deschênes et al. 2005) and Neanthes virens (Ouellette et al. 2004), with a general increase in activity levels occurring with increasing temperature. Ouellette et al. (2004) found the highest Db values at intermediate temperatures (13°C), but they also found increasing biotransport levels with increasing temperature. Similarly, Deschênes et al. (2005) found changes in behaviour with temperature, with locomotion being strongly inversely correlated with temperature.

    The change in behaviour of H. diversicolor and the resultant consequences for ecosystem process and function in our experiment has important implications for the understanding of BEF relationships under future climate scenarios. There was a change in functional response, indicating that functional groups (Petchey & Gaston 2002; Naeem & Wright 2003) may be fluid in nature, with some species changing their behaviour in response to shifting environmental conditions. Functional groupings defined under current conditions may not be a reflection of functional groups in the future, suggesting that current levels of functional redundancy may not be a reliable representation of redundancy in the future, a problem that is compounded when multiple functions are considered. This fluidity has the potential to significantly alter BEF relationships, and suggests that the precautionary principle (Myers 1993) may have an increased relevance when managing such systems for the delivery of ecosystem services for the future under changing climatic conditions.

    In similar marine systems, Emmerson et al. (2001) found variability in the contributions of particular species to ecosystem function resulting in a consistent, but idiosyncratic, increase in function with increasing richness. Results here suggest that this concept can be extended to incorporate extra dimensions describing climate change variables. Overall patterns between species richness and function were consistent despite strong species composition effects, but differences in temperature and/or CO2 concentration could elicit idiosyncratic changes in functioning levels due to the variability in the behavioural responses of the species to these changes in climatic variables.

    To date, the majority of research on the effects of global climate change has been focused on single drivers (Fabry et al. 2008; Przeslawski et al. 2008; Byrne et al. 2009), despite recognition that multiple drivers may act in non-additive ways (Christensen et al. 2006; Byrne et al. 2009; Laurance & Useche 2009). Reich et al. (2001) found interactive effects between atmospheric CO2 and biodiversity, and between nitrogen deposition and biodiversity on biomass levels in a terrestrial plant system. However, they did not find significant interactions involving both variables CO2 and nitrogen. Our findings strongly suggest that there will be significant interaction effects between climate change variables and biodiversity on levels of future ecosystem functioning and services. We found significant interactions between temperature, atmospheric CO2 concentration and species richness and/or species composition determining concentrations of NH4-N. There was a general trend of increasing NH4-N concentration (ecosystem functioning) with increasing species richness but this tendency was modified by the two climate change variables. Our results also suggest that the effects of climate change variables on the BEF relationship may be nonlinear for certain ecosystem functions. We found a distinct peak in the steepness of the species richness—NH4-N concentration relationship at the intermediate temperature. It is notable that there was variability in the complexity of the models developed. Those for NH4-N concentrations involved three two-way interaction terms. Of the two models for Db, only the model incorporating species identity contained a three-way interaction term, and PO4-P concentrations were driven by the climatic variables in an additive manner or via a single two-way interaction.

    Given the already high uncertainty associated with predictions of future climatic conditions and the expected regional variations (IPCC 2007), the presence of interactive effects between climatic drivers and biodiversity, as well as contrasting and nonlinear effects on ecosystem functioning, must be major concerns in terms of prediction and management of the consequences of climate change. In assessing the probable ecological consequences of anthropogenic and environmental change, we need to embrace the concept that multiple drivers, and their interactions and feedbacks, are likely to modify the natural environment in ways that are difficult to predict. There is therefore a requirement for a more holistic framework that aims to understand how multiple drivers, alone and in concert, affect multiple functions (Elmqvist et al. 2003; Hector & Bagchi 2007; Gamfeldt et al. 2008).

    The orders of overall influence of independent variables within a model were consistent across all models for NH4-N and Db, with species richness or species composition being markedly the most influential, followed by temperature, and finally CO2 concentration. However, for the species richness model for PO4-P, temperature was most influential followed by CO2 concentration, with species richness not having a significant influence. The corresponding model with species combination as an independent variable, however, incorporated both climatic drivers and biodiversity (temperature being most influential, followed by species composition and CO2 concentration), highlighting the importance of the roles of individual species. Differences in the order and strength of influence of the climatic and species variables for the different ecosystem functions and processes considered here also emphasize the need for a more holistic approach for developing future management strategies in the light of specific mitigation priorities.

    It is important to note that under two levels of CO2 concentration the fauna was exposed to conditions that are likely to be found 100 years into the future. This longer time scale would allow other mechanisms to play an influential role, in particular changes in species' distributions and evolutionary adaptation processes. The study was also run for a short time (7 days) which limits the available time for species to acclimate. In this context it is important to recognize that this experiment was not designed to examine and predict specific physiological and behavioural responses of organisms to climate change and the resulting changes in ecosystem functioning. Rather, our model system is a simplification of a natural community and is therefore amenable to manipulation (Benton et al. 2007), allowing us to explicitly test the sensitivity of the BEF relationship to climate change variables and potential interaction effects.

    Our results, based on benthic marine systems, have shown that the relationships between ecosystem functioning, biodiversity and climate change are not straightforward, and that they are mediated by interactions between the different drivers within the system. Future BEF research needs to focus on the mechanistic interactions of multiple drivers and multiple functions, and on the development of more holistic approaches to underpin future management strategies.

    Acknowledgements

    The authors would like to thank O. McPherson for technical assistance, and T. Fujii for help in the field. They would also like to thank two anonymous reviewers for helpful and constructive comments. This work was supported by NERC grant NE/E006 795/1.

    Footnotes