Tree Species Diversity and Stand Attributes Differently Inﬂuence the Ecosystem Functions of Pinus yunnanensis Secondary Forests under the Climate Context

: It has been widely reported that biodiversity, ecosystems, and functional traits are positively interrelated in natural forest ecosystems. However, it remains unclear whether these relationships should be expected in secondary forests. In this study, we hypothesized that the multifunctionality (EMF) is affected by the climate dependency of tree-species diversity and stand attribute diversity in a secondary forest dominated by Pinus yunnanensis . By using forest inventory data from a wide range of areas, we quantiﬁed the aboveground biomass, soil organic carbon, ratio of soil carbon and nitrogen, total soil nitrogen, total soil phosphorus, total soil potassium, tree-species diversity, and stand attribute diversity (i.e., individual tree-size variations). We also quantiﬁed the climate data, including the mean annual temperature (MAT), and mean annual precipitation (MAP). We found that a higher MAT directly constrains all the ecosystem multifunctionalities (EMFs) and three of the ﬁve single functions. A higher MAP was negatively correlated with all the EMFs and four of the ﬁve single functions, but indirectly through diversity indices. Stand attribute diversity better explained the EMFs rather than tree species diversity. Meanwhile, most of the single functions were highly correlated with stand attribute diversity rather than tree species diversity. These results highlight the importance of diversity in promoting forest multifunctionality and underscore the importance of the climate context in deﬁning EMF and shaping the relationship between diversity and ecosystem functions. We argue that the climate context should be taken into account when maximizing forest complexity, so as to enhance the multifunctionality of Pinus yunnanensis secondary forests.


Introduction
Forests harbor much of the world's biodiversity and are providers of essential functions and ecological services, including carbon stocks, nutrient cycling, water maintenance, non-timber forest production, and so on [1]. However, human activities including those associated with socioeconomic development, such as slash-and-burn agriculture, habitat destruction, and resource over-exploitation, have caused a drastic loss of terrestrial biodiversity and ecosystem changes [2,3]. The need to understand the interrelationship between biodiversity and ecosystem functions for ecology and biodiversity conservation has become urgent as we seek to forecast and prepare for the consequences of global biodiversity change [4]. The past few decades have witnessed a growing number of studies on the biodiversity-ecosystem-function relationship, ranging from experimental studies mostly in temperate grasslands [5,6] and forests [7] to observational studies in natural ecosystems [8]. A positive relationship between biodiversity and ecosystem function has among diversity, the climate context, and forest multifunctionality. Specifically, we investigated tree species diversity and stand structural diversity, overall, aboveground, and belowground multifunctionality, and several important individual functions. We employed six widely used forest functions of overall forest multifunctionality, consisting of aboveand belowground functions and processes [26,[38][39][40].
We hypothesized that the multifunctionality, including above-and belowground multiple and single functions of Pinus yunnanensis secondary forest can be promoted by both tree species diversity and stand structural attributes, while climate factors, namely, the mean annual temperate (MAT) and mean annual precipitation (MAP), have both direct and indirect effects on the functions mentioned above. To address the proposed hypothesis, we asked, based on the multiple paths in the conceptual model (see Figure 1 for an explanation), the following two research questions: (1) how are tree-species diversity, stand structural attributes, and forest functioning affected by MAT and MAP? And (2) how do tree species diversity and stand structure diversity interact and/or mediate the effect of the climate context on overall forest multifunctionality, as well as aboveground and belowground components?

Study Forests
Our study sites were located in four counties of Yunnan province in southwest China. The mean annual temperature of the study areas varies between 10.9 and 18.2 • C, and the mean annual rainfall ranges from 701 to 1203 mm, at an elevation of 1065 and 3003 m above sea level. Field surveys were conducted in 2020. In total, 122 Pinus yunnanensis-dominated forest plots (20 m × 20 m each) were established, representing the main forest types in all four counties (Table 1). Within the plots, all woody stems with diameters at breast height (DBHs) above 1 cm were identified, and the tree size, including the DBH and height, was also measured. The field surveys covered four different types of Pinus yunnanensis-dominated forest zones, i.e., Southeast Yunnan karst mountain valley forest zone (Qiubei, QB), Southwest Yunnan middle mountain forest zone (Shidian, SD), Basin valley forest zone of central Yunnan Plateau (Shuangbai, SB), and P. yunnanensis primary forest zone in Northwest Yunnan (Yunlong, YL).

Quantification of Forest Single Functions and Multifunctionality
Seven individual forest ecosystem functions were measured, including the aboveground biomass (AGB), the soil organic carbon (SOC), the total soil nitrogen (TSN), the total soil phosphorus (TSP), the total soil potassium (TSK), and the ratio of soil carbon and nitrogen (Table S2). These variables (or functions) were selected because of their direct or indirect relations to forest functioning, biological productivity, and soil nutrient cycling or soil fertility. Since the support capability of ecosystem services highly depends on these single functions, they could act as a proxy for forest functioning on the whole [39][40][41][42].
Using the global allometric equations, we calculated the total aboveground biomass (AGB) of each individual tree (DBH ≥ 1 cm) from their corresponding DBH, height, and woody density [43]. By using a soil auger 5 cm in diameter, we collected the soil samples at a depth of 0-20 cm from the four quarters and middle of each plot. Five soil samples were collected from each plot. After they were uniformly mixed, the soil samples were brought back to the laboratory and naturally dried, passed through a 0.25 mm sieve, and put into bottles for use. The contents of available carbon, total nitrogen (TN), total phosphorus (TP), hydrolytic nitrogen, and available phosphorus and the pH value were analyzed using the Kjeldahl method (HJ717-2014), sodium hydroxide alkali fusion, the molybdenum antimony colorimetric method, the alkali distillation potentiometric titration method, the liquor extraction method, and a pH meter (Thermo orion-868), respectively. The measurements of each sample were repeated three times.
The meteorological data for the study forests were calculated using ClimateChina v4. 40, a climate simulation software [44]. Combining the altitude and latitude of each plot, we measured the meteorological data for the annual mean temperature, annual precipitation, and effective accumulated temperature higher than 5 • C, etc.
We used multiple approaches including an overall averaging and multi-threshold approach to quantify the forest multifunctionality based on the six abovementioned functions [13], the aboveground and belowground forest multifunctionality, and the individual forest functions. To obtain the averaged forest multifunctionality index [40], we standardized all the variables, Z-scored, and then averaged before using the averaging approach. A similar approach was adopted to calculate both aboveground and belowground forest multifunctionality indices, with one index for aboveground and five for belowground functions. We used the multi-threshold approach (i.e., 25%, 50%, and 75% of the maximum observed functioning) developed by Byrnes et al. [13], taking the tradeoffs and synergies between different functions into account and to overcome some limitations in the averaged approach.

Quantification of Forest Diversity Attributes
We quantified four forest diversity attributes: the tree species diversity (Diversity), individual tree DBH variation (CV D ), individual tree height variation (CV H ), and individual stem abundance (Stem). By counting the number of species within each plot, we determined the species richness. To quantify the stand structural variations [23], we used the coefficient of variation (i.e., CV D and CV H ) for individual trees' DBHs and heights within each subplot.
We summed up the total number of stems in each plot as the stem abundance (Stem). CV D and CV H are highly correlated with Stem ( Figure S1), so we used only CV D and CV H in the final analysis as a latent variable (Size) responsible for quantifying both the horizontal and vertical niche space.

Statistical Analyses
To assess how predictor variables and response variables (single functions) were correlated with each other ( Figure S1), we calculated pairwise Pearson correlations. By using simple linear regression, we showed the bivariate relationships on the basis of the hypothesized paths in the conceptual model ( Figure 1), so as to complement the results from structural equation modeling (SEM). All the variables (MAT, MAP, CVdbh, CVheight, diversity, OMF, AMF, BMF, EMF 25%, EMF 50%, EMF 75%) were natural-log transformed and then standardized ( Table 2), in order to avoid collinearity issues and to obtain comparable standardized direct, indirect, and total effects (i.e., varying between 0 and 1 in either the positive or negative direction). Employing the multifunc [13] and Latent Variable Analysis (lavaan) [45] packages in R, respectively, we performed the analysis of the forest multifunctionality and single ecosystem multifunctionality. Using chi-square tests and associated p-values (i.e., p > 0.05 indicated acceptance of the SEM), we evaluated the SEM model fit [46]. R version 3.6.1 was used in all the analyses.

Results
The tested SEMs showed that MAT had a consistent negative direct effect on the overall (OMF), aboveground (AMF), and belowground multifunctionality (BMF), while MAP promoted AMF but not OMF or BMF (Figures 2, 3 and S14). We did not find a positive relationship between tree species diversity and ecosystem functions in all three facets. However, there was a strong positive relationship between tree diversity (Diversity) and stand attribute diversity in terms of the individual size variation (CV). Stand attribute diversity promoted all the considered ecosystem multifunctionality (EMF), and the effect was the strongest for AMF, followed by OMF and BMF (Figures 2, 3 and S14). Contrary to our expectations, a higher MAP constrained both tree species diversity and stand attribute diversity, and the effects of MAP on EMF were mediated by both tree species diversity and stand attribute diversity, whereas there was no such significant indirect relationship between MAT and EMF (Figures 2, 3 and S14). This indirect effect was greater in the relationship between MAP and AMF ( Figure 3) than in the relationship between either MAP and OMF or MAP and BMF (Figures 2 and S14).
Similar to the relationship between MAT and EMF, MAT showed a consistent strong negative effect on all the considered single functions (Figures S15-S19), while MAP had a significant positive relationship with the function of total soil phosphorus ( Figure S19). Contrary to the effect of tree species diversity on EMF, diversity increased both the function of soil carbon content ( Figure S15) and ratio between soil carbon ( Figure S16) and soil total nitrogen ( Figure S18) but constrained the function of soil total potassium ( Figure S17). Among all the single functions, stand attribute diversity only had a positive relationship with soil total nitrogen ( Figure S18). Meanwhile, the effects of MAP on single functions were mediated mainly by tree species diversity rather than stand attribute diversity (Figures S15-S19).

Discussion
Our study showed that secondary forest biodiversity promotes forest functioning for overall, aboveground, and belowground forest multifunctionality, and single functions (i.e., forest aboveground biomass and most of the soil nutrient functions), possibly due to the resource utilization complementarity (niche complementarity hypothesis) and decreased competition effects [47]. These results are consistent with previous findings in natural forest ecosystems, showing that species diversity positively contributes to multifunctionality in forests [1,48,49]. However, we found that a higher tree species diversity only directly increased the single functions of secondary forests but not multifunctionality, while stand attributes indirectly enhanced multifunctionality through a higher species diversity. The forest multifunctionality was promoted by tree species diversity better than stand attribute diversity, reflecting the idea that soil nutrients might be contributed to by different plant species via litterfall inputs through top-down processes in more diverse stands [50]. Forest multifunctionally could be enhanced by a higher species diversity in more diverse stands, which can coexist due to the increased available niches and confirm the complementarity effects [22]. Furthermore, the contribution of CV to individual aboveground functions and aboveground forest multifunctionality was higher than that to belowground functions, probably because light capture and use can directly influence the aboveground functions through stand structural complexity, reaching maximum crown occupancy [23].
Our secondary forest results underscore the fact that the abiotic environment (i.e., MAP and MAT) influence EMF both directly and indirectly through shaping the tree species diversity and stand attribute diversity at the community level. We found that both a higher MAT and MAP might constrain most of the forest functioning, both above-and belowground functionality, and most single functions. The finding of the constraint placed on EMFs by higher MAT was consistent with a study conducted in global dryland ecosystems [40], implying that the supporting ability related to aboveground biomass accumulation and soil nutrient cycling of secondary ecosystems will be decreased by an increasing temperature. A higher temperature is supposed to accelerate the activity and interaction among soil microbes and other invertebrate decomposers in soil, thus influencing soil carbon, nitrogen, and phosphorus cycling [51,52]. However, a higher temperature was co-varied with low rainfall in our studied ecosystem [34], which may cause higher aridity in soil and reduce the abundance and activity of microbes [53]. Previous empirical studies suggested that MAP might affect EMF positively [54] or neutrally [55], which is inconsistent with the overall negative relationship between MAP and EMF found in our secondary forest. Further studies are needed to examine the interaction of MAP and MAT, especially in dry and wet seasons, to better understand how the within-year climate may affect EMF in secondary forests.
The mediating effect of MAP on EMF in our study further confirmed that abiotic factors might indirectly affect ecosystem functions by altering the composition of communities [56]. There have been relatively fewer studies on the ecological mechanisms linking multiple climate factors and aboveground biomass or productivity [19,57]. However, there has been less discussion on the forest ecosystems across large-scale ecological gradients in relation to the direct effects of climate on aboveground biomass versus the indirect effects which are mediated by species diversity and stand structural complexity [58]. Here, we found that MAP constrained EMF mostly through the indirect effects mediated by stand attribute diversity, and negatively affected most single functions by the indirect effects mediated by tree species diversity. Our indirect climate hypothesis that variations in EMF across large-scale ecological gradients are best explained by an indirect influence of climate [59] is well supported by the consistent strong indirect effect of MAP on EMF via tree species diversity or stand structural complexity.
This study may have important research implications for biodiversity conservation, carbon stocks, and forest management in secondary forests. For instance, our results demonstrate that the climatic context affects ecosystem functioning not only directly, but also indirectly through diversity attributes, implying that future climate change may bring much influence to the widespread Pinus yunnanensis secondary forests. More specifically, the negative effect of temperature on ecosystem functioning suggests that global warming might constrain long-term carbon storage and soil nutrient cycling in this ecosystem; however, future increases in precipitation might affect the secondary forest in an opposite manner. In addition, both tree species diversity and stand attribute diversity are highly correlated with ecosystem functioning, implying that the future management of this secondary forest via increasing diversity is plausible. The ecosystem multifunctionality can be enhanced by maintaining high levels of stand structural complexity, at both the tree DBH and height levels, and the single functions may be improved by maintaining a higher species diversity. Overall, all the management regimes should be realized under a favorable climate context.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/su14148332/s1, Figure S1: Correlation relationships between all variables, Table S1: Direct, indirect, and total standardized effects on forest multifunctionality, multithreshold level, and each single function based on structural equation models (SEMs). Significant effects are at p < 0.05 (*), < 0.01 (**), and < 0.001 (***), Table S2: Model fit statistics summary of the tested structural equation models (SEMs) for overall (OMFA), above-ground (AMFA), and belowground (BMFA) forest multifunctionality, multi-threshold level (i.e., MFT25, MFT50, and MFT75) and each single function. Abbreviations: df, degrees of freedom; CFI, comparative fit index; SRMR, standardized root mean square residual; AIC, Akaike information criterion; R2 indicates the total variation in multifunctionality and also each single function that is explained by the combined independent variables, Figure S2: Variability of Diversity in four sites. Significant effects are at p < 0.001 (***) and < 0.0001 (****), Figure S3: Variability of CV DBH in four sites. Significant effects are at p < 0.001 (***) and < 0.0001 (****), Figure S4: Variability of Size height in four sites. Significant effects are at p < 0.0001 (****), Figure S5: Variability of Soil C:N in four sites. Significant effects are at p < 0.05 (*), < 0.01 (**), and < 0.0001 (****), Figure S6: Variability of Soil organic carbon in four sites. Significant effects are at p < 0.0001 (****), Figure S7: Variability of Total soil nitrogen in four sites. Significant effects are at p < 0.05 (*), < 0.001 (***), and < 0.0001 (****), Figure S8: Variability of AMF in four sites. Significant effects are at p < 0.05 (*), < 0.001 (***), and < 0.0001 (****), Figure S9: Variability of BMF in four sites. Significant effects are at p < 0.05 (*) and < 0.0001 (****), Figure S10: Variability of OMF in four sites. Significant effects are at p < 0.05 (*) and < 0.0001 (****), Figure S11: Location of research sites and relevant values of AMF, Figure S12: Location of research sites and relevant values of BMF, Figure S13: Location of research sites and relevant values of OMF, Figure S14: The structural equation models (SEMs) show that the belowground forest multifunctionality (BMF AVE ) interacts with tree species diversity (Diversity), individual tree size variation (Size), mean annual temperature (MAT), and mean annual precipitation (MAP). The significant paths (p < 0.05) are represented by solid arrows, and the nonsignificant paths (p > 0.05) are indicated by dashed arrows. The positive relationships are indicated by black paths, while the negative relationships are indicated by red paths. The direct and indirect effects are shown by bar charts, and relative contributions are shown by pie charts. The names of all variables match the colors provided in SEMs. Tables S1 and S2 present the model fit statistics for each SEM, Figure S15: The structural equation models (SEMs) show that the soil organic carbon is interrelated with tree species diversity (Diversity), individual tree size variation (Size), mean annual temperature (MAT), and mean annual precipitation (MAP). The positive relationships are indicated by black paths, while the negative relationships are indicated by red paths. The direct and indirect effects are shown by bar charts, and relative contributions are shown by pie charts. The names of all variables match the colors provided in SEMs. Tables S1 and S2 present the model fit statistics for each SEM, Figure S16: The structural equation models (SEMs) show that the ratio of soil carbon and nitrogen is linked with tree species diversity (Diversity), individual tree size variation (Size), mean annual temperature (MAT), and mean annual precipitation (MAP). The positive relationships are indicated by black paths, while the negative relationships are indicated by red paths. The direct and indirect effects are shown by bar charts, and relative contributions are shown by pie charts. The names of all variables match the colors provided in SEMs. Tables S1 and S2 present the model fit statistics for each SEM, Figure S17: The structural equation models (SEMs) show that the total potassium of soil is interrelated with tree species diversity (Diversity), individual tree size variation (Size), mean annual temperature (MAT), and mean annual precipitation (MAP). The positive relationships are indicated by black paths, while the negative relationships are indicated by red paths. The direct and indirect effects are shown by bar charts, and relative contributions are shown by pie charts. The names of all variables match the colors provided in SEMs. Tables S1 and S2 present the model fit statistics for each SEM, Figure S18: The structural equation models (SEMs) show that the total nitrogen of soil is linked with tree species diversity (Diversity), individual tree size variation (Size), mean annual temperature (MAT), and mean annual precipitation (MAP). The positive relationships are indicated by black paths, while the negative relationships are indicated by red paths. The direct and indirect effects are shown by bar charts, and relative contributions are shown by pie charts. The names of all variables match the colors provided in SEMs. Tables S1 and S2 present the model fit statistics for each SEM, Figure S19: The structural equation models (SEMs) show that the total phosphorus of soil interrelates with tree species diversity (Diversity), individual tree size variation (Size), mean annual temperature (MAT), and mean annual precipitation (MAP). The positive relationships are indicated by black paths, while the negative relationships are indicated by red paths. The direct and indirect effects are shown by bar charts, and relative contributions are shown by pie charts. The names of all variables match the colors provided in SEMs. Tables S1 and S2 present the model fit statistics for each SEM.