Functional trait diversity maximizes ecosystem multifunctionality.

Understanding the relationship between biodiversity and ecosystem functioning has been a core ecological research topic over the last decades. Although a key hypothesis is that the diversity of functional traits determines ecosystem functioning, we do not know how much trait diversity is needed to maintain multiple ecosystem functions simultaneously (multifunctionality). Here, we uncovered a scaling relationship between the abundance distribution of two key plant functional traits (specific leaf area, maximum plant height) and multifunctionality in 124 dryland plant communities spread over all continents except Antarctica. For each trait, we found a strong empirical relationship between the skewness and the kurtosis of the trait distributions that cannot be explained by chance. This relationship predicted a strikingly high trait diversity within dryland plant communities, which was associated with a local maximization of multifunctionality. Skewness and kurtosis had a much stronger impact on multifunctionality than other important multifunctionality drivers such as species richness and aridity. The scaling relationship identified here quantifies how much trait diversity is required to maximize multifunctionality locally. Trait distributions can be used to predict the functional consequences of biodiversity loss in terrestrial ecosystems.

Global threats to biodiversity have motivated ecologists to better understand the relationship between biodiversity and ecosystem functioning1 (BEF hereafter). Pioneering BEF experiments conducted in grasslands have shown a positive effect of species richness on ecosystem functioning2,3, which has been later confirmed in many ecosystems4–7. However, there is a general consensus that it is not only the number of species per se that influences ecosystem functioning, but also the diversity and the abundance of their functional traits within communities8–10. Because functional traits relate to how species acquire, share, and conserve resources11, they are often invoked to explain how species assemble within communities12–14 and impact ecosystem functioning8–10,15.

Higher trait diversity has been hypothesized to enhance ecosystem functioning because co-occurring species with contrasting trait values may exploit different resources or the same resources at different spatial or temporal scales2, and then may increase overall resource utilization16. This hypothesis has received support from BEF experiments17 and empirical surveys10,18,19. However, most trait-based BEF studies have been conducted at a local or regional scales17,20, limiting our ability to generalize the effect of trait diversity on ecosystem functioning across ecosystems. Also, many studies have focused on single ecosystem functions (e.g. productivity10), something that may provide a partial measure or even biased estimate of overall ecosystem functioning if trade-offs or synergies among functions are ignored4,21,22. Emerging empirical evidence suggests that higher trait diversity is especially important to sustain multiple ecosystem functions simultaneously19,23 (multifunctionality hereafter). However, the extent to which trait diversity matters for multifunctionality remains unknown at a global scale.

Coupling a global survey of multifunctionality conducted in 124 perennial plant dryland communities from all continents except Antarctica5 (Fig. 1) to an extensive trait database24, we uncovered a scaling relationship between the abundance distribution of functional traits within communities and maximized multifunctionality. We quantified the 124 abundance distributions for specific leaf area (SLA) and maximum plant height (trait-abundance distributions hereafter)25. These two traits capture the global spectrum of plant form and function26, and have been largely used to define species resource utilization strategies in terrestrial ecosystems27. Multifunctionality was quantified using five key ecosystem variables: plant productivity, soil enzymatic activities (phosphatase and β-glucosidase), ammonification and potential N transformation rate. These variables were uncorrelated (Supplementary Table 1) and are good proxies for nutrient cycling, biological productivity, and soil fertility in global drylands5,28 (see Methods for details).

Figure 1

Map of the 124 sampled drylands used in this study.

Global map (a) and close-up maps of the sites located in North America (b), the Mediterranean Basin (c), China (d), South America (e), Africa and the Middle East (f) and Australia (g).

We developed an analysis, inspired by multivariate optimization procedures that are used in various scientific fields (e.g. physics, geomorphology, climatology and particularly in economy, where it is employed for analysing budget portfolios maximizing benefits and minimizing costs in multiple investments29; see Supplementary Note 1 for more details), to link the trait-abundance distributions to multifunctionality. Specifically, our approach focused on the skewness and the kurtosis of the trait-abundance distributions24,25,30–32. Whereas the mean and the variance reflect the location and the scale of a distribution (i.e., the dispersion of trait values within a community), the skewness and the kurtosis describe its shape. The skewness and the kurtosis complement the mean and the variance, particularly when assessing complex distributions that strongly deviate from the normal distribution24,25,32(Fig. 2).

Figure 2

The Skewness-Kurtosis Relationship (SKR) represented in the Cullen and Frey graph37.

This graph represents the location of common distributions in the skewness-kurtosis space. The normal distribution represents a unique combination of skewness and kurtosis values of 3 and 0, respectively (see other example in the right panels). Families of trait distributions can be represented by a line for the gamma or lognormal distributions, and by a surface for the beta distribution. Note that the normal distribution represents a special case of the gamma and the lognormal distributions. The skewness-kurtosis relationship has a lower boundary below which no kurtosis value can be observed for any degree of skewness. This lower boundary corresponds to the Bernoulli distribution. We also indicate the location of the empirical SKRs for specific leaf area and maximum plant height (green and red dashed lines, respectively, see Main Text). Low values of kurtosis indicate high trait diversity25.

The skewness represents the asymmetry of the distributions. High negative or positive values of skewness occur when trait-abundance distributions are strongly left-, or right-tailed, with a few abundant species that have extreme trait values compared to the bulk of the distribution. Skewed trait distributions typically result from phenomena such as asymmetric competition33,34 or rapid environmental change25,32. Kurtosis represents the relative peakiness of the trait-abundance distribution and the heaviness of its tails, and is a measure of trait diversity25,30,31. High kurtosis values characterize peaked distributions, indicating a high abundance of species with similar trait values in a given community. High kurtosis values therefore reflect a low trait diversity25, which might occur with strong environmental filtering32. Low kurtosis values reflect an even abundance distribution of trait values within the community, which implies the co-dominance of functionally contrasted species32. Low kurtosis reflects a high trait diversity25, for example in case of limiting similarity35 or storage effects36.

By analysing the shape of the 124 trait-abundance distributions, we first showed that trait diversity in dryland plant communities is generally higher than expected by chance. Then we showed that observed trait distributions in global drylands maximize multifunctionality, in agreement with the hypothesis that species with contrasting trait values collectively exploit a greater diversity of resources2,16.

Results and Discussion

The Skewness-Kurtosis relationship

Skewness (S) and kurtosis (K) are mathematically related according to the following Skewness-Kurtosis Relationship37, SKR hereafter (Fig. 2): This inequality generates a mathematical constraint triangle on the possible values of skewness and kurtosis that can be observed for a community (Fig. 2, see also Supplementary Note 1 for detailed information). For a given set of distributions, an increase in (skewness)2 will generate an increase in kurtosis all other things being equal. In other words, as the trait-abundance distributions become more skewed, they also become more peaked, indicating a decrease in trait diversity25. The slope β of the SKR measures the strength of the relationship, i.e., the extent to which trait diversity decreases as trait distributions become more skewed. The Y-intercept α indicates the lowest kurtosis value at skewness = 0, which corresponds to the highest trait diversity predicted by a given SKR. Importantly, the inequality has a lower boundary that sets a limit to the minimal kurtosis value predicted for any degree of skewness, i.e. the potential maximum trait diversity (K = S2 + 1, Fig. 2). Although skewness and kurtosis separately provide valuable information on how trait diversity and abundance is distributed within communities, the SKR reveals the extent to which trait diversity is maximized.

Observed trait distributions in global drylands

The observed 124 trait-abundance distributions were highly heterogeneous, and deviated from a symmetric normal distribution in most cases (Fig. 3). Consistent with equation (1), we observed a strong positive relationship between skewness (S) and kurtosis (K) of the distributions for both SLA and height (H) (Fig. 3): However, within the constraint triangle, there are a variety of possible SKRs that can be generated (Fig. 2). Empirical SKRs must then be compared to the null distribution of relationships that might arise by sampling random values within the constraint triangle. We used a Monte-Carlo analysis to test whether the empirical SKRs differed from random expectations. These randomization tests used three null models that imposed increasing constraints on the species and trait pools.

Figure 3

Observed skewness-kurtosis relationships (SKRs) and deviation of from null expectations.

Observed and null SKRs for a, specific leaf area (SLA) and b, maximum plant height. Blue dots represent the skewness and kurtosis values of random communities (2000 randomizations × 124 communities); purple lines show the 2000 random SKRs; green and red dots represent the observed 124 communities; the green and red dashed lines represent the observed SKRs for SLA and height. We indicate the equation for each observed SKR, and the conditional pseudo-P values from null model “Richness” for the slope β, P(β|α), the Y-intercept α, P(α|β), the whole model, P(β∩α) and the distance to the lower boundary, P(dis) (see Table 1 and Supplementary Figure 2 for details). Right panels indicate the relationship between the slopes β of the random SKRs and the random Y-intercepts α (null parameters, blue dots); the green and red dots indicate the position of the observed parameters. The observed SKRs for SLA and height had significantly different Y-intercepts (Slope β: Fratio1,244 = 5308.16, P<0.001; Y-intercept α: Fratio1,244 = 16.8, P<0.0001; Slope β×Y-intercept α, Fratio1,244 = 3.68, P= 0.06).

In the first null model, we simulated the assembly of local plant communities from global pools of species and traits. In this null model, the local species number and trait values are randomly assigned to the null communities according to the hypothesis that the local environment may modify both local trait composition and species richness38. In the second null model, we fixed species richness at the site level, but allowed traits to vary across communities to decouple the effect of species richness and trait composition on local trait abundance distributions. In the last null model, we additionally constrained the range of trait value observed within community to be fixed. This third randomization procedure tested whether the distribution of trait values between co-occurring species maximized trait diversity locally (see Methods).

The empirical SKRs for SLA and height deviated strongly from the predictions of all three null models (Fig.3). The slopes β and the Y-intercepts α were all lower than expected by chance (Fig. 3 right panels, Table 1a; see detailed results in Supplementary Note 2, Supplementary Figure 2). Consequently, observed kurtosis values for each of the 124 trait-abundance distributions for both SLA and height were significantly closer than expected by chance to the lower boundary of the mathematical constraint triangle (Table 1b, see the analysis on the distance to the lower boundary in Methods). In other words, observed kurtosis within communities was minimal after controlling for the degree of skewness in the data. Deviations from null expectations were consistent for all three null models (Table 1), so the results are robust to the details of how species and trait pools were constructed and randomized.

Table 1

Results from Null Model for specific leaf area (SLA) and maximum plant height (H).

a, Direct assessments of probabilities derived from the null models for the skewness-kurtosis relationships (SKRs). b, Distance to the lower boundary. In Table 1a, the parameters β and α were tested separately, conditionally, and jointly under the three randomization procedures. P(β): probability of finding a lower random slope β during the randomization procedure than the observed slope β; P(α): probability of finding a lower random Y-intercept α than the observed α; P(β|α): conditional probability of finding a lower random slope than the observed slope; P(α|β): conditional probability of finding a lower intercept in the randomizations than observed; P(β∩α): probability of finding in the randomizations both a lower intercept and slope than observed. We indicate the regression parameters (estimate ± SE, r2) for observed SKRs. Pseudo P-values below 0.05 indicate that the observed parameters are lower than the parameters expected by chance. In Table 1b, we compared the distance of observed data and the predictions of the three Null models for specific leaf area (SLA) and maximum plant height (H). P-values below 0.05 indicates that observed values are closer to the lower boundary of the SKR than expected by chance (see Methods for more information on each null model; Supplementary Information S2).

a. Direct assessments										b. Distance To the Lower Boundary
Trait	Observed coefficient			Null model	P-values					Trait	Null model	probability
			r2		P(β)	P(α)	P(β|α)	P(α|β)	P(β∩α)			Pvalue
SLA	β	1.236 ± 0.049	0.839	Frequency	0.279	0.313	0.002	0.002	0.001	SLA	Frequency	0.121
	α	2.149 ± 0.716		Richness	0.470	0.149	0.013	0.004	0.002		Richness	0.022
				Swap	0.766	0.118	0.047	0.007	0.006		Swap	0.119
log (SLA)	β	1.010 ± 0.019	0.956	Frequency	0.022	0.505	0.001	0.023	0.001	log (SLA)	Frequency	0.023
	α	1.961 ± 0.217		Richness	0.003	0.084	<0.001	<0.001	<0.001		Richness	<0.001
				Swap	0.371	0.196	0.031	0.016	0.006		Swap	0.011
H	β	1.482 ± 0.059	0.835	Frequency	0.711	0.504	0.455	0.322	0.229	H	Frequency	0.604
	α	2.52 ± 1.629		Richness	0.874	0.285	0.574	0.187	0.164		Richness	0.011
				Swap	0.830	0.379	0.562	0.257	0.213		Swap	0.663
log (H)	β	1.065 ± 0.020	0.955	Frequency	0.045	0.406	0.001	0.011	0.001	log (H)	Frequency	0.051
	α	2.921 ± 0.236		Richness	0.035	0.235	<0.001	<0.001	<0.001		Richness	0.004
				Swap	0.344	0.201	0.010	0.006	0.002		Swap	0.024

A general assembly rule in global drylands

Our findings have important implications for understanding community assembly in global drylands. The abundance and the diversity of trait values observed within the 124 communities evaluated cannot be observed by chance, and can be summarized by two empirical SKRs, one for height and another for SLA (Figs. 2, 3). Each of these SKRs defined a family of trait abundance-distributions in which the trait distribution of each studied community represents a particular instance of a more general distribution quantified at the biome scale (Fig. 2). These two biome-scale distributions predict that the kurtosis of trait distributions is always minimized (Fig. 3). Our study therefore uncovered a general “assembly rule”12,39: across a variety of continents, species pools, and environmental conditions, trait diversity is maximized within dryland plant communities.

For asymmetrical distributions (skewness2 > 0, Fig. 2), lower kurtosis than expected by chance implies heavier-tailed distributions, with a higher relative abundance of rare species with extreme trait values. This pattern is consistent with the hypothesis that trait differences between rare and common species promote coexistence trough stabilizing mechanisms of niche differentiation13,14,36, e.g. limiting similarity35 and associated density-dependent mechanisms40 or facilitation41. For symmetrical distributions (skewness = 0, Fig. 2), observed kurtoses were also lower than expected by chance. For the null assemblages, the Y-intercept α at kurtosis ≈ 5 corresponds to a hyperbolic-secant distribution (Fig. 2). In contrast, the observed Y-intercept α at kurtosis ≈ 2 for SLA, and kurtosis ≈ 3 for height correspond to uniform and normal distributions, respectively (Fig. 2). A uniform distribution implies strikingly high levels of trait diversity within dryland communities and the co-dominance of functionally contrasting species.

Trait distributions and ecosystem multifunctionality

The skewness and kurtosis of the trait-abundance distribution were the main biotic drivers of multifunctionality (Fig. 4). The model that included all abiotic and biotic predictors explained up to 66 % (adj. R) of the total variation observed in multifunctionality. Together, the skewness and kurtosis for SLA and plant height explained 38% of the variation in multifunctionality, which was greater than that explained by the mean and variance of these two traits (15%). Species richness explained only 5% of the variation in multifunctionality, which is comparable in magnitude to that reported in earlier studies5 (see Supplementary Note 3 for more details).

Figure 4

Effect of abiotic and biotic factors on ecosystem multifunctionality.

a, average parameter estimates (standardized regression coefficients) of model predictors, associated 95% confidence intervals and relative importance of each factor, expressed as the percentage of explained variance. The adj. R of the averaged model and the p-value of each predictor are given as: *, p < 0.05; ** p < 0.01; *** p < 0.001. See Supplementary Table 3 for model selection. b, Frequency distribution of adj. R from the null model analyses testing the relationship between trait-abundance distributions and multifunctionality. Black line is the 0.95 quantile of the null model prediction. The red line is the observed adj. R. We provide pseudo-P values for each null model (see Supplementary Note 3 for detailed results).

To investigate the BEF relationship, ecologists have usually analysed the effect of biotic attributes such as species richness1–3,5, community-weighted means of the trait-abundance distribution8,11,15,25, or diversity indexes related to trait dispersion within communities (e.g. the variance of trait distribution)10,17,23. Our results do not contradict the findings of previous BEF studies, because we also found that species richness, the mean and variance of the trait distribution were all significant predictors of variation in multifunctionality. However, we showed that models including the skewness/kurtosis of the trait-abundance distributions had a much higher explanatory power for multifunctionality than models including only the mean and the variance of the trait abundance distribution (see Supplementary Note 3). Considering the skewness-kurtosis relationship represents a significant step forward in our ability to predict multifunctionality in global drylands (~ additional 38% explained variance compare to previous studies5).

Linking a general assembly rule to ecosystem multifunctionality

We conducted a sensitivity analysis of the best-fitting model to explore how coordinated changes in the skewness and the kurtosis of the trait-abundance distributions affected multifunctionality (model 4, Fig. 4, see Supplementary Note 4 for more details on the sensitivity analysis). The two empirical SKRs predict a peak of multifunctionality to occur at their Y-intercept α (skewness fixed at = 0), i.e., at the highest trait diversity predicted by the two empirical SKRs (green & red lines in Fig. 5). This peak is predicted to occur at the uniform distribution for SLA (SKRSLA predicts a kurtosis value ≈ 2 at the Y-intercept) and the normal distribution for plant height (SKRH predicts a kurtosis value ≈ 3 at the Y-intercept) described above (Fig. 2). Each empirical SKR also describes a ridge that minimizes the loss of multifunctionality for any change in the degree of skewness (see the green & red lines in Fig. 5). The existence of this ridge implies a tight coupling between the empirical SKRs and a maximization of multifunctionality in drylands.

Figure 5

Scaling relationship between the empirical SKRs and ecosystem multifunctionality.

Predictions of the sensitivity analysis of the best selected multifunctionality model (model 4, Fig. 4, see Supplementary Note 4 for the detailed methods). Multifunctionality was predicted on a skewness axis for SLA (a) and height (b) for: (i) the two empirical SKRs where KSLA = 1.01 SSLA2 + 1.96 (green line on panel a) and KHeight = 1.06 SHeight 2 + 2.92 (red line on panel b); (ii) the lower boundary where KBoundary = SBoundary2 + 1 (grey area: 95% confidence band); and (iii) and the 124 communities studied (green and red dots). We located in the figure the position of emblematic dryland ecosystems from around the globe and dominated by contrasted flora (diamond symbols): a, Northern African steppes dominates by Stipa tenacissima (yellow diamonds) and Mediterranean shrublands dominated by Rosmarinus officinalis, Thymus vulgaris, Quercus coccifera and Q. ilex (orange diamonds); b, American xeric shrublands dominated by Tetracoccus halli, Larrea tridentata and Ambroisia dumosa (blue diamonds) and Australian woodlands dominated by Casuarina cristata, Alectryon oleifolius and Eucalyptus sp. (purple diamonds). Within each ecosystem, changes in the shape of the trait abundance distributions can be used to monitor changes in multifunctionality. Photos credit: David Eldridge, Nicolas Gross, Fernando Maestre.

Using the same randomization test applied across communities (see Methods, Supplementary Note 5), we found that the relationship between trait-abundance distributions observed within communities and multifunctionality also cannot be explained by chance (Fig. 4b). The two empirical SKRs corresponded to the unique sets of trait distributions able to explain multifunctionality with such a high predictive power (Fig. 4b). Results were robust to the three randomization scenarios employed, including the most highly constrained null model based on the local species composition and trait pool observed within communities (Supplementary Note 3). Our analysis thus supports the existence of a unique scaling relationship linking plant functional traits (SLA and plant height), the distribution of their abundance within communities, and multiple ecosystem functions.

It is noteworthy that increasing trait diversity to its highest potential values defined by the lower boundary of the constraint triangle does not necessarily increase multifunctionality (Fig. 5b). For example, decreasing the kurtosis value for height beyond the normal distribution toward a kurtosis value of 1 (i.e., the lower boundary) had a strong detrimental effect on the multifunctionality of xeric shrublands and dry forests (blue and purple diamonds in Fig. 5b, respectively, see also Supplementary Note 4). We observed: (i) a maximized multifunctionality state closer to a normal distribution (y-intercept α ≈ 3, along red line in Fig. 5b); and (ii) a reduced multifunctionality state closer to the lower boundary (see the shaded area in Fig. 5b). The two empirical SKRs not only show that high trait diversity is needed to sustain multifunctionality, but also precisely quantify how much diversity is required to maximize it locally.

Our study suggests that SLA and maximum plant height – two key plant attributes26 – can be used as efficient functional markers (sensu ref.15) to scale up trait diversity to the ecosystem level. However, the mechanistic linkages between these two traits and multifunctionality remain to be explored. In this context, other proximal traits known to impact ecosystem functioning should be considered in future studies, such as below-ground traits42 (e.g. rooting depth), element stoichiometry43 (e.g. C:N ratio), and litter traits44 (e.g. lignin content). Also, intraspecific trait variability has been shown to modulate community assembly45,46 and ecosystem functioning47 and was not included in our study. How intraspecific trait variability modulates multifunctionality at local and global scale remains largely unexplored. Lastly, our study focuses only on perennial plants. The role of functional diversity of annual plants, which can account for a large proportion of plant species diversity in drylands, certainly deserves additional attention.

Conclusion

We uncovered a general family of trait abundance-distributions operating at the biome-scale that quantify the maxima of trait diversity needed to sustain multifunctionality in global drylands (Fig. 5). Biome-scale distributions emerged from local-scale processes of maximization operating within communities, and highlight a tight coupling between community assembly and ecosystem functioning in global drylands. By quantifying a scaling relationship between the abundance distribution of functional traits and maximized multifunctionality, our analysis suggests that trait distributions could be used to assess the functional consequences of biodiversity loss in terrestrial ecosystems and to guide management efforts aimed at maintaining key ecosystem services linked to productivity and soil fertility.

Methods

Characteristics of the study sites

We obtained field data from 124 sites located in 13 countries (Argentina, Australia, Chile, China, Ecuador, Israel, Kenya, Mexico, Morocco, Spain, Tunisia, USA and Venezuela; Fig. 1). These sites are a subset of the global network of sites from ref.5. Our dataset includes representative sites from the major vegetation types found in drylands, which differ widely in plant species richness (5 to 44 species per site, average = 16.44) and environmental conditions32 (mean annual temperature and precipitation ranged from -1.8°C to 27.8°C, and from 79 mm to 1177 mm, respectively).

Trait-abundance distributions

We focused on two key plant functional traits26: specific leaf area (SLA), which indexes leaf-level carbon gain strategies48, and maximum plant height, which reflects a trade-off for biophysical constraints in determining water fluxes within the plant49, and is related to competitive ability33. We quantified the abundance distributions of these traits, which describe the relative abundance of each trait value within each of the 124 communities, by using two independent datasets: (i) a detailed dataset containing the percent cover of each perennial plant species measured in 80 quadrats of 2.25 m2 within each site, where the sum of the cover for each species is used as a proxy of species abundance at each site5; and (ii) data for mean trait value per species for SLA and maximum plant height retrieved from the TRY database24. The number of sites selected was based on the availability of trait data. We selected those sites for which trait data were available for all the most common perennial plant species that together accounted for a cumulative relative abundance >80%50 (averaged cumulative relative abundance for SLA = 94.3 ± 6.5; for plant height = 96.6 ± 5.24). Therefore, our approach focused on the dominant and subordinate species, known to impact the most ecosystem functioning, according to the mass-ratio hypothesis51. We acknowledge that our analyses do not account for the effect of species with very low abundance. In total, trait data were available for 316 and 526 plant species out of 622 species for SLA and maximum plant height respectively.

We calculated the mean, variance, skewness, and kurtosis of the 124 trait-abundance distributions25,32 for SLA and maximum plant height separately as follows: where pi and Ti are the relative abundance and the trait value of the species i respectively, n is the total number of species in a community with available trait values. For each community, the sum of relative abundance is equal to 100%, i.e. Σ = 1. In contrast to the mean and the variance, the skewness and the kurtosis are unitless. The weighted abundance trait distribution quantifies how much relative abundance exhibits each trait value present in a given community, and provides the shape of the distribution (see Supplementary Note 1).

Assessing multifunctionality

We estimated multifunctionality for the 124 selected sites using plant productivity, soil enzymatic activities (i.e., activity of phosphatase and β-glucosidase), ammonification and potential N transformation rate. These variables (hereafter functions) were uncorrelated among themselves (Supplementary Table 1), and measure either “true” ecosystem functions (sensu ref.52; e.g. potential N transformation rate, productivity) or are key properties/processes (sensu ref.53; e.g., soil enzymatic activities), which together constitute a good proxy for nutrient cycling, biological productivity, and buildup of nutrient pools. The ecosystem functions evaluated are also considered to be critical determinants of ecosystem functioning in drylands5,28,43,54,55. We used a conservative approach to quantifying multifunctionality by minimizing the number of ecosystem functions considered. However, including additional ecosystem functions related to carbon, nitrogen and phosphorus stocks and cycling when assessing multifunctionality did not affect our results (see Supplementary Note 3).

Plant productivity was quantified using the Normalized Difference Vegetation Index (NDVI) at the site level, which acts as a proxy of photosynthetic activity of terrestrial ecosystems for a given composite period28,56. These data were obtained with a resolution of 250 m from the Moderate Resolution Imaging Spectroradiometer aboard NASA's Terra satellites (http://daac.ornl.gov/index.shtml). Here, we used NDVI values obtained for the months before, during and after sampling at each of the surveyed plots averaged across a period of 13 years (2000 to 2013). Soil variables were measured from soil samples (0 and 7.5 cm depth) collected under the canopy of the dominant perennial plants, and in open areas devoid of vascular vegetation as described in refs.5,43. The different ecosystem functions measured were used to obtain a multifunctionality index (M), as described by Maestre et al.5. First, we standardized separately the five ecosystem functions (F) using the Z-score transformation: where Fij is the value of an ecosystem function i in the community j , and Mean F and SD F are the mean and the standard deviation of the ecosystem function F calculated for the 124 studied communities, respectively.

We then obtained M for each community j as the average of the Z-scores of the five ecosystem functions i assessed.

M has good statistical properties5,21 and is a straightforward and easily interpretable measure of multifunctionality21. All selected single-functions were positively correlated with M (Supplementary Table 1a).

Environmental and spatial variables measured

We summarized the climatic features of the studied sites using aridity, a major determinant of ecosystem structure and functioning in drylands worldwide54,55. We obtained values of the aridity index (AI, precipitation/potential evapotranspiration) from ref.57, which uses the data interpolations provided by Worldclim58. To facilitate the interpretation of results, we calculated the aridity level for each site as 1 – AI, so higher values of this aridity level indicate drier conditions43. Aridity is strongly correlated with annual mean precipitation in our dataset (r = -0.84).

We summarized local topo-edaphic parameters at each site using slope angle (as a surrogate of topography), soil sand content, and soil pH. These variables, measured as described in ref.5, play key roles in the availability of water and nutrients in drylands57, and are major drivers of the composition and diversity of plant and microbial communities28,32. Clay and silt content were not used in our analyses due to their correlation with sand content (r = -0.52 and -0.55, respectively). We also considered the elevation, latitude, and longitude of the study sites in our analyses to account for potential effects of spatial autocorrelation among them5,32. We used the sinus and cosinus of the longitude to avoid any bias due to intrinsic circularity of longitude in the models.

Community level analyses

The mathematical relationship between the skewness (S) and kurtosis (K) of the trait-abundance distributions (SKR) takes the form of an inequality37 (1), where K = S defines the lower boundary to the SKR (Fig. 2, see Supplementary Note 1 for a full mathematical demonstration of the SKR and for more information on its implications and use in various scientific fields).

We used a null model approach to test whether observed SKRs differ from random expectations. Using the package PICANTE59 in the R statistical software, version 3.1.160, we constructed three null models with increasing constraints on the randomization procedure. We performed 2000 randomizations for each null model. We assumed in the first null model that local communities reflect random distributions of species drawn from the global species pool. In this null model, species abundance was shuffled across communities using the function frequency in PICANTE. This procedure kept the abundance of each species constant at the global scale, but allowed the distributions of species richness and abundance to randomly vary within communities. We assumed in the second null model that these distributions were fixed within communities (function richness in PICANTE). This null model randomized the trait values across species to break out the observed relationship between species traits and abundance. The randomization of trait values avoided any bias in null predictions promoted by local variations in the number and cover of species observed at the community scale61. In the third null model, we restricted species richness and the range of trait values to those already filtered by the environment within communities (function independentswap in PICANTE). Therefore, trait values are randomly shuffled across the species occurring in each community.

We calculated the mean, variance, skewness, and kurtosis of the trait-abundance distributions for SLA and maximum plant height for each of the 2000 randomizations and for each of the null models used. These were calculated using both raw and log-transformed trait values, as log-transformation may impact the level of detection of non-random processes62. We calculated the 2000 random relationships between skewness and kurtosis (SKRrandom) (Eqn 1). Then we extracted both the slope (β) and the Y-intercept (α) of the 2000 SKRrandom. We compared then β and α to observed parameters of the SKR (β and αobs) for each of the three independent null models used. We used two different approaches to assess whether the observed SKR significantly differed from SKRrandom:

: we compared observed parameters of the SKR (i.e., the slope β and Y-intercept α) to those generated by the three null models. Importantly, null models predicted a negative relationship between β and α, indicating a dependency between random parameters (Supplementary Figure 1). Therefore, we used the Bayes theorem to calculate conditional probabilities for β and α (Supplementary Figure 1). Specifically, we calculated:

an initial pseudo P-value for β (P(β)), generated by comparing β with the probability distribution of β (Supplementary Figure 1a). We counted the frequency with which we generated a random slope (β) smaller than the observed slope (β). We then calculated P(β) as the number of randomizations with β ≤ β divided by the total number of randomizations. An observed slope β significantly lower than the probability distribution of β indicates that the slope of the observed SKR is less steep than expected by chance. Any change in the skewness of the distributions leads to slower increased in kurtosis than expected by chance.

a second pseudo P-value for α (P(α)) (Supplementary Figure 1b). For doing so, we counted the frequency of observing a Y-intercept expected by chance (α) being lower than the observed Y-intercept (α). We then calculated P(α) as the number of randomizations with α ≤ α divided by the total number of randomizations. An observed Y-intercept (α) significantly lower than expected by chance indicates a lower kurtosis at the Y-intercept than expected by chance.

a third pseudo-P value for α taking into account the observed slope β (α if β: conditional P-value = P(α|β); Supplementary Figure 1c). We obtained P(α|β) by counting the number of observations with an Y-intercept expected by chance (α) lower than the observed Y-intercept (α), within the subset of randomizations in which β is lower than the slope expected by chance (β).

a fourth pseudo-P value for β taking into account the observed Y-intercept α (β if α: conditional P-value = P(β|α); Supplementary Figure 1d). We calculated P(β|α) by using the Bayes theorem:

a fifth pseudo-P value including β and α together (β+α, P(β∩α); Supplementary Figure 1e). We followed the Bayes equation to calculate P(β∩α): P(β∩α) represents the probability to observe a whole model with both a slope and an Y-intercept (β and α) being lower than expected by chance.

Distance to the lower boundary: we used a complementary approach to assess whether the observed SKR differs from SKRrandom by testing whether observed kurtoses of each of the 124 communities sampled were significantly closer to the lower boundary than the random data. The lower boundary defines the lowest level of kurtosis that can be expected for any degree of skewness (Fig. 2). If the 124 observed communities are closer on average to the lower boundary than expected by chance, we can conclude that trait diversity within communities is maximized for any degree of skewness. In each randomization event, we calculated the distances to the kurtosis at the lower boundary (K) for both random and observed kurtosis values for each community j, where K is defined as: where S is the Skewness square value of the community j. We then calculated the distance (D) of the community j to the lower boundary as: We fitted an ANOVA model using the function aov() in R to compare the observed distances to the lower boundary and those obtained for each randomization. A positive estimate of the ANOVA comparison indicates that the observed kurtosis was closer to the lower boundary than expected by chance. We performed 2000 randomizations for each of the three null models, and stored all resulting estimates from the ANOVA. Using these data, we calculated a pseudo P-value based on the frequency of observing a slope ≤ 0 within the distribution of 2000 slopes.

Ecosystem level analyses

We used multiple regression models to assess the effect of trait-abundance distribution on multifunctionality. We built four competing models including as predictors: (i) spatial (latitude, cos-longitude, sin-longitude, and elevation) and abiotic (aridity, soil sand content, slope angle, and pH considering a quadratic term for pH only) variables (“Abiotic” model); (ii) all variables of the model “Abiotic” plus the mean and the variance of trait-abundance distributions for SLA and plant height (“Biotic Low Moments” model); (iii) all variables of the model “Abiotic” plus the skewness and the kurtosis of trait-abundance distributions for SLA and plant height (“Biotic High Moments” model); (iv) all variables of the model “Abiotic” plus the mean, the variance, the skewness and the kurtosis of trait-abundance distributions for SLA and plant height (“Biotic All Moments” model). We also included in the four models species richness as predictor to control for its potential effects on multifunctionality5. Species richness was not correlated with the four moments of the trait-abundance distributions for both maximum plant height and SLA (Supplementary Table 2) except in the case of the variance, where a positive trend was observed (r < 0.40 for maximum plant height and SLA). We also considered quadratic terms for the mean, the variance and the kurtosis both for height and SLA to assess potential non-linear effects of these variables on multifunctionality23. Note that we never included a quadratic term for skewness because K= βS, equation (1). We did not include in the final models quadratic terms for the abiotic variables (i.e., aridity, sand) with the exception of soil pH for the sake of simplicity.

We used a model selection procedure based on AICc (ΔAICc < 2)63 to select the best predictors of multifunctionality. This procedure was performed using the function dredge in the R package MuMIn64. Model averaging was performed based on AICc weights when multiple models were selected. Model residuals were inspected for constant variance and normality. All predictors and response variables were standardized before analyses using the Z-score to interpret parameter estimates on a comparable scale. Predictors were log-transformed when necessary before analysis to meet the assumptions of the tests used.

We evaluated the relative importance of the predictors under consideration as drivers of multifunctionality. To do so, we calculated the relative effect of the parameter estimates for each of the predictors compared to the effect of all parameter estimates in the model. This method is similar to a variance decomposition analysis because we transformed all predictors to Z-scores prior to analysis. The following five identifiable variance fractions were examined: i) spatial variables, ii) abiotic variables, iii) species richness, iv) mean and variance, and v) skewness and kurtosis of trait-abundance distributions.

We evaluated how the observed SKRs explained variations in multifunctionality using the final best-fitting models. Whenever kurtosis and kurtosis2 for SLA and plant height were selected in the averaged model, we used their averaged parameter estimates to model its effects on multifunctionality. All other predictors selected in the averaged model were treated as constant and fixed to their mean (i.e. 0 since all predictors were transformed to Z-scores, see Supplementary Note 4 for more details). Finally, we used the same null model approach developed at the community level to test whether the relationship between trait-abundance distributions and multifunctionality could be observed by chance (see Supplementary Note 5 for more details).