Does the use of biological traits predict a smooth landscape of ecosystem functioning?

Abstract The biodiversity crisis has increased interest in understanding the role of biodiversity for ecosystem functioning. Functional traits are often used to infer ecosystem functions to increase our understanding of these relationships over larger spatial scales. The links between specific traits and ecosystem functioning are, however, not always well established. We investigated how the choice of analyzing either individual species, selected modalities, or trait combinations affected the spatial patterns observed on a sandflat and how this was related to the natural variability in ecosystem functioning. A large dataset of 400 benthic macrofauna samples was used to explore distribution patterns. We hypothesized that (1) if multiple species (redundancy) represent a trait combination or a modality their spatial patterns would be smoothed out, and (2) the lost spatial variability within a trait combination or modality, due to the smoothing effect, would potentially affect their utility for predicting ecosystem functioning (tested on a dataset of 24 samples). We predicted that species would show heterogeneous small spatial patterns, while modalities and trait combinations would show larger and more homogeneous patterns because they would represent a collection of many distributions. If modalities and trait combinations are better predictors of ecosystem functioning than species, then the smoother spatial patterns of modalities and trait combinations would result in a more homogeneous landscape of ecosystem function and the number of species exhibiting specific traits would provide functional redundancy. Our results showed some smoothing of spatial patterns progressing from species through modalities to trait combinations, but generally spatial patterns reflected a few dominant key species. Moreover, some individual modalities and species explained more or equal proportions of the variance in the ecosystem functioning than the combined traits. The findings thus suggest that only some spatial variability is lost when species are combined into modalities and trait combinations and that a homogeneous landscape of ecosystem function is not likely.

well established. We investigated how the choice of analyzing either individual species, selected modalities, or trait combinations affected the spatial patterns observed on a sandflat and how this was related to the natural variability in ecosystem functioning. A large dataset of 400 benthic macrofauna samples was used to explore distribution patterns. We hypothesized that (1) if multiple species (redundancy) represent a trait combination or a modality their spatial patterns would be smoothed out, and (2) the lost spatial variability within a trait combination or modality, due to the smoothing effect, would potentially affect their utility for predicting ecosystem functioning (tested on a dataset of 24 samples). We predicted that species would show heterogeneous small spatial patterns, while modalities and trait combinations would show larger and more homogeneous patterns because they would represent a collection of many distributions. If modalities and trait combinations are better predictors of ecosystem functioning than species, then the smoother spatial patterns of modalities and trait combinations would result in a more homogeneous landscape of ecosystem function and the number of species exhibiting specific traits would provide functional redundancy. Our results showed some smoothing of spatial patterns progressing from species through modalities to trait combinations, but generally spatial patterns reflected a few dominant key species. Moreover, some individual modalities and species explained more or equal proportions of the variance in the ecosystem functioning than the combined traits. The findings thus suggest that only some spatial variability is lost when species are combined into modalities and trait combinations and that a homogeneous landscape of ecosystem function is not likely.

K E Y W O R D S
benthic macrofauna, biodiversity-ecosystem functioning, functional traits, key species, redundancy, spatial patterns

| INTRODUC TI ON
Environmental change and biodiversity loss have increased interest in the role of biodiversity for ecosystem functioning. The overall goal of biodiversity-ecosystem function (BEF) research has been to understand why biodiversity matters and how ecosystems may be able to maintain the functions that support ecosystem services despite environmental degradation (Srivastava & Vellend, 2005). Evidence of positive relationships between biodiversity and functioning is piling up, but the patterns are often context-dependent and equivocal (Gamfeldt et al., 2015;Reiss, Bridle, Montoya, & Woodward, 2009).
The question has, however, moved on from whether biodiversity has an effect, to how BEF relationships change in space, time, or under specific environmental conditions. There has also been a demand for studies on larger scales in order to support ecosystem-based management decisions (Gamfeldt et al., 2015;Reiss et al., 2009;Snelgrove, Thrush, Wall, & Norkko, 2014). In order to enhance the understanding of BEF relationships, studies have not only started to consider multiple functions (Byrnes et al., 2014;Gagic et al., 2015;Hiddink, Davies, Perkins, Machairopoulou, & Neill, 2009;Villnäs et al., 2013), but, more importantly, the scope of biodiversity components analyzed has been widened (Purvis & Hector, 2000;Reiss et al., 2009;Thrush et al., 2017).
The use of functional traits-defined as any morphological, physiological, phenological, or behavioral feature that can be measured on the level of an individual and can be used to describe its performance (Violle et al., 2007)-has advanced our understanding of BEF relationships. While species richness as a biodiversity metric assumes that all species are equal with respect to function, functional traits demonstrate that species differ in their roles in ecosystem functioning (Bengtsson, 1998;Walker, 1992). If many species in a community express the same traits (redundancy), they might be complementary and support continued function even if a species is lost (i.e., the insurance hypothesis, Yachi & Loreau, 1999). Functional redundancy is thus an important aspect of resilience of ecosystems (e.g., Walker, 1992).
Redundancy in a functional group may also depend on the number and specificity of the traits that are used to form the group, as this will define the number of species that contribute (Micheli & Halpern, 2005). Thus, spatial scales of heterogeneity in the distribution of species with functionally similar traits become important. The variability in species abundance and occurrence within a functional group can be high in heterogeneous landscapes (Hewitt, Thrush, & Dayton, 2008;Walker, 1992;Wellnitz & Poff, 2001), with redundancy within functional groups affected by the spatial variation in species composition (Naeem, Duffy, & Zavaleta, 2012). Additionally, a functional group containing low species richness and low abundance would not necessarily be considered to contain redundancy, but if the group occurs widely over a landscape, it might still be important for the system's resilience (Greenfield, Kraan, Pilditch, & Thrush, 2016).
Scale is thus important in all aspects of biodiversity-ecosystem function relationships. The relationship between biodiversity and ecosystem functioning depends on the investigated temporal (Cardinale et al., 2007;Stachowicz, Graham, Bracken, & Szoboszlai, 2008;Tilman et al., 2001) and spatial scales of heterogeneity (Cardinale et al., 2011;Dimitrakopoulos & Schmid, 2004;Dyson et al., 2007;Godbold, Bulling, & Solan, 2011;Griffin et al., 2009;Raffaelli, 2006). Such strong scale-dependence makes it challenging to validate how well functional traits work as surrogates for functioning, and how the heterogeneity of landscapes is reflected in terms of functional traits.
Estuarine and coastal marine benthic ecosystems are very suitable for BEF studies due to their high biodiversity, encompassing trophic levels, ease of sampling, and the naturally occurring environmental gradients and diverse set of habitats within a relatively restricted distance (e.g., Snelgrove et al., 2014). We used a large dataset of 400 benthic macrofauna samples from an extensive sandflat to explore spatial patterns of species distributions, trait modalities, and trait combinations. From a smaller dataset of 24 samples in a subset of locations, we investigated the ability of individual species, trait modalities, and trait combinations to explain the ecosystem multi-functionality related to nutrient recycling. In this study, the investigated trait combinations were selected based on their importance for nutrient recycling processes at the sediment-water interface (Kristensen et al., 2012;Norkko, Villnäs, Norkko, Valanko, & Pilditch, 2013;Solan et al., 2004), which were investigated on the sandflat through solute flux measurements (Thrush et al., 2017).
More specifically, we investigated how the choice of analyzing either individual species, selected modalities, or trait combinations affected the spatial patterns observed on the sandflat and how this was related to the natural variability in ecosystem functioning.
We hypothesized that (1) if multiple species represent a trait combination or a modality, their individual spatial patterns would be smoothed over. That is, the species would show heterogeneous small spatial patterns across the sampled area, while the modalities and the trait combinations would show larger and more homogeneous patterns because they would represent a collection of many partially overlapping distributions ( Figure 1). We further hypothesized (2) that this lost spatial variability within a trait combination or modality would not affect their ability to predict ecosystem functioning, that is, they would be better predictors than individual species. If both hypotheses were found to be true, then the smoother spatial patterns of the modalities and trait combinations would result in a smoother landscape of ecosystem function and the number of species exhibiting traits would provide functional redundancy.
From the extensive survey data, 24 experimental locations with varying abundance and richness of macrofauna species with functional traits likely to affect nutrient processes in the sediments were selected. These locations correspond to the control plots in the experiments described in Thrush et al. (2017) and shown in their Appendix S1: Figure S1. Solute fluxes (oxygen, ammonium, and phosphate) across the sediment-water interface were measured at each location in March 2014 through chamber incubations and a multivariate response matrix of the normalized solute fluxes were used as a measure of ecosystem multifunction. At the same time, the macrobenthic community was sampled after the flux measurements had been made, using 2 replicate cores (13 cm diameter, 15 cm deep) at each location. The macrofauna cores were sieved (500 µm mesh) and preserved in 50%-70% isopropyl alcohol, before being identified to the lowest taxonomic level possible and counted.

| Flux measurements
Chamber bases (50 × 50 × 10 cm height) were pressed approximately 5 cm down into the sediment during low tide, and when the water level during incoming tide reached approximately 30 cm, the chambers were sealed with Perspex domes. Opaque shade cloths were used to ensure darkness within the chambers throughout the incubation period, which occurred during a midday high tide of approximately 4 hr. Dark incubations were used to control for photosynthetic activities and nutrient uptake by microphytobenthos and seagrass. To measure solute concentrations, water samples (60 ml) from the chambers were collected through sampling ports at the be-

| Trait data
We focused on functional traits that are known to have an impact on nutrient cycling and fluxes at the sediment-water interface by using traits related to vertical movement of sediment particles and solutes, creation of sediment topographic features, body size, and degree of motility (Norkko et al., , 2015Solan et al., 2004;Villnäs et al., 2013). These traits were expected to affect solute fluxes in the sediment by moving sediment particles and organic material, pumping pore water, and changing sediment topography (Thrush et al., 2017;Volkenborn et al., 2012;Woodin et al., 2016).
Where species exhibited attributes of several modalities within one trait, fuzzy probabilities were used to assign species to modalities (Chevenet, Dolédec, & Chessel, 1994), with allocation across traits summing to 1. First, four trait combinations were created from a fuzzy coded modality by species matrix. Three different ways of producing a vertical mixing combination trait were explored and one for surface modification (Appendix S1). The trait combination mixing contained the sum of the modalities surface to depth and depth to surface, describing the direction of particle and solute movement. The trait mixing*L contained the same modalities but only the large taxa (species > 20 mm), and the trait mixing*size*motility contained the same modalities but weighted by size (1 = small <5 mm, 2 = medium 5-20 mm, 3 = large >20 mm) and motility (1 = sedentary/movement in a fixed tube, 2 = limited movement, 3 = freely motile). Body size, F I G U R E 1 Schematic illustration of expected spatial patterns in case of smoothing or no smoothing effects. Each colour depicts a species, modality, or functional group respective to each level, whereas the shades illustrate the abundance of each species, modality, and functional group at each level. A smoothing effect would be indicated by larger homogeneous patches at the modality level and subsequently lead to a homogeneous distribution of the functional group. Whereas, the nonsmoothed patterns would stay more heterogeneous with smaller patches and variable abundances on the modality level, and the distribution of the functional group would be more variable across the seafloor in addition to direction of particle movement, was included since size has been shown to play an important role for nutrient dynamics at the sediment-water interface (e.g., Norkko et al., 2013;Sandnes, Forbes, Hansen, Sandnes, & Rygg, 2000;Thrush, Hewitt, Gibbs, Lundquist, & Norkko, 2006). Larger organisms are for example likely to move more sediment and water, thus creating stronger gradients in pore water pressures within sediments (Volkenborn et al., 2012).
The trait surface modification, describing fauna-produced structures in or on the sediment, contained the modalities permanent burrow, tube structure, simple hole or pit, mound, and trough. After every species was assigned a trait value, the values were abundance weighted and a sum across species was calculated to result in a single value for each trait combination and included modality in each sample (see Appendix S2 for species allocations). The abundance weighted traits were calculated in the same way for both datasets.

| Trait combination selection
Initially The two trait combinations (mixing*L and surface modification) that best explained multifunction (  (Dormann et al., 2007;Legendre & Fortin, 1989), and values range from +1 (strong spatial correlation), 0 (no correlation), to −1 (negative correlation). The analyses were run with equal distance in each distance class, that is, average correlations between samples laying within 0-20 m of each other, 20-40 m of each other, and so on were determined. Before each correlogram of the Moran's I coefficients were explored, a global test was performed to account for the multiple tests conducted. The global tests were performed through checking that each correlogram contained at least one value that was significant at a level of aʹ = a/v, where a was 0.05 and v was the number of tests performed (Oden, 1984).
Autocorrelograms provide information on average spatial patterns and do not necessarily indicate similar specific spatial location (Sokal & Oden, 1978). Therefore, to further investigate the

| Trait selection
The

| Hypothesis 1-spatial patterns of the trait combinations, modalities, and species (autocorrelation)
The correlograms of the Moran's I coefficients indicated that the two investigated trait combinations, the modalities and species, showed

| Hypothesis 1-spatial correlations of trait combinations, modalities, and species abundances
Spearman correlations were used to investigate the spatial location of the patches in relation to each other across the sandflat. The modalities surface to depth and depth to surface, contributing to the trait mixing*L, were highly correlated to the trait combination (i.e., Spearman rho 0.74 and 0.95, respectively; Table 2). Furthermore, the modalities also exhibited a correlation of > 0.5 (0.59), indicating these patches overlapped even if they showed slightly different spatial patterns (see also Figure 2). Thus, these modalities contributed to the trait combination mixing*L in an additive manner.
The species distributions were not strongly correlated (all < 0.2; liliana was also correlated (0.66) with the trait mixing*L, which again together with the matching scales of the patches, suggested that this species could be a driver of the distribution of the trait mixing*L.
There were no strong correlations between the modalities contributing to the trait combination surface modification, and no modality was strongly correlated with the trait combination (Table 3) There were no strong correlations among the species distributions belonging to the trait surface modification (all < 0.45; Appendix S5), indicating that the species were scattered across the sandflat. No species was directly correlated with the trait surface modification, and correlations > 0.5 only occurred between the modality tube structure F I G U R E 3 Maps of the spatial distribution (on the left) and correlograms (on the right) of (a) the trait combination surface modification, the modalities (b) permanent burrow, (c) tube structure, (d) simple hole or pit, (e) mound, and (f) trough. The data in the maps are normalized to run from 0 to 1, to illustrate the patch patterns and locations of the patches across the sandflat (1 km × 0.3 km).

| Hypothesis 2-predictions of multifunction by trait, modalities, and species
Marginal and sequential tests in the DistLM analysis were used to investigate how useful trait combinations, modalities, and species were in predicting multifunction, in order to further elucidate if smoothing of the spatial patterns on the trait level would result in decreased ability to predict functionality. The results indicated that some modalities and individual species accounted for a larger proportion of the variance in the multifunction than the actual trait combinations mixing*L and surface modification (Table 4). For instance, the abundance of the bivalve Austrovenus stutchburyi (26%) could alone explain a larger proportion of the variance in the multifunction than either of the two trait combinations (cf. 23% and 12%). Three of the five modalities included in surface modification also explained more variance than the trait combination surface modification (see Table 4; simple hole or pit, tube structure, trough).

| D ISCUSS I ON
Through exploring spatial distributions of two trait combinations, the included modalities and species, we investigated the potential for functional redundancy to occur over space, and additionally, the degree to which variability in functional traits may predict spatial variability in ecosystem multifunctioning across an extensive sandflat area. We found some evidence for our first hypothesis; modalities and trait combinations would show larger and more homogeneous spatial patterns than individual species because they would represent a collection of many distributions, allowing functional redundancy to occur. However, we did not find clear evidence for our second hypothesis; that this loss of spatial variability would not prevent the trait combinations predicting ecosystem functioning better than individual species.
We chose to investigate two trait combinations that have been proven to be important for nutrient recycling processes (measured as solute fluxes across the sediment-water interface). One trait Instead of finding a relationship between the amount of smoothing and number of species (redundancy) contributing to a modality or trait combination, as was hypothesized (1), we observed a reliance on the dominance patterns of the species included in each trait combination or modality. For example, the trait combination mixing*L displayed less patchiness across the landscape than the trait combination surface modification (Figures 2 and 3), despite higher number of species within surface modification. This was likely due to some overlap between the two modalities within mixing*L and that they had many species in common of the few species in total. The two most abundant species, Austrovenus stutchburyi and Macomona liliana, drove the spatial patterns of the modalities, and thus, also the spatial pattern of the trait combination mixing*L. Similar indications could also be seen in some modalities included in the trait combination surface modification. In two of the modalities (mound and trough), the smoothed patterns (indicated by large-patch patterns, Figure 3 and Appendix S3) were driven by dominant species (Orbinia papillosa and M. liliana, respectively), whereas in two other modalities (permanent burrow and tube structure), all species expressed a high patchiness and many species were included and/or the species had more even abundances, which then resulted in patchy spatial patterns in  (Figure 3 and Appendix S3). Thus, the results indicated that if a modality includes an abundant species with a high contribution to the modality and which is distributed in large patches, this species will drive the distribution of the modality and mask the patchiness of the other less abundant species (i.e., smoothed distribution). Whereas, if a modality includes many species with smaller patch patterns and even abundances, the modality will also show a more heterogeneous distribution across the seafloor.
The ability of a few species and individual modalities to explain more of the variability in our multifunction than the trait combinations (Table 4) was most likely due to the favorable combination of traits some species express for the explored multifunction (particularly Austrovenus stutchburyi in this case). Thus, unlike hypothesized (2), the ability to predict the ecosystem functioning was actually smaller when using the larger-scale, more homogeneous patterns provided by modalities and trait combinations, TA B L E 4 Results from the marginal and sequential tests in DistLM analyses, marginal tests indicate the proportion of the variance each predictor accounts for in terms of multifunction, that is, the phosphate (PO 4 3− ), ammonium (NH 4 + ), and oxygen (O 2 ) fluxes (n = 24). Sequential tests indicate which combination of predictors within the different levels (trait combination, modalities, and species) best explains the variance of the multifunction and stronger relationships were found when using the presence of the individual dominant species. The use of trait combinations or functional groups has a relatively long tradition within the field of ecology (e.g., Fauchald & Jumars, 1979;Wilson, 1999), but the concept and the methods to form these groups are still being developed (Butterfield & Suding, 2013;Murray, Douglas, & Solan, 2014;Törnroos & Bonsdorff, 2012). Earlier studies often based the classification on trophic groups (e.g., Cardinale et al., 2006;Hairston, Smith, & Slobodkin, 1960;Hunt, 1925), whereas today, the classifications are often based on a variety of functional traits that addresses specific research questions (e.g., Harris, Pilditch, Greenfield, Moon, & Kröncke, 2016;Thrush et al., 2017). Functional groups are suggested to be especially useful when describing potential functionality and resilience across larger scales and environmental gradients (e.g., Greenfield et al., 2016;Villnäs, Hewitt, Snickars, Westerbom, & Norkko, 2018). However, in this study, the trait combinations did not work well overall (i.e., low explanation of the multifunction), which might be more common than generally expected. A similar result was reported by Norkko et al. (2015)  Current evidence suggests that even if there is functional redundancy among functionally similar species, the identity of some species may play a pivotal role for the functioning of our ecosystems (e.g., Sandwell, Pilditch, & Lohrer, 2009;Smith & Knapp, 2003). In this study, especially Austrovenus stutchburyi and also Macomona liliana are typically among the most abundant and biomass dominant species living in New Zealand sandflats, and they have been shown in several studies to have an effect on their surrounding environment and ecosystem functioning (Jones, Pilditch, Bruesewitz, & Lohrer, 2011;Pratt, Pilditch, Lohrer, Thrush, & Kraan, 2015;Sandwell et al., 2009;Thrush et al., 2006). Importantly, M. liliana although abundant was not a significant predictor for this particular combination of solute fluxes. The high abundance of this shellfish is driven by high juvenile density, skewing the importance of abundance relative to size and adult living position in the sediment. High abundance of a species does not automatically translate into a large influence on ecosystem functioning. Many studies have, however, demonstrated significant identity effects on functioning (Stachowicz, Bruno, & Duffy, 2007), with the composition of species equally or even more important than species richness (Bruno, Boyer, Duffy, Lee, & Kertesz, 2005;Gammal et al., 2019;Gammal, Norkko, Pilditch, & Norkko, 2017;Gustafsson & Boström, 2009).
Different functions, even if closely related, similar to the functions included in the multifunction of this study, have been indicated to be dependent on different aspects of biodiversity (Murray et al., 2014).
Functional redundancy when multiple functions are considered on larger and longer scales is thus likely much lower than the redundancy for single functions (Gamfeldt et al., 2008;Rosenfeld, 2002), highlighting that ongoing biodiversity loss could have far graver consequences for the overall functioning of ecosystems than earlier presumed (Reiss et al., 2009).

ACK N OWLED G M ENTS
This study was supported by Walter and Andrée de Nottbeck Marine Laboratory for the great hospitality.

CO N FLI C T O F I NTE R E S T
The authors have declared that no conflict of interest exists.