Research ArticleECOLOGY

Habitat diversity and ecosystem multifunctionality—The importance of direct and indirect effects

See allHide authors and affiliations

Science Advances  08 Feb 2017:
Vol. 3, no. 2, e1601475
DOI: 10.1126/sciadv.1601475


Ecosystems worldwide are facing habitat homogenization due to human activities. Although it is commonly proposed that such habitat homogenization can have negative repercussions for ecosystem functioning, this question has yet to receive explicit scientific attention. We expand on the framework for evaluating the functional consequences of biodiversity loss by scaling up from the level of species to the level of the entire habitats. Just as species diversity generally fosters ecosystem functioning through positive interspecies interactions, we hypothesize that different habitats within ecosystems can facilitate each other through structural complementarity and through exchange of material and energy across habitats. We show that experimental ecosystems comprised of a diversity of habitats show higher levels of multiple ecosystem functions than ecosystems with low habitat diversity. Our results demonstrate that the effect of habitat diversity on multifunctionality varies with season; it has direct effects on ecosystem functioning in summer and indirect effects, via changes in species diversity, in autumn, but no effect in spring. We propose that joint consideration of habitat diversity and species diversity will prove valuable for both environmental management and basic research.

  • Habitat diversity
  • multifunctionality
  • direct and indirect effects
  • Marine
  • Sediment
  • microbial diversity
  • SEM


Ecosystems around the globe are facing habitat homogenization due to human activities (13). Homogenized ecosystems will reduce the diversity of species and consequently diminish valuable ecosystem functions and services (4). In addition to affecting species diversity, reduced habitat diversity may also directly impair ecosystem functioning because of reduced cross-habitat exchanges of material and energy. Understanding the functional consequences of biodiversity loss across multiple scales of organization, that is, at both habitats and species levels, is thus critical. However, studies on the functional consequences of changes in biodiversity have been solely designed to investigate the effects of diversity among and within species (taxonomic or functional) (4). The effects of changes in habitat diversity on ecosystem functioning remain unexplored.

Here, we scale up the framework for consideration of biodiversity and ecosystem functioning by investigating whether ecosystems comprised of a diversity of habitats have higher levels of functionality than ecosystems with low habitat diversity (Fig. 1). Our experiment directly addresses the issue of potential effects of habitat homogenization on ecosystem functioning. It merges theory on biodiversity and ecosystem functioning (4) with landscape ecology (5) and thus incorporates the concept of meta-ecosystems (6). We propose that, like interspecies interactions, habitats can facilitate each other, and that ecosystem-wide functioning is promoted by habitat complementarity. Complementarity among habitats includes the exchange of material and energy, such as various forms of oxidants and reductants (for example, oxygen and organic matter). The effect of habitat diversity can be direct but also indirect via changes in species diversity (Fig. 1). One example of interactions among habitats is the interplay between mangroves and coral reefs in tropical ecosystems (7). Mangroves serve as nursery habitats that affect the community structure and biomass of fish on neighboring coral reefs (8), which, in turn, may protect mangrove by functioning as a breakwater (9). Similar positive interactions may be common in ecosystems that are composed of a mosaic of different habitats.

Fig. 1 Conceptual diagram of our framework.

Ecosystem homogenization (caused by, for example, human disturbance) results in a change in habitat diversity (A). Because habitats have different physical and chemical characteristics, they are likely associated with different sets of species. Loss of habitat diversity thus potentially leads to loss in species diversity (the union of the species in all habitats, indicated by different symbols) (B). Changes in habitat diversity can affect ecosystem functioning not only directly through changes in structural complexity and the cross-habitat exchange of nutrients and other resources (C) but also indirectly via changes in species diversity.

Natural ecosystems perform many ecosystem functions simultaneously; all of these functions have the potential to be positively or negatively affected by biodiversity (10). To quantify the overall effect of habitat diversity on ecosystem functioning across a range of functions, we therefore used an index of multifunctionality (10), which summarizes the focal ecosystem functions that are examined in our study. Because habitat diversity generally favors species diversity (11), the net effect of habitat diversity on ecosystem multifunctionality was partitioned into direct and indirect effects using structural equation modeling (SEM). Moreover, the relationships between habitat diversity, species diversity, and ecosystem multifunctionality may differ across seasons. Therefore, we also tested how changes in habitat diversity affect ecosystem multifunctionality during different seasons.

As model systems, we used experimental ecosystems consisting of four different natural habitats common in coastal marine ecosystems: cyanobacterial mats, plant meadows (the seagrass Ruppia maritima), silty mud, and sandy beach. Cores with sediment and overlying water from each habitat were sampled in the field and arranged randomly into ecosystems to form a diversity gradient, including one, two, three, or four habitat types. This setup allowed for interactions between and within the habitats via a common water column. Given the importance of shallow-water coastal systems in terms of productivity and nutrient cycling (12), we measured four key biogeochemical processes that describe the productivity and nitrogen cycling in these ecosystems: gross primary production (GPP), nitrogen fixation, denitrification, and uptake of dissolved inorganic nitrogen (DIN). We hypothesized that higher habitat diversity within an ecosystem (i) directly enhances ecosystem multifunctionality and (ii) increases bacterial and microalgal diversity, and that (iii) an increase in bacterial and microalgal diversity also increases multifunctionality. To examine whether these hypotheses are valid across seasons, we repeated the experiment in spring, summer, and autumn, representing different combinations of light, temperature, and concentrations of organic matter and inorganic nutrients.


Habitat and bacterial and algal diversity

On the basis of habitat descriptors (see Materials and Methods), a permutational multivariate analysis of variance (PERMANOVA) showed that the degree of interhabitat difference varied with season (season × habitat interaction; P = 0.001) such that habitat types differed from each other during summer and autumn (Psummer = 0.001, Pautumn = 0.002), but not in spring (Pspring = 0.214) (Fig. 2A). The bacterial community structure of all four habitats was different between seasons (P = 0.006), but not between habitats within each season (Pspring = 0.506, Psummer = 0.961, Pautumn = 0.261; Fig. 2B). The habitat diversity of a given ecosystem was calculated with Euclidean distances for each habitat (based on chemical, physical, and biological characteristics) and was thus defined as the sum of all pairwise distances between the constituting habitats. Habitat diversity was then scaled between 0 and 1 (turning the habitat diversity into a continuous variable) by dividing all values (across habitat diversity levels and seasons) by the maximum total Euclidean distance (see Materials and Methods). The diversity of the bacterial communities (fig. S2) increased with increasing habitat diversity during summer (P = 0.027) and autumn (P = 0.005), but decreased in spring (P = 0.031) (Fig. 3A and table S1). Benthic microalgal diversity (fig. S3) did not correlate with habitat diversity during any season (summer, P = 0.196; autumn, P = 0.944; spring, P = 0.340; table S2). We also provide additional analyses with habitat diversity as a categorical variable with four levels (1, 2, 3, and 4), which gave similar results as those described for habitat diversity as a continuous variable (see the Supplementary Materials, “Results with habitat diversity as a categorical variable” section).

Fig. 2 Principal components analysis and nonmetric multidimensional scaling plots on habitat descriptors and OTUs.

The intra- and interhabitat variability and variability of bacterial community structure during spring, summer, and autumn are shown. (A) Ordination of habitat samples based on habitat descriptors displayed in a principal components analysis (PCA) plot with Euclidean distances; n = 48. (B) Bacterial community structure (OTU-based) displayed in a nonmetric multidimensional scaling (NMDS) plot based on weighted UniFrac distances; n = 48. Color codes indicate the habitat types Sandy beach (light brown), Silty mud (dark brown), Cyanobacterial mats (blue), and Ruppia maritima meadows (green).

Fig. 3 Linear functions of relationships between habitat diversity, microbial diversity, and multifunctionality across seasons.

(A) Relationship between habitat diversity and bacterial diversity, (B) habitat diversity and index of multifunctionality (weighted average value of standardized functions), (C) bacterial diversity and index of multifunctionality, and (D) benthic microalgal diversity (effective number of species) and index of multifunctionality. Shaded areas indicate ±95% confidence interval; ntot = 84, nlevel 1 = 16, nlevel 2,3,4 = 4 per season.

Temporal importance of habitat, bacterial, and microalgal diversity for multifunctionality

Both habitat diversity and bacterial diversity significantly affected ecosystem multifunctionality, but the effects varied with season (Fig. 3). During spring, neither habitat diversity (P = 0.208; Fig. 3B; fig. S4, A and B; and table S3) nor bacterial diversity (P = 0.680; Fig. 3C and table S4) nor algal diversity (P = 0.926; Fig. 3D and table S5) affected ecosystem multifunctionality, which was also supported by the results from the SEM (Fig. 4A and table S6). In summer, ecosystem multifunctionality increased with habitat diversity (P = 0.002; Fig. 3B; fig. S4, A and B; and table S3) and bacterial diversity (P = 0.028; Fig. 3C and table S4), but not with algal diversity (P = 0.773; Fig. 3D and table S5). In summer, the SEM identified a direct effect of habitat diversity on multifunctionality (Phabitat diversity→multifunctionality = 0.0001), that is, not meditated through increased bacterial diversity (Fig. 4B and table S7). In autumn, we observed no direct effects of habitat diversity (P = 0.304; Fig. 3B; fig. S4, A and B; and table S3), bacterial diversity (P = 0.942; Fig. 3C and table S4), or algal diversity on multifunctionality (P = 0.821; Fig. 3D and table S5). However, there was an indirect effect of habitat diversity through increased bacterial diversity (Phabitat diversity→bacterial diversity = 0.010, Pbacterial diversity→multifunctionality = 0.039; Fig. 4C and table S8). Among the four individual functions supporting multifunctionality, GPP (P = 0.008; table S9 and Fig. 5A), nitrogen fixation (P = 1.6 × 10−9; table S10 and Fig. 5B), and uptake of DIN (P = 0.041; table S11 and Fig. 5C) significantly increased with habitat diversity during summer, but not in autumn or spring. Denitrification showed no trend of habitat diversity effects (table S12 and Fig. 5D). Analyses with habitat diversity as a categorical variable show similar results (see the Supplementary Materials, “Results with habitat diversity as a categorical variable” section), and the effects of habitat diversity on multifunctionality were independent of the method used to calculate multifunctionality (fig. S4, A and B).

Fig. 4 Structural equation models.

Path diagrams based on SEM showing how habitat and bacterial diversity affect ecosystem multifunctionality during (A) spring, (B) summer, and (C) autumn. Solid paths (blue) are statistically significant (P < 0.05) with standardized path coefficients indicated, whereas the dashed gray lines are not. Percentages indicate the variance explained by the model; ntot = 84, nlevel 1 = 16, nlevel 2,3,4 = 4 per season.

Fig. 5 Linear models of individual functions.

Linear models of individual functions used to calculate multifunctionality against habitat diversity during spring, summer, and autumn are shown. (A) GPP (mmol O2 m−2 day−1), (B) nitrogen fixation (mmol N2 m−2 day−1), (C) uptake of DIN (ammonium and nitrate + nitrite) (mmol DIN m−2 day−1), and (D) denitrification (nmol N g wet sediment−1 hour−1). Shaded areas indicate ±95% confidence interval; ntot = 84, nlevel 1 = 16, nlevel 2,3,4 = 4 per season.

Net effects of habitat diversity on individual functions and multifunctionality

A net effect [sensu Loreau and Hector (13)] of habitat diversity was observed in summer for GPP (P = 0.0005; Fig. 6D), nitrogen fixation (P = 0.01; Fig. 6B), and multifunctionality (P = 0.008; Fig. 6E). Adjusting the P values for multiple comparisons within each season left GPP (P = 0.002) and multifunctionality (P = 0.034) as significant factors but weakened support for nitrogen fixation (P = 0.057) (Fig. 6, A to E). In spring, a negative net effect on GPP was observed (P = 0.027, Padjusted = 0.137), but no other significant net diversity effect was observed in spring or autumn. A positive net effect means that the mixture of four habitat types has higher functioning than what could be expected from the functioning of single habitats. For example, in summer, the mean observed multifunctionality was 0.58, and the expected multifunctionality based on the individual habitats was 0.35. Thus, the net diversity effect was 0.23, or 66% higher than expected.

Fig. 6 The net habitat diversity effect on individual ecosystem functions and multifunctionality.

“Expected” is the expected functionality in the treatment based on each of the single habitats, and “observed” is the observed functionality. (A) Denitrification (nmol N g wet sediment−1 hour−1). (B) Nitrogen fixation (mmol N2 m−2 day−1). (C) Uptake of DIN (ammonium and nitrate + nitrite) (mmol DIN m−2 day−1). (D) GPP (mmol O2 m−2 day−1). (E) Index of multifunctionality (weighted average value of standardized functions). The points are slightly spread along the x axis (grouped by season) and jittered (within season) for clarity. The triangles represent group means.


In landscape ecology, spatial heterogeneity is a central and causal factor of ecological systems. The spatial arrangement of different elements in a landscape matrix affects fluxes of energy and material (5), but the field of biodiversity–ecosystem functioning has not yet incorporated this organizational aspect into experimental investigations. The importance of scale and spatial heterogeneity has been discussed in verbal (14, 15) and theoretical (6) frameworks, but a concept for experimental tests has not been described. Our experiment demonstrates that habitat diversity directly and indirectly drives ecosystem multifunctionality and that habitat homogenization can threaten ecosystem multifunctionality beyond the loss of species diversity. We propose that our general framework (Fig. 1), in which we scale up from the level of species to the level of the entire habitats, provides a profitable way forward if net ecosystem consequences of environmental change are to be modeled and managed.

In our experiment, multifunctionality during summer was larger than the expected sum of all individual habitats, with the highest habitat diversity showing an observed level of multifunctionality that was 66% higher than expected from single habitat treatments (Fig. 6E). This is the equivalent to overyielding at the level of species; that is, polycultures are more productive than monocultures due to positive interactions, such as complementarity (13). Habitats can be complementary in terms of both their species composition and their physical structure (58). Biogeochemical and structural characteristics of shallow-water ecosystems can differ significantly from each other, a circumstance that allows for potentially positive interactions (16). For example, if nitrogen fixation is favored in one habitat, organismal growth may be supported in adjacent habitats with less available nitrogen. Moreover, Ruppia meadows affect oxygen availability in the surrounding sediment and overlying water (17), which connected the habitats in our experiment. Oxygen originally produced by Ruppia, for example, may therefore stimulate mineralization of organic nitrogen and further oxidation to nitrate in habitats that are typically less oxygen-rich, such as silty mud. However, in contrast to species diversity, diversity effects of habitats can only be attributed to positive interactions between habitats and not to niche complementarity or selection effects [sensu Loreau and Hector (13)]. While species in habitat mixtures can increase in relative abundance and thereby drive diversity effects, the proportion of each habitat in this experiment was constant over time. Each type of habitat constituted 25% of the total surface in the four-habitat mixture over the course of the experiment. Therefore, the observed diversity effect cannot be explained by high-performing single habitats. Furthermore, there cannot be any niche complementarity among habitats, as there can be for species, because there are no “habitats” for habitats. Although a habitat can harbor different numbers of species, total habitat cover cannot deviate from 100%. Therefore, we attribute ecosystem-level overyielding in our study to positive interactions among habitats. As shown in fig. S5, GPP is considerably higher in mixtures containing three and four habitat types compared to habitats containing a single habitat type. GPP must therefore be influenced by strong positive interactions. The flux of inorganic nitrogen and nitrogen fixation shows similar, although weaker, patterns. In contrast, denitrification in mixtures is the average of what it is in habitat “monocultures,” signaling no positive interactions. The fact that multifunctionality and single functions in the habitat mixtures deviate from what would be expected on the basis of the functioning in the single habitat treatments (Fig. 6 and fig. S5) provides evidence that the four habitats did exchange material and energy via the common water column.

The structure of the four habitat types in our experiment was physically, chemically, and biologically different during the three seasons. Temperature and light (fig. S6) and concentrations of inorganic nutrients (fig. S7) and organic material (fig. S8, A and B) all displayed seasonal differences, which affected the habitat-defining properties. In spring, the lack of effects on multifunctionality was probably a consequence of the high within- and low between-habitat dissimilarity (Fig. 2A) in combination with low water temperature and bacterial diversity due to the hyperdominance of a single operational taxonomic unit (OTU) most closely affiliated with Pseudomonas sp. The structural dissimilarity of the different habitats in spring indicates physical sediment disruption, causing habitat homogenization. This is a direct consequence of the winter season conditions, such as erosion by freezing, ice, storms, low light, and low availability of organic nutrients. All these environmental factors induce variation in the biogeochemical properties in sediments of the individual habitats (18). By contrast, the four habitats were clearly separated in summer due to increased growth of autotrophic components, such as Ruppia meadows, cyanobacterial mats and well-developed diatom mats, higher temperature, and more stable weather conditions. Thus, structural differences could underlie the observed direct relationship between habitat diversity and multifunctionality.

In summer, habitat diversity directly affected multifunctionality, while in autumn multifunctionality was mediated via bacterial diversity. Hence, our study also adds to a growing body of evidence showing the importance of indirect effects in mediating ecosystem processes (1921). Indirect effects have previously been found to be at least as important as direct effects in structuring communities (20, 21). Comparing the results of the linear models (estimating net effects) and SEM (partitioning net effects into direct and indirect effects) provided insights into the importance of the indirect effects in our experiment. For example, both bacterial community diversity and habitat diversity were significantly correlated to ecosystem multifunctionality in summer, suggesting that both aspects of diversity were driving the relationship. However, when bacterial diversity was analyzed simultaneously with habitat diversity in a SEM framework, only habitat diversity directly drove multifunctionality. This is in accordance with recent literature, which shows that biogeochemistry can be the single most important driver for the functioning of bacteria-dominated systems (22) and that changes in bacterial diversity generally have weak effects on ecosystem functioning; only around 25% of the dilution-to-extinction experiments reviewed by Roger et al. (23) found positive relationships [although there are also recent studies showing strong effects of bacterial diversity; for example, Delgado-Baquerizo et al. (24)]. However, in autumn, the effect of habitat diversity on ecosystem multifunctionality was indirect and mediated through changes in bacterial diversity. Thus, both structural complementarity and bacterial diversity were linked to ecosystem multifunctionality, but seasonal changes in ecosystem components, such as biogeochemistry and physical structure (18, 25), may have affected their relative importance. This finding illustrates the value of using statistical methods (for example, SEM) that partition between direct and indirect effects if the consequences of changes in biodiversity within ecosystems are to be explained and predicted.

Benthic microalgal diversity was unrelated to ecosystem multifunctionality during all three seasons. Because it was unfeasible to identify all benthic diatoms visually to species level, identification methods might have been insufficiently precise to detect potential differences in species diversity. A general positive effect of species richness of primary producers on production has previously been shown (4, 26), but the consequences of changes in benthic microalgal diversity are not well understood. An observational study that related benthic microalgal diversity to functioning found both positive and negative relationships (27). The only study so far in which the species richness of benthic diatoms was manipulated reported positive effects of diversity on production (28). However, the maximum richness was only eight species, and natural communities of benthic diatoms are far more diverse (27). Diversity and community composition of other types of microbiota, such as meiofauna and protozoa, might play a role in ecosystem multifunctionality. However, a full mechanistic elucidation of all relationships between the biotic and abiotic components of habitat diversity on ecosystem-scale multifunctionality was beyond the scope of our study and warrants further research.

The present study demonstrates a direct link between habitat diversity and ecosystem-scale multifunctionality, an association that was partly mediated by effects associated with bacterial diversity. From a management perspective, our results support the concept (Fig. 1) that habitat homogenization can have negative consequences for ecosystem functioning and the ecosystem services these functions underpin. Environmental management often focuses on the ecosystem and habitat level (29). However, the continuing focus in biodiversity-functioning research on species, at the expense of other levels of diversity, is unlikely to provide managers with the full range of information that they need (30). In managing complex ecosystems efficiently, simplification of ecosystem knowledge is often necessary. Our paper thus provides a framework for scaling up biodiversity–ecosystem functioning research to simplify ecological information. Studying the diversity and composition of habitats instead of, or in concert with, species diversity has the potential to provide more relevant data for management decisions.

In conclusion, our findings advance the core finding of biodiversity–ecosystem functioning research that species diversity often plays a central role in influencing the magnitude and stability of ecosystem functioning. We demonstrate that the diversity of habitats within an ecosystem can complement species diversity and even independently influence multifunctionality. An aspect of habitat diversity that is not considered in our study is the dispersal and migration of organisms between the habitats. Examples of these phenomena include the observed large-scale synergistic effects between coral reefs and mangrove ecosystems (8) and the small-scale dispersal that affects the biodiversity-functioning relationship within metacommunities (31). Nonetheless, to the extent that our study can be extrapolated to other systems, our results suggest that functions may be negatively affected if habitat diversity declines, even though species diversity is sustained.


Experimental design

To experimentally investigate the direct and indirect effects of habitat diversity on ecosystem multifunctionality, natural, intact sediment habitats of “Sandy beach,” “Silty mud,” “Cyanobacterial mats,” and “Ruppia maritima meadows” were sampled on the west coast of Sweden during the summer and autumn of 2013 and spring 2014 using plastic cylinders [inner diameter (ID) = 8 cm, height (h) = 11 cm]. In total, 112 habitat cores were sampled per season (Sandy = 29, Silty = 29, Cyano = 26, Ruppia = 28), and randomly assembled into four diversity levels (1 to 4), by placing four habitat cores in one larger cylinder (ID = 25 cm, h = 25 cm), representing one replicate. For diversity level 1, one ecosystem consisted of four cores from the same habitat type, and each ecosystem was replicated four times (4 × 4 = 16 ecosystems in total). The ecosystems with habitat diversity levels 2 to 4 consisted of a mixture of different habitat types: 2 + 2 for diversity 2, 1 + 1 + 2 for diversity 3, and 1 + 1 + 1 + 1 for diversity 4. The cylinders containing the experimental ecosystems were placed in a greenhouse with a continuous flow of surface water pumped from an adjacent bay (fig. S9). Each experiment was run for 2 weeks, which is enough time for the sediment habitats to interact and recover from eventual disturbance caused by sampling, but short enough to minimize changes in the microbial community under the experimental conditions used (32).

Habitat descriptors

Ten environmental factors were determined in sediments sampled in situ in four replicates per habitat and season. Porosity was calculated as the percentage weight loss following drying of ~30 g of wet sediment to constant weight (24 hours). The density of wet sediment was measured on 20 ml of sediment. The organic carbon and nitrogen content was analyzed in ~30 mg of dried sediment through vapor-phase acidification using a Carlo Erba NA 1500 elemental analyzer. To extract sediment pore water, 50 ml of sediment was centrifuged (32,000g) for 30 min (Eppendorf Centrifuge 5810 R), and the pore water was filtered (0.45-μm surfactant-free cellulose acetate filters, MinisartH) and immediately frozen at −80°C. Pore-water nutrients were analyzed for ammonium and nitrate + nitrite according to standard colorimetric procedures (33). R. maritima was rinsed, dried, and weighed (g dry weight m−2). Microalgal pigments (chlorophyll b; proxy for green algae, echinenone; proxy for cyanobacteria and fucoxanthin; proxy for diatoms) were analyzed using high-performance liquid chromatography according to the method of Wright and Jeffrey (34) (all environmental data are shown in tables S13 to S15).

Habitat diversity

The habitat diversity was used both as a categorical (levels 1 to 4) and as a calculated continuous variable. For the latter, the within- and between-habitat variability can be taken into consideration. The 10 habitat descriptor variables were z-transformed, and multivariate Euclidean distances were calculated for all pairwise habitat combinations and among the replicates within each habitat. The habitat diversity of a given ecosystem was defined as the sum of all pairwise distances between the constituting habitats. The resulting values for habitat diversity were scaled between 0 and 1 by dividing all values (across habitat diversity levels and seasons) by the maximum total Euclidean distance, which was the summer experiment with diversity level 4. Differences in habitat diversity among ecosystems were statistically assessed using multivariate homogeneity of groups’ dispersion (function “betadisper” within the “vegan” package in R) (35). To determine whether there were significant differences among the habitat types, we used PERMANOVA implemented in the “adonis” function within the vegan package in R (35). The clustering of the samples was displayed with Euclidean distances in a PCA plot (Fig. 2A).

Ecosystem functions

We measured four functions—GPP, nitrogen fixation, uptake of DIN, and denitrification—using fresh sediment samples. All sediment–water fluxes (oxygen and uptake of DIN) were measured during day and night by start-stop batch incubations on the ecosystem level and were calculated as daily fluxes. To calculate GPP, light and dark fluxes of oxygen were measured with Unisense oxygen optodes to estimate net primary production and community respiration, respectively. DIN concentrations were assessed from concentrations of ammonium and nitrate + nitrite. For potential denitrification, 2 ml of homogenized sediment from all replicates of each habitat type was added to 10-ml gas-tight glass Exetainers (Labco Limited) flushed with helium and preincubated for 24 hours at in situ temperature to remove eventual oxygen and nitrate present in the sediment. Eighty-milliliter stock solution of 15NO3 (Na15NO3, 99.6 atom %, Europa Science Ltd.) was added to each Exetainer, resulting in a concentration of 50 μM in the sediment pore water. Samples were shaken and incubated in the dark for 2.5 hours. The incubation was terminated by the addition of 6.1 M ZnCl2. Sediment samples were analyzed at the stable isotope laboratory at the University of California, Davis, and potential denitrification was calculated according to Thamdrup and Dalsgaard (36).

Rates of nitrogen fixation were measured using the acetylene reduction assay (37) modified according to Capone (38). Sediment from each habitat type was homogenized, and four subsamples (3 ml each) were added to four 5.8-ml cylindrical glass Exetainers fitted with gas-tight rubber septa (Labco Limited). In each Exetainer, 1.5 ml of filtered seawater was added, leaving approximately 1.3 ml of untreated air as headspace. The headspace was then enriched with ~20% (v/v) acetylene (C2H2) gas. Incubations were terminated after 0, 24, 48, and 72 hours by addition of 0.1 ml of 6.1 M ZnCl2. Nitrogen fixation was assessed by extracting 250 μl of headspace and analyzing the concentration of ethylene using a gas chromatograph (Hewlett Packard 5890, series I). The ethylene production was recalculated to atmospheric N fixation per square meter and day using a ratio of 1 N2 fixed per three ethylene molecules produced (38, 39).

Bacterial community diversity metric and structure

From the cores in the experiments, one sediment sample (~5 g) was taken for each unique habitat such that the number of samples corresponded to the habitat diversity level (1 to 4). For sediment cores of the same habitat within one ecosystem, the samples were pooled. The sediment samples were stored at −20°C. For each sample, DNA was extracted from 0.3 g of sediment using the FastDNA SPIN Kit for Soil (MP Biomedicals) according to the manufacturer’s instructions. The amount of extracted DNA was quantified using the Qubit fluorometer and reagents (Life Technologies Corporation). Bacterial community structure and diversity were assessed by amplification and sequencing of the V3-V4 region within the 16S ribosomal RNA (rRNA) gene using a two-step protocol (40) and the universal primers pro341F and pro805R (41) with Nextera adapter sequences (Illumina Inc.). A total of 33 cycles were used (25 + 8), and the annealing temperature in both polymerase chain reaction cycles was 55°C. Sequencing was performed by Microsynth AG on the MiSeq platform (Illumina) using 2 × 300–base pair (bp) paired-end chemistry.

The resulting sequences were trimmed using FASTX-Toolkit (, and paired-end sequences were merged using PEAR (paired-end read merger) (42) with a minimum overlap of 30 bases, quality score threshold of 30 for the two consecutive bases, and a minimum length of 300 bp for the assembled sequences. Quality filtering was performed with USEARCH version 8.0 (Drive5) (43) using a maximum expected error (“maxee”) value of 1. The final set of sequences was then assigned to OTUs at the 97% sequence similarity level using the USEARCH61 algorithm within QIIME version 1.8.0 (44). We aligned the representative sequences of all OTUs in QIIME using the PyNAST algorithm and the prealigned Greengenes database (version 13_8) as template (44, 45). The phylogenetic tree was built with FastTree and midpoint rooting and made ultrametric using the program PATHd8 (46).

We calculated bacterial diversity for each experimental ecosystem as effective phylogenetic diversity of order q = 1, following the method developed by Marcon and Hérault (47) based on the study by Chao et al. (48). This metric corresponds to the number of species in an equally diverse community, where all species have the same abundance and are equally related to each other. To calculate ecosystem-wide bacterial diversity, the reads of the single habitats within each ecosystem were first rarefied and subsequently summed. At habitat diversity level 3, where one habitat was present twice, the weighted sum was calculated. We excluded five ecosystems, with missing samples. Bacterial community dissimilarity was calculated as abundance-weighted UniFrac distance between the communities, based on OTU data. All data manipulation was performed in R (Development Core Team R) (49), and bacterial diversity was calculated with the entropart package (47).

Benthic microalgal diversity metrics

From the cores in the experiments, sediment samples were taken from the top 5-mm sediment using a 2-ml cutoff syringe. Live cells (with fluorescing chlorophyll) were counted in a Bürker counting chamber using epifluorescence microscopy at ×500 magnification. Cells were identified to the nearest taxon level possible (species or genus), measured, and allocated to size groups. We calculated diversity as the effective number of taxonomical units of order q = 1 (50).

Statistical analysis and calculations of multifunctionality

We used three methods to describe multifunctionality: (i) a weighted average function approach (51), (ii) the threshold approach, and (iii) the multiple-threshold approach (52). Before multifunctionality was calculated, the correlation between the individual functions was also calculated to ensure that none of the functions demonstrates strong correlation (figs. S10 to S12 and tables S16 to S18). All functions were recalculated and standardized to be positive by adding the lowest value to all data and then standardized by the maximum observed value. The weighted average index was calculated by taking the average of all functions and subtracting the SD (51). The relationship between habitat diversity, bacterial diversity, microalgal diversity, and the weighted average index (as well as all individual functions) was fitted using a linear model, with multifunctionality as the dependent variable, season as a categorical independent variable, habitat diversity as an independent continuous variable, and the interaction between season and habitat diversity. For the threshold approach, we selected thresholds of 20, 40, 60, and 80% of the observed maximum multifunctionality and plotted these values against the habitat diversity (52). The multiple-threshold approach created an index of the number of functions surpassing the full range of functional thresholds (0 to 100% of maximum functional values). After identifying the maximum value for each function, we tallied the number of functions above these thresholds for habitat diversity (52).

We also compared the observed multifunctionality and the individual functions in the communities consisting of four habitat types with what we would expect on the basis of the multifunctionality and individual functions in each single habitat type. This value is commonly referred to as the net diversity effect (Divnet) (13). It is defined as the difference between the observed effect in a mixture and the expected effect (Divnet = Divobs − Divexp), where Divobs is the observed multifunctionality/individual function (MF) at the focal diversity level (in our case, the community with all four habitat types) and Divexp is the expected multifunctionality/individual function in the focal diversity treatment based on the single habitat communities (cyano, ruppia, sand, silt). With a substitutive design like ours, Divexp = 0.25 × (MFcyano + MFruppia + MFsand + MFsilt).

Structural equation modeling

To test the relative contribution of habitat and species diversity on ecosystem multifunctionality, we analyzed our data within a SEM framework. We only used bacterial diversity because we did not have sufficient amount of data on benthic microalgal diversity to run a SEM. First, data were separated into three groups—spring, summer, and autumn—and analyzed with a multigroup SEM (52), because there were seasonal differences in the effects of habitat diversity on bacterial diversity and multifunctionality. Second, we performed a comparative fit evaluation between models with or without a direct path between habitat diversity and multifunctionality. The difference in Akaike information criterion (59.978 for the model with no direct path to multifunctionality and 54.000 for the model with a direct path to multifunctionality) indicated that the fully mediated model with a direct path to multifunctionality was the best model for further analysis. Because we have little doubt about the causal structure in the model, the evaluation of our SEM is a strictly confirmatory analysis, meaning that data are compared to only a single model and no model statistics (χ2, df, and P value) are available because the model is fully saturated. Significance levels for individual paths between variables were set at α = 0.05. SEMs were run in Amos (version 20).


Supplementary material for this article is available at

fig. S1. Phylogenetic bacterial diversity.

fig. S2. Benthic microalgal diversity.

fig. S3. Linear model of benthic microalgal diversity and habitat diversity.

fig. S4. Multifunctionality and multiple thresholds.

fig. S5. Ecosystem functions and multifunctionality for individual habitats and habitat diversity.

fig. S6. Temperature and light.

fig. S7. Concentrations of inorganic nutrients.

fig. S8. Total nitrogen and organic content.

fig. S9. Schematic figure illustrating the experimental design.

fig. S10. Ecosystem functions that were used to calculate the multifunctionality index plotted against each other during spring.

fig. S11. Ecosystem functions that were used to calculate the multifunctionality index plotted against each other during summer.

fig. S12. Ecosystem functions that were used to calculate the multifunctionality index plotted against each other during autumn.

table S1. Linear model of bacterial diversity − habitat diversity × season.

table S2. Linear model of microalgal diversity − habitat diversity × season.

table S3. Linear model of multifunctionality − habitat diversity × season.

table S4. Linear model of multifunctionality − bacterial diversity × season.

table S5. Linear model of multifunctionality − algal diversity × season.

table S6. Standardized total, direct, and indirect effects for the group spring.

table S7. Standardized total, direct, and indirect effects for the group summer.

table S8. Standardized total, direct, and indirect effects for the group autumn.

table S9. Linear model of GPP − habitat diversity × season.

table S10. Linear model of N2 fixation − habitat diversity × season.

table S11. Linear model of DIN uptake − habitat diversity × season.

table S12. Linear model of denitrification − habitat diversity × season.

table S13. Environmental data for each habitat during spring.

table S14. Environmental data for each habitat during summer.

table S15. Environmental data for each habitat during autumn.

table S16. Correlation coefficients for the individual ecosystem functions used to calculate multifunctionality during spring.

table S17. Correlation coefficients for the individual ecosystem functions used to calculate multifunctionality during summer.

table S18. Correlation coefficients for the individual ecosystem functions used to calculate multifunctionality during autumn.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank C. M. Jones for helping with the initial analyses of the 16S rRNA gene sequences and F. Käll, S. Tytor, M. Hedblom, O. Bäckman, L. Svanberg, and C. Alsterberg for their contributions during the experiments. Funding: Funding was provided by Swedish Research Council Formas (grant 2012-695) and the project COCOA funded by BONUS (Art 185), funded jointly from the European Union’s Seventh Framework Programme for research, technological development, and demonstration and from Baltic Sea national funding institutions. Author contributions: C.A., K.S., S. Hulth, S. Hallin, and L.G. designed research; C.A., L.G., and F.R. analyzed data; C.A. wrote the first manuscript draft; and all authors performed research and wrote the paper. Competing interests: The authors declare that they have no competing interests. Data and materials availability: Raw data are available through figshare ( and sequencing data through

Stay Connected to Science Advances

Navigate This Article