Size matters: implications of the loss of large individuals for ecosystem function

Size is a fundamental organismal trait and an important driver of ecosystem functions. Although large individuals may dominate some functions and provide important habitat structuring effects, intra-specific body size effects are rarely investigated in the context of BEF relationships. We used an in situ density manipulation experiment to explore the contribution of large, deep-burrowing bivalves to oxygen and nutrient fluxes across the sediment-water interface. By manipulating bivalve size structure through the removal of large individuals, we held species identity constant, but altered the trait characteristics of the community. The number of large bivalves was the best predictor of ecosystem functioning. Our results highlight that (a) accounting for body size provides important insights into the mechanisms underpinning biodiversity effects on ecosystem function, and (b) if local disturbances are recurrent, preventing individuals from reaching large sizes, the contribution of large adults may be lost, with largely unknown implications for ecosystem functionality.

disturbance are rarely random 31 and large organisms are often vulnerable. The life-stages of individual species differ in their potential recovery following disturbance 32 . Bivalves have decreasinging mobility with increasing size and it is common for the small early life-stages to dominate recovery, while adult stages take considerable time to establish through growth. While complete extinctions of regional species pools are comparatively rare, compositional changes and reductions in abundance and biomass in the degradation process are common, so that recovering populations, while contributing to species richness contribute little to ecosystem function 33 . Increased mortality of one species to below its ecologically effective population size (EEP), while not making this species go extinct, may indeed have functional effects resulting in the extinction of other species instead 34 .
In marine systems very little attention has been directed towards changes in ecosystem function in the community assembly process following disturbance (but see 35 ). Importantly, while speciesabundance patterns may exhibit comparatively fast recovery 36 , communities may take considerable time to develop populations with undisturbed demographic characteristics 32 . Hence, if local disturbances are recurrent the contribution of large adult stages may be lost, with largely unknown implications for ecosystem functionality.
In the Baltic Sea, structural and functional biodiversity is naturally reduced due to low salinity, and the critical role of the few functional groups is apparent as losses of any species may entail a loss of the only representative of a function, such as suspension feeding 37,38 . It thus provides an ideal environment for empirical testing of key traits for ecosystem function. In addition very few benthic species in the Baltic Sea are long-lived or large, i.e. with traits that are predicted to have important influences on ecosystem function. The shallow soft-sediment communities are typically comprised of only a handful of species; the Baltic clam Macoma balthica and the soft-shell clam Mya arenaria typically make up an average of 15% of total community abundance and 75%, or more, of community biomass. Observations from recent manipulative field experiments in subtidal soft-sediment habitats suggest that community assembly processes following disturbance may result in substantial transient dominance shifts, with relatively quick recovery in terms of both species numbers and abundances 36 . Nevertheless, our observations also suggest that the recovery of mature and large-sized components of the bivalve populations may take considerable time (several years).
We conducted a field experiment to test the overall hypothesis that large adult bivalves are foundation species in soft-sediment communities with profound influences on ecosystem function. After disturbance it takes a long time for these adult bivalves to re-establish. Our prediction was that the contribution of bivalves to ecosystem function would mirror their relative dominance in terms of biomass. We tested this prediction by conducting an in situ density manipulation experiment where we (1) disturbed a community to eliminate all fauna to initiate a community assembly process where species-abundance patterns would have recovered (i.e. after 12 mo), but where large, mature life-stages would be lacking, and (2) seeded large individuals of bivalves (Macoma and Mya) to undisturbed control communities to obtain elevated densities (still within their natural range). We then incubated the sediment in situ to examine the contribution of deep-burrowing adult stages of large bivalves to measures of ecosystem functioning: ammonium and phosphate fluxes at the sediment-water interface, and community respiration. These measures are key ecosystem functions in soft-sediment habitats. We demonstrate that body-size is a key organism trait with important implications for understanding BEF relationships.

Results
The two different sampling occasions were combined in our analysis of both macrofauna and fluxes to increase replication of our study. Bottom water temperature was 19uC at T1 and 13uC at T2. Water column concentrations of oxygen and nutrients differed slightly : 0.24 6 0.05 vs. 0.30 6 0.10 mmol l 21 , respectively). We also observed variability in macrofauna community structure between sampling occasions. This variability, however, merely added strength of inference to our findings.
Treatment effects on benthic community structure. After the experimental disturbance (Fig. 1), the manipulated plots were left to recover for a year. After 12 months we observed no differences in diversity and the average numbers of taxa were more or less identical across treatments (Table 1).There were, however, differences in total community abundance and biomass, and in the distribution of bivalves (Table 1). Interestingly, the highest community abundance values were observed in the disturbed plots (D), which was also the case for bivalves in general and Macoma in particular. This may seem counterintuitive as bivalves were added to the elevated plots (E); however these differences could be explained by higher densities of post-settlement juvenile bivalves and polychaetes in plots with a disturbance history. Nevertheless, differences in total biomass and bivalve biomass were very clear, with the lowest biomasses observed in the recovering community (D) and the highest in the treatment to which bivalves had been added (E), with biomass values ranging from 56 to 460 g wwt m 22 (Table 1).
Multivariate analyses showed that overall community structure was not different between treatments for abundance, as no clear groupings were detected (Fig. 2a). In contrast, patterns of community structure in terms of biomass were distinctly different between treatments, forming clear groupings and a gradient from D, C to E (Fig. 2b). Interestingly the disturbed plots (D) exhibited the largest variability between replicates. The multivariate PERMANOVA analysis confirmed the significance of these patterns and showed that abundance variations indeed were non-significant (Pseudo-F 5 1.53, p 5 0.126), while groupings in biomass (Fig. 2b) were highly significant overall (Pseudo-F 5 14.68, p 5 0.001) and also between all treatments ( Table 2).
The SIMPER analyses identified bivalves as contributing most to differences between sample clusters observed for biomass. The average group dissimilarity between C and D was 73% and the bivalves contributed . 60% of this difference. The average group dissimilarity between C and E was 72% and here bivalves (Macoma and Mya) contributed 94% of the difference between treatments. Overall community dissimilarity between D and E was, as expected, the largest at 93% and also here bivalves contributed most to the difference, 90%. Other common taxa were hydrobid gastropods and the errant polychaete Hediste diversicolor, which were important for the within-group variability in, especially, the disturbed treatment (D). The average within-group dissimilarity was largest in the disturbed treatment (D, 77%), and smallest in the Elevated treatment (E, 20%).
Relationship between macrofauna and nutrient fluxes -and the contribution of bivalves. Dark chambers were used to measure the net flux of oxygen and nutrients across the sediment-water interface in the absence of primary production.
Bivalves were the dominant drivers of community biomass patterns and explained 98% of biomass variability (r 2 5 0.98; p , 0.0001, linear regression) and were hence expected to drive ecosystem function relationships. Indeed, bivalve biomass and the number of large bivalves explained more of the variability in O 2 consumption, NH 4 1 and PO 4 32 -fluxes than total residual community biomass (Fig. 3, Table 3). PO 4 32 -fluxes were exceedingly low, as is typical for sandy sediments low in organic matter.
In DistLM marginal tests, solute fluxes were correlated most strongly with the number of adult bivalves ( Table 3). The only other significant predictors were the number of juvenile bivalves and number of species, which were weakly correlated with O 2 flux. However in partial tests (i.e. after correcting for the influence of large bivalves) neither of these two variables was significant. The number of large bivalves was the best linear predictor of PO 4 32 and NH 4 1 flux, explaining 53 and 79% of the variability, respectively. Including other variables in the model only explained an additional 3-5% of the variation. For the O 2 flux, a combination of adult and juvenile bivalves and the number of species explained 19% more of the variation (cumulative r 2 5 0.56) than a model containing just the number of adult bivalves. Interestingly, total community abundance had no significant effect on these fluxes.

Discussion
We have demonstrated that intraspecific variations in body-size can be a key predictor of ecosystem functioning. We used a density manipulation experiment to explore the contribution of large deepburrowing bivalves to oxygen and nutrient fluxes across the sediment-water interface, important measures of ecosystem function. We defaunated the seafloor a year in advance to initiate the community assembly process and observed recovery in terms of speciesabundance distribution, but as expected observed only very limited recovery of macrofaunal biomass (Fig. 2, Table 2). Important members of the benthic community, such as adult polychaetes and gastropods colonized the disturbed plots, and juvenile bivalves were also observed in high numbers, but adult bivalves remained more or less absent (Table 1). To increase density variations in our experiment, we also added adult bivalves to undisturbed plots, and showed that bivalves dominate measures of ecosystem function. While the Elevated treatment (E) made an important contribution to the observed response, the bivalve densities in these plots were still within the natural density variation observed in the area. Importantly, we also observed distinct differences between disturbed and control plots (D vs. C), and the distribution of biomass across treatments formed a clear gradient (Fig. 2b). In contrast, species numbers and abundances were not significantly different between treatments.
An important goal in ecology, and for successful restoration and conservation, is to understand how species contribute to ecosystem processes, such as the rate and stability of nutrient cycling 1 . In softsediment habitats, shifts in ecosystem performance are often associated with changes in species influencing organic matter recycling and nutrient regeneration 22,39,40 . We used in situ flux chambers to examine the role of large individuals and their contribution to community respiration and sediment nutrient fluxes. In the analysis we  www.nature.com/scientificreports ignored the categorical treatments and simply explored the relationship between biomass and large individuals and measures of ecosystem function across treatments (Fig. 3). Again, neither species diversity nor abundance could explain much of the variability in the observed responses. Our results show that the number of large bivalves were the strongest predictors of ecosystem function (Fig. 3).
These results support earlier studies, reporting that presence of bivalves enhances benthic respiration and the release of ammonium through bioturbation and excretion 23,26 . Generally bioturbation can enhance the amount of fresh organic material transported into the   1 fluxes are, however, in addition to sediment reworking also due to bivalve excretion 23,26,44 . Still, species-specific traits are likely to affect sediment redox-dependent processes in different ways and result in complex biogeochemical interactions. For example, both Macoma and Mya are sediment biodiffusers. The more shallow-burrowing Macoma is positioned in the sediment nitrification zone, and may enhance NO 3 2 efflux to the overlying water. In contrast, the deeper-dwelling Mya arenaria, transfers oxygen into the reduced zone of the sediment and may enhance nitrification-denitrification rates and thus cause an uptake of NO 3 2 45 . Hence animal-sediment interactions are complex and might result in different impacts on nutrient regeneration processes, depending on the biology and trait-composition expressed by the resident species. Nevertheless, our results unequivocally demonstrate that large bivalve individuals are strong predictors of ecosystem function. Although it would be intresting to conduct additional experiments where the same high biomass is made up of a large number of small indiviudals, such a situation is not likely to exist in natural bivalve beds. In addition, large bivalves often bury deeper in the sediment and can be expected to displace more sediment, pump more water and create stronger pore-water pressure gradients 17 .
Our study is one of the first to partition the contribution of large individuals to important measures of ecosystem function through an in situ manipulation of a real community. In fact, by manipulating bivalve size structure through the removal and addition of large individuals, we held species identity more or less constant, but altered the trait characteristics and functional diversity of the community 9 . Other taxa (e.g. polychaetes) at our study site are more fast growing than bivalves and were able to attain normal biomasses over the oneyear recovery period. Our results highlight that without the presence of large adults, ecosystem functionality is radically changed. Bivalve species such as Macoma balthica and Mya arenaria have life-spans of 6-10 and 10-20 years, respectively, and maximum life-spans of 30 years have been reported for both species 46,47 . This indicates that while the regional supply of bivalve larvae and post-settlement juveniles may result in rapid colonization into disturbed habitat patches 36,48 , mature stages will take years to recover, especially since adult infaunal bivalves have limited mobility and recovery is thus largely dependent on individual growth. Inter-and intra-specific traits such as longevity and large size are disproportionately affected by habitat loss and too frequent disturbance regimes 49 . Increases in disturbance-regimes are hence of particular concern as the community assembly processes may be interrupted before bivalves reach full size and are able to contribute to important ecosystem processes. Indeed, historical reconstructions have highlighted that losses of suspension-feeding bivalves have profoundly influenced food webs and ecosystem function 50 .
In soft-sediment systems the degradation of macrobenthic communities as a result of disturbance has been shown to result in the loss of deep-burrowing large taxa and is predicted to reduce bioturbation 31 . Eutrophication-induced hypoxia and anoxia has spread widely across the world 51 and is particularly common in the Baltic Sea in both coastal and open-sea waters 52 . The consequent loss of deep-burrowing and bioturbating taxa, and particularly their adult life-stages is of concern because of their major influence on all oxygen-dependent biogeochemical processes 41,43 . As highlighted by Ellison et al. (2005) 16 , the dynamics of communities shaped by foundation species, such as the bivalves in our system, may be dominated by a small number of strong interactions, which makes these types of communities fragile to switching between alternative stable states. In such communities disturbances have the potential to flip the ecosystem across a threshold into a different stability domain, and the probability for this to happen increases as foundation species are driven to regional functional extinction 53 .
Changes in biodiversity affect ecosystem functioning and can thus disrupt the way ecosystems contribute to valuable ecosystem services (e.g. nutrient regeneration processes 2 ). Biodiversity losses typically involve declines in both abundance and biomass of common species, thus shifting dominance patterns of communities 1,54 . BEF studies have, however, mostly focused on species richness effects even though reported species identity or ''sampling'' effects are common, indicating that particular dominant traits may be underpinning ecosystem function 6,9,54 . We show that individual body-size is important for ecosystem functionality. As highlighted by Bengtsson (1998) 55 , body-size distributions have mechanistic links to ecosystem functions (e.g. energy flow and nutrient cycling), because most rates of ecosystem processes are mechanistically related to biomass through uptake, feeding and physiology. In contrast, the mechanistic link between species diversity and process rates is less clear. Our in situ findings support hypotheses put forward 9,10 and experimental findings 20 that evaluation of body size, not only between but also within species, provides important insights into the mechanisms behind biodiversity effects on ecosystem function. Importantly, we show that large individuals in natural communities may have a major influence on ecosystem function. If local disturbances are recurrent, preventing individuals from reaching large sizes, the contribution of large adult stages may be lost, with severe implications for ecosystem functionality 34 . The characterization of body size and its importance for communities and ecosystem function has ramifications for conservation and restoration efforts, because it facilitates the interpretation of how disturbances, through the functional elimination of species (i.e. the loss of large-sized individuals), might propagate through the ecosystem 56 .

Methods
Study area. There are few areas where the loss or degradation of both habitat and species diversity is as evident as in the Baltic Sea. Hypoxic zones cover up to 70.000 km 2 that are largely devoid of all benthic macrofauna 52,57 and it is clear that the reduction in the distribution and diversity of Baltic Sea benthos due to hypoxic events has already altered the way benthic ecosystems contribute to key ecosystem processes (i.e. nutrient cycling). Recently, the problem of seasonal hypoxia in shallower, nearshore areas has also been highlighted 37,51,52 . Our experiment was conducted in the northern Baltic Sea, near Tvärminne Zoological Station (59u 509 440 N, 23u 149 960 E), Finland. The site was at 4 m depth and had sandy sediments, with a median sediment grain size of 0.29 mm, organic matter content of 0.5 6 0.03% (SD), and total carbonand nitrogen content of 0.18 6 0.02 and 0.02 6 0.01%, respectively. Salinity is around 6 and there are no tides in the area.
Experimental setup. The experiment included three treatments; the disturbed (D), the control (C) and one treatment with an elevated number of adult bivalves (E). One year before the flux measurements, three blocks were established along a 50 m transect line. For the disturbed treatment, a 16 m 2 plot was defaunated in each block by covering the sediment surface with black LDPE plastic to induce anoxia to underlying sediments. The disturbance manipulation simulated patchy hypoxia, for example induced naturally by drifting algal mats 37 . Plots were covered for a period of 16 days, to ensure complete defaunation and the plastic was removed in late July 2006 (details in 36 ). Following the experimental manipulation, large dead bivalves were observed to have emerged to the sediment surface, a common escape response to hypoxia (Fig. 1). One year later, the control and elevated treatments plots were placed . Measurements of sediment nutrient and oxygen fluxes. Measurements of sediment nutrient and oxygen fluxes were performed with one dark benthic chamber (504 cm 2 , volume 6 l) for each treatment replicate (i.e. 18 chambers in total; see 43 for methodological details). Chamber frames were deployed (pushed 6.5 cm into the sediment) one day before the incubation started. Simultaneously, large bivalves were added to the elevated treatment and allowed to rebury over night, with nets keeping potential predators from foraging on the bivalves. Incubation started once chamber lids were installed and ended 6 hours later. Water samples were taken from the chambers at start and end of the incubation. The chamber water was manually stirred with a paddle from the outside before water samples were withdrawn with syringes (in total 200 ml) using a sampling port in the chamber lid. The initial water volume of the tube was discarded, before the chamber water was collected. Replacement water was supplied through a port that was placed distant from the sampling port. To correct for water column effects, three dark 1 l LPDE bottles were used for incubation of ambient water at 4 m depth during the experiment. Water samples were processed immediately on the boat. For determination of dissolved oxygen, 60 ml of each sample was fixed with 0.5 ml Mn(OH) 2 and 0.5 ml KI. The rest of each sample was filtered through a Whatman GF/F filter, directly into a 250 ml Nalgene bottle for nutrient analyses. All samples were stored on ice in the dark during transport to the laboratory and nutrient samples were then frozen (220uC) until further analysis. Dissolved oxygen concentrations were determined by the Winkler procedure, while NH 4 1 and PO 4 32 were measured with Lachat flow injection analysis.
Sampling of benthic fauna and sediment. Samples for sediment organic matter, total carbon (TC) and total nitrogen (TN) were taken with a 2.1 cm diameter core from the control and disturbed sediments. In the laboratory, the surface sediment (upper 1 cm) was stored at -20uC until analysis. Sediment organic matter was determined by loss on ignition (3 h at 500uC). Sediment for nutrient analyses was freeze-dried and thoroughly homogenized. Analyses of TC and TN in sediments were performed with a Carlo Erba high temperature combustion elemental analyzer.
After the incubation, chamber lids were removed and one macrofaunal core (Ø 5.6 cm, depth 15 cm) sample was taken from the middle of each chamber. The area enclosed by each chamber was then excavated (to 30 cm depth), in order to account for any deeper-burrowing bivalves. Large excavated bivalves (. 5 mm) were counted and measured (shell length and wet weight, including shell).
The macrofaunal core samples were preserved in 70% ethanol and stained with Rose Bengal. To account for small recruits, the samples were elutriated by first suspending the sediment in a bucket of spinning water and decanting off the supernatant through a 200 mm sieve (repeated five times), and then checking the remaining sediment for any larger animals. The fauna was identified, counted and measured at 10 3 magnification. To obtain the proportions of juveniles and adults of dominant taxa, shell lengths of bivalves and gastropods, and the width of the 10th setiger of Hediste diversicolor and Marenzelleria sp. were measured. Gastropods , 1 mm shell length were only identified to family. The total weight of each species was determined (precision 0.1 mg blotted wet weight, including shell).
Data analysis. Multivariate analyses of benthic community data were performed with the PRIMER software 59 . Multidimensional scaling (MDS) and PERMANOVA were used to identify differences in abundance and biomass between treatments, while the SIMPER analysis was used to identify species contributions to (dis)similarities within and between treatments. Bray-Curtis similarity coefficient was based on untransformed data. Distance-based linear models (DistLMs) in PERMANOVA1 59 were used to determine if macrofauna variables were significant predictors of O 2 consumption, PO 4 32 and NH 4 1 fluxes. Predictor variables included the abundance of juvenile (, 5 mm shell length) and adult bivalves, residual community abundance and biomass (i.e. less the contribution of bivalves) and the number of species. Marginal tests were run to identify strong, significant predictors irrespective of other variables, then partial tests were performed to assess the explanatory value of a predictor variable after other variables had been accounted for. The number of large bivalves and bivalve biomass were strongly collinear and therefore only the former was included. Remaining predictor variables were weakly correlated with each other (Pearson's r , 0.4). Similarity matrices of untransformed fluxes were constructed using Euclidean distances and the model was run using the step-wise selection procedure and R 2 selection criterion. P values were obtained for predictor variables by 99999 permutations.