Functional diversity drives ecosystem multifunctionality in a Pinus yunnanensis natural secondary forest

It is essential to understand how the loss of biodiversity impacts both ecosystem function (EF) and multifunctionality (EMF). Previous studies have mostly focused on predicting how species richness (SR) impacts EMF, while the effect of functional diversity (FD) on EMF remains unclear. Specifically, we know little about the primary functional drivers impacting EMF compared with SR. Therefore, we analysed 8 ecosystem functions within 58 natural secondary forest plots to investigate the effect of FD on both individual EF and EMF. Our results suggest that SR and FD had very significant positive effects on plant phosphorus, soil available phosphorus, and soil total nitrogen. FD explained significantly more variations in these functional responses than SR for individual ecosystem functioning. We also used a multiple threshold approach to test the effect of SR and FD on EMF. We found that FD and SR were positively related to EMF regardless of whether low-level function or high-level function was desired, but FD had a larger effect than SR. Based on the averaging approach, OLS regression, multivariate linear regression model and random forest analysis, we found that SR and FD were both drivers of EMF but that FD had a stronger effect and could explain more variation. As such, we conclude that FD drives ecosystem multifunctionality more than SR.

It is essential to understand how the loss of biodiversity impacts both ecosystem function (EF) and multifunctionality (EMF). Previous studies have mostly focused on predicting how species richness (SR) impacts EMF, while the effect of functional diversity (FD) on EMF remains unclear. Specifically, we know little about the primary functional drivers impacting EMF compared with SR. Therefore, we analysed 8 ecosystem functions within 58 natural secondary forest plots to investigate the effect of FD on both individual EF and EMF. Our results suggest that SR and FD had very significant positive effects on plant phosphorus, soil available phosphorus, and soil total nitrogen. FD explained significantly more variations in these functional responses than SR for individual ecosystem functioning. We also used a multiple threshold approach to test the effect of SR and FD on EMF. We found that FD and SR were positively related to EMF regardless of whether low-level function or high-level function was desired, but FD had a larger effect than SR. Based on the averaging approach, OLS regression, multivariate linear regression model and random forest analysis, we found that SR and FD were both drivers of EMF but that FD had a stronger effect and could explain more variation. As such, we conclude that FD drives ecosystem multifunctionality more than SR.
Biodiversity is crucial for sustaining ecosystem processes and functioning 1 . Biodiversity is declining at an unprecedented rate and will continue to decline over the 21 st century 2 . As such, understanding how the loss of biodiversity impacts ecosystem functioning has become an important goal in biology and conservation 3 . The study of the relationship between biodiversity and ecosystem functioning (hereafter, BEF) has been ongoing for over twenty years 4,5 . Many empirical studies have demonstrated that a loss of biodiversity can lead to reductions in ecosystem functioning 6 and the ecosystem services that provide a multitude of benefits to humans 7,8 . Progress in studies focused on BEF has been spurred with the deepening of understanding, and scientists found that ecosystems can perform multiple functions or services simultaneously 1,9 . Thus, multifunctionality (hereafter, EMF) has been generally accepted in recent years and is defined as "the simultaneous provision of multiple functions" 10 . Because different species have different functional traits, they could contribute different ecosystem functions. Thus, biodiversity might have an important effect on the multifunctionality 11 .
Most experiments conducted on the relationship between biodiversity and ecosystem multifunctionality (hereafter, BEMF) have focused on species richness (hereafter, SR) as the primary metric of biodiversity 1,[9][10][11][12][13][14] because SR is defined as the number of species and is easily measured. However, some studies suggest that functional diversity (hereafter, FD) is a stronger predictor of ecosystem functioning [15][16][17] . The effect is most immediate mainly because the effect of FD is the result of interactions between species and their environment relative to the species traits 15 . The reason for FD as a stronger predictor may be due to high FD causing an increase in the utilization efficiency of environmental resources, which can subsequently promote ecosystem productivity and strengthen defences against diseases, insect pests, and disturbances 18,19 .
Developing a better understanding of the links between diversity and ecosystem functioning is especially important in forest ecosystems. Forests provide innumerable benefits for humans, including timber, nutrient and water cycling, and recreational opportunities, among many others 8,13 . Despite this, previous studies on BEMF have mostly focused on experimental grassland ecosystems 1,20-23 . Regarding forests, more attention has been paid to temperate and boreal forests. Most of studies proved that higher biodiversity was necessary to maintain multiple ecosystem functions among different ecosystems, but little attention has been paid to subtropical forests 8,[24][25][26] . We noticed that some scientists, such as Mouillot et al. 27 , Valencia-Gómez et al. 28 and Finney & Kaye 29 , conducted many studies on the influence of FD on EMF. These studies found that FD could increase multifunctionality. However, this result was found in grasslands, dry lands or agricultural systems, and we did not identify related studies on forests to date.
Here, we seek to help fill this gap using data collected from a Pinus yunnanensis forest in southwest China 30 . We used data collected from 58 natural secondary forest plots located in Qiubei and Shuangbai County in Yunnan Province, China ( Fig. 1) and examined eight individual ecosystem functions (plant nitrogen, plant phosphorus, soil hydrolysable nitrogen, soil available phosphorus, soil total nitrogen, soil total phosphorus, soil total carbon, and woody plant biomass). Aboveground biomass was the primary key function in the BEMF study 8,31 . In the soil samples, soil total carbon, soil total nitrogen, soil hydrolysable nitrogen, soil total phosphorus and soil available phosphorus represent good substitutions related to carbon (C), nitrogen (N) and phosphorus (P) cycling 9,21 . Nitrogen and phosphorus are two important nutrient elements for plants. Nitrogen mainly determined photosynthetic C fixation and plant productivity 32 . Phosphorus is considered necessary for the storage and reproduction of genetic information and is mainly involved in energy-related processes of living organisms 33 . Regardless of whether nitrogen or phosphorus was limited, plant growth was affected. Plant nitrogen and plant phosphorus are considered nutrient pools in aboveground biomass 21 . We used this unique dataset to address two primary questions: (1) How do SR and FD influence individual ecosystem functions?; (2) Does FD predict EMF better than SR? We expected to identify the relative importance of biotic and abiotic variables on EMF in a Pinus yunnanensis natural secondary forest.
BEMF relationships based on the averaging approach. Both SR and FRic were related to EMF. SR explained 16.11% of the variation in EMF, while FRic explained 20.72% of the variation. SR was significantly positively correlated with FRic (R 2 = 0.7847, P < 0.001) (Fig. 2).

Effect of SR and FRic on EMF.
SR had a positive effect on EMF at thresholds between 12% (minimum threshold, T min ) and 91% (maximum threshold, T max ) (Fig. 3a), and FRic had a positive effect on EMF at thresholds between 8% (minimum threshold, T min ) and 97% (Fig. 3b). The relationship between SR and EMF peaked at www.nature.com/scientificreports www.nature.com/scientificreports/ a threshold of 55% (threshold of maximum diversity effect, T mde ), where each additional species was associated with 0.22 functions provided at levels above the threshold (realized maximum effect of diversity, R mde ) (Fig. 3a). The relationship between FRic and EMF also peaked at a threshold of 55%, but each increase in FRic was associated with 0.31 functions (Fig. 3b). The maximum possible effect of SR (percentage of maximum possible diversity effect, P mde ) for our experiment was 36.75%, and the P mde of FRic was 41.85% (Fig. 3).
Estimation of the relative contribution of the biotic and abiotic factors to EMF. The results of the multivariate linear regression model suggested FRic of the biotic factor had a significant positive effect on EMF.   www.nature.com/scientificreports www.nature.com/scientificreports/ The absolute value of the standardized regression coefficients of FRic (0.6535) was the biggest among these variables, indicating that FRic had the strongest impact on EMF compared with other factors. SR of the biotic factor and the abiotic factors, including soil pH, MAP and MAT, had no significant effect on EMF (Table 2).
According to random forest analysis, we found that FRic had the highest importance for EMF in all variables. In contrast, SR had the lowest importance among all selected variables. For the abiotic factors, the rank of importance was MAP, MAT, and soil pH (Fig. 4).  (Table 3).

Discussion
The primary goal of this work was to determine the relationship between biodiversity and multiple measures of ecosystem function. We found that both SR and FRic were positively correlated with three of the eight potential components of ecosystem functioning: plant phosphorous, soil available phosphorus and soil total nitrogen. Furthermore, these effects were observed when combining different measures and calculating ecosystem multifunctionality. Importantly, in all cases, we found that FRic was a better predictor than SR.    www.nature.com/scientificreports www.nature.com/scientificreports/ Functional diversity is thought to be declining worldwide 34 . The number of studies about the relationships between functional diversity and ecosystem functioning is increasing. However, the role of functional diversity in EMF has not been fully explored. This notion is particularly true for forest communities, which have been less well studied within the context of BEMF. In this study, we found that several individual ecosystem functions, such as plant phosphorus, soil available phosphorus and soil total nitrogen, were sensitive to changes in SR and FRic. For the two metrics of biodiversity, SR and FRic, the relationship between them and the individual ecosystem functions was relatively consistent. However, FRic could explain more variations in these three ecosystem functions than SR. These results were consistent with most previous studies 17,35,36 . It is worth mentioning that SR and FRic had no effect on woody plant biomass in this study, and the result was not consistent with previous conclusions of BEF studies 37 . This difference could be due to the structure of the Pinus yunnanensis natural secondary forest. In this forest type, Pinus yunnanensis held an absolutely dominant position relative to other species 30 . Pinus yunnanensis contributed 91.73% of the biomass in our study, so increasing or decreasing one species might not have a significant effect on the biomass. When considering multiple ecosystem functions simultaneously, biodiversity likely has a stronger effect 22 . In our study, this conclusion was confirmed regardless of whether SR or FRic served as the metric of biodiversity.

Factors
Our results also indicate no strong trade-offs between different functions in the Pinus yunnanensis natural secondary forest (Table S2). Normally, trade-offs between different functions are commonplace within local communities 1,20,38 . However, at larger spatial scales, heterogeneity in a community composition causes different parts of the landscape to provide different ecosystem functions under the precondition that different species provide different ecosystem functions 1,13,20 . In theory, this implies that all ecosystem functions could be maximized simultaneously if information on their causal relationships is known. These results were consistent with the results of van der Plas et al. 26 , who found that trade-offs were weak among different functions at larger scales.
We used the multiple threshold approach to evaluate whether and when SR or FRic were important drivers of EMF. We obtained a complete picture of how SR or FRic drove EMF by combining some metrics; FRic had an obvious positive effect on EMF at thresholds between 8% (minimum threshold) and 97% (maximum threshold). This range was wider than that for SR. Byrnes et al. 10 indicated that if the T min was low, the T max was high, and both the T mde and the R mde remained high, diversity would have a strong effect on EMF. Compared with that observed for SR, these parameters all had the strongest effect at a moderate threshold (T mde = 55%). However, R mde and P mde values of FRic were greater than those of SR; therefore, our results at least showed that the effect of FRic on EMF was greater than that of SR.
In the multivariate linear regression model and random forest analysis, we all found that FRic was the main driver in BEMF. Our result is consistent with previous studies [27][28][29] . Compared with taxonomic diversity for grassland, dry land, agricultural or subtropical forest ecosystems, functional diversity is the strongest predictor of EMF and has a positive effect. For EFs, in our study, functional diversity has a positive effect with the expectation of WPB (Table 3) because functional diversity is based on functional traits. Functional traits are the key mechanism by which single species and groups of species affect ecosystem functions 39 .
It should be noted that there is an inevitable connection between SR and FRic 40 . In this study, we observed a strong correlation between these variables (P < 0.001, R 2 = 0.7847). Because FRic is influenced by SR, it is somewhat difficult to tease apart the fraction of variation that is uniquely explained by functional diversity 41 . This problem could potentially be resolved using experimental studies in combination with variance partitioning analyses.

Conclusions
We found that 3 of 8 ecosystem functions, PP, SAP and STN, were strongly correlated with SR and FRic. However, FRic could explain more variations in these three ecosystem functions. Using a multiple threshold approach, we proved that the effect of FRic on EMF was greater than the effect of SR. In addition, we found no strong trade-offs among ecosystem functions in the Pinus yunnanensis natural secondary forest. Thus, all ecosystem functions may be maximized simultaneously. Moreover, our findings provide strong evidence that FRic drives EMF more than SR. For ecosystem managers, our findings are useful. Functional diversity was based on multiple functional traits. Thus, to increase EMF, the managers could start with species assemblages and chose species with functional traits as different as possible. To conserve biodiversity, we reinforce the view that if conserving or promoting biodiversity is the target in an ecosystem, ensuring a moderate level of EMF is at least a precondition 25 .

Materials and Methods
Plot selection and sampling. Field surveys were conducted in April 2015. In total, 58 forest plots measuring 20 × 20 m were established under the most representative vegetation over an extensive area (approximately 9, 000 km 2 ). Across all plots, the mean annual temperature ranged from 15.7 to 18.2 °C, the mean annual precipitation ranged from 954 to 1202 mm, and the altitude ranged from 1065 to 2125 m. For our study, the number of plots was sufficient to statistically detect the relationship between biodiversity and ecosystem functioning 42 . We selected plots with a gradient of plant species to detect the relationship between biodiversity and ecosystem functions and tried to minimize other confounding factors, for example, selected the forest of the same age, undisturbed soil and plant community structure. Within the plots, we collected information on species identity, height, DBH, and spatial location for all individuals >1 cm DBH. For each forest plot, we obtained the MAT and MAP data from ClimateAP (University of British Columbia, Vancouver, British Columbia, Canada) developed by Wang et al. 43 .
For all species, we also collected 50-100 sun-exposed mature leaves from the middle or upper part of the trees from five to ten individuals. All samples were carefully placed into paper bags and marked. We then measured leaf morphological characteristics using a laser area metre (LI-COR 3100C Area Meter, LI-COR, USA). The leaf area of Pinus yunnanensis was estimated as a cylinder, so we measured the diameter (d) and the length (L) by a Vernier www.nature.com/scientificreports www.nature.com/scientificreports/ calliper (precision: 0.02 mm). The leaf area of Pinus yunnanensis was calculated with the following formula: LA = 2π 3 /9·dL 44 .
All leaf samples were dried for 72 h at 60 °C, and the dry weight for each leaf was measured immediately with an electronic balance (precision: 0.0001 g). Afterwards, all samples were ground to a fine powder using a ball mill (NM200, Retsch, Haan, Germany) to measure plant N and P content. We also obtained 5-10 branches (1 cm ≤ DBH ≤ 2 cm) from some shrub species (DBH ≤ 5 cm). We used the water replacement method to measure branch volume after removing the peel. Then, all branches were dried for 72 h at 103 °C and subsequently weighed. The density of the branches was calculated as the ratio of the dry weight to the volume. For trees and shrubs >5 cm DBH, we sampled the wooden core by the growth cone. The volume was measured using the diameter and the length of the wooden core, and the wood density was calculated as the ratio of the dry weight to the volume of the wooden core. All functional traits were measured following standardized protocols 45 . For each plot, we collected 5 soil samples (0-20 cm) and mixed them evenly. The soil samples were sieved by a 2-mm mesh and air-dried for one month prior to physiochemical analyses.
Defining biodiversity. We used species richness (SR) and functional richness (FRic) as our measures of biodiversity. To correspond to species richness, FRic was used to measure functional diversity. We chose three commonly used plant functional traits to calculate functional diversity, including leaf area (LA, mm 2 ), specific leaf area (SLA = leaf area/dry weight, mm 2 ·mg −1 ), and wood density (WD = wood dry weight/wood volume, g·cm −3 ). These traits were expected to be related to the acquisition and utilization of resources in plants, including trees 46 .

Measurement of individual ecosystem functions.
We measured eight ecosystem properties or characteristics (used as proxies of functions) linked to total nutrient pools, nutrient cycling, and biological productivity. These ecosystem functions were as follows: plant nitrogen, plant phosphorus, soil hydrolysable nitrogen, soil available phosphorus, soil total nitrogen, soil total phosphorus, soil total carbon, and woody plant biomass. These ecosystem functions were consistent with previous studies 21,47 . Soil total carbon, plant nitrogen and soil total nitrogen were measured by combustion in a CHN elemental analyser (2400II CHN elemental analyser, PerkinElmer, Boston, MA, USA). Soil available phosphorus was measured by the Olsen method. Soil total phosphorus and plant phosphorus were measured by molybdenum antimony blue colourimetry. Soil hydrolysable nitrogen was measured by alkaline hydrolysis diffusion. Based on the DBH and height of every individual plant, we estimated the biomass of the 54 investigated species of woody plants by the growth equation (Table S1). If some species had no growth equation, we chose the growth equation of a similar species instead 48 . For less frequently occurring species, we chose the growth equation proposed by Ali et al. 49 to estimate the biomass following the suggestion of Lambert et al. 50 .
Prior to calculating the metrics of multifunctionality, we first analysed the correlations between the eight functional variables to determine whether they were strongly correlated with one another (Table S2). Among the 28 pairs of relationships, 12 pairs showed significant relationships. However, one exception was noted in that the relationship between soil total carbon and plant nitrogen was slightly negative (r = −0.282, p < 0.05). Therefore, this finding suggested the lack of no strong trade-offs among the different functions 51 .
Assessing EMF. We used two complementary approaches to evaluate the role of SR and FRic in driving EMF. The first was the averaging approach 9 . To obtain an averaged multifunctionality index for each site, we first calculated the standardized Z scores for the eight variables. The Z scores were then averaged to obtain a multifunctionality index 9 .
We next used the multiple threshold approach 10 . Here, we tallied the number of measured functions that simultaneously exceeded multiple critical thresholds, and the threshold was defined as a given percentage of the highest performance for each function. To eliminate potential outliers, the highest performing was defined as an average of the top five values for each function 52 . For our analysis, we selected threshold values ranging from 1 to 99%.
For these two approaches, the averaging approach has been widely used in analyses of multifunctionality 12,28,53 , and the index has good statistical properties. Specifically, it allows us to directly estimate the ability of a community to sustain multiple functions simultaneously 9,21 . It is more suitable for linear model analysis 10 . To date, in the reported approaches, the multiple threshold approach is the most comprehensive method to evaluate BEMF relationships 10,25 . This approach can provide more powerful information and flexibility than other approaches 10 .
Statistical analyses. We first evaluated the relationships between biodiversity (i.e., SR and FRic) and individual ecosystem functions using Pearson correlation and OLS regression analysis. The relationships among SR, FRic and EMF were evaluated using simple OLS regressions. We then used the multiple threshold method to evaluate the effect of SR and FRic on EMF. T min is the minimum threshold where SR or FRic begins to have an effect. T max is the threshold beyond which the slope first declined and was not significantly different from zero. T mde is the threshold where SR or FRic had its strongest effect. R mde is the realized maximum effect of SR and FRic, which is the strength of the linear relationships of SR or FRic-multifunctionality where SR or FRic has the strongest effects. P mde is the percentage of the maximum possible SR or FRic effect 10 . According to these parameters, the effect of SR on EMF could be compared with the effect of FRic on EMF 38 . We conducted a multivariate linear regression model to evaluate the comprehensive effects of biotic factors (FRic, SR) and abiotic factors (Soil pH, MAP and MAT) on EMF in the Pinus yunnanensis natural secondary forest. To select the best model, we employed different alternative models, such as abiotic factors entered first or species richness entered first into the multivariate linear regression model (Table S3). Then, we used random forest analysis to tease apart the relative importance of these variables on ecosystem multifunctionality. Random forest analysis was used widely because it could determine the effect of each predictor variable individually and could also interact with other predictors in a multivariable