Robust leaf trait relationships across species under global environmental changes

Recent studies show coordinated relationships between plant leaf traits and their capacity to predict ecosystem functions. However, how leaf traits will change within species and whether interspecific trait relationships will shift under future environmental changes both remain unclear. Here, we examine the bivariate correlations between leaf economic traits of 515 species in 210 experiments which mimic climate warming, drought, elevated CO2, and nitrogen deposition. We find divergent directions of changes in trait-pairs between species, and the directions mostly do not follow the interspecific trait relationships. However, the slopes in the logarithmic transformed interspecific trait relationships hold stable under environmental changes, while only their elevations vary. The elevation changes of trait relationship are mainly driven by asymmetrically interspecific responses contrary to the direction of the leaf economic spectrum. These findings suggest robust interspecific trait relationships under global changes, and call for linking within-species responses to interspecific coordination of plant traits.

L eaf traits represent plant functional strategies and have fundamental effects on vegetation properties and ecosystem functions [1][2][3] . Divergent leaf traits between species present strong trade-offs that called the leaf economic spectrum (LES) 2 . For example, plants with low specific leaf area (SLA) and leaf nitrogen content (N m ) have slow photosynthetic returns, while plants with the opposite traits have fast photosynthetic returns [4][5][6] . Large trait variations between species typically represent the evolutionary divergences resulting from genotypes or species turnover. The plant LES integrates trade-offs of leaf traits and thus provides a useful framework for elucidating leaf-to-ecosystem scaling and for modeling vegetation functional diversity and dynamics in a changing climate [7][8][9][10] . In addition, correlations among leaf traits can provide significant constraints on the estimates of vegetation-atmosphere carbon exchange [11][12][13] .
Leaf traits of a given species can also vary widely in response to environmental changes through a diverse array of physiological, behavioral, and ecological mechanisms 14,15 , which is defined as leaf trait plasticity 16 . Understanding leaf trait plasticity is a major challenge for predicting plant responses to global environmental changes 17,18 . Various empirical trait-environment relationships 9,19,20 and trait-trait correlations 21,22 have been presented to incorporate trait variations into mechanistic vegetation models. However, the extant patterns are the results of historical evolutionary selections from their specific environmental conditions 23 . Leaf trait plasticity may differ in magnitude and even direction from the existing trait-environment relationships 24,25 . Thus, it is still unclear whether the plasticity-caused trait changes would follow the principles derived from current LES analyses. In fact, the uncertain changes of trait-pairs in response to environmental changes have long been limiting the predictive application of the interspecific trait relationships 26 . Furthermore, the projected novel climate may lead to nonanalog plant functional types and/or trait combinations [27][28][29] . A better understanding of how leaf traits and their relationships in response to future environmental changes is crucial for predicting vegetation distribution and ecosystem function 30 .
With the advent of global change manipulative experiments which characterize the effects of environmental changes on leaf traits, directly quantifying leaf trait plasticity is possible 31,32 . Here, we combine leaf traits data from 102 field experiments and 108 environmentally controlled experiments (Supplementary Table 1 and Fig. 1a). The field experiments in this database cover extensive climatic regions, with mean annual temperature (MAT) ranging from −9 to 27°C, and precipitation from 40 to 4158 mm (Supplementary Fig. 1). The environmental temperature of environmentally controlled experiments ranges from 5 to 33°C and relative humidity from 27 to 97%. In total, we compile records of 515 species from 90 families, 60% of which from environmentally controlled experiments and 40% from field experiments (Inset in Fig. 1a). Based on these data, we quantify the plasticity of net photosynthetic rate (A m ), N m , and SLA to multiple global environmental changes. In addition, we analyze if the changes in traitpairs tend to be follow the directions of the interspecific trait relationships. Furthermore, we test whether the intrinsically interspecific trait relationships hold under future environmental changes. Leaf traits show allometric relationships with the form of power function: 1,2 y ¼ cx k . In the power law, c is the coefficient and k the dimensionless scaling exponent 33 . On the logarithmic axes of trait relationships log y ¼ k log x þ b; b log c ð Þ , slope (k) represents the fundamental mechanisms of the bivariate relationship and slope elevation (intercept, b) indicates the resource-use efficiency 34 . As illustrated in Fig. 1b, if both the slope and elevation of the log-log axes have no response to environmental changes, then the empirical interspecific trait relationship stays unchanged. Alternatively, if the slope keeps constant but the elevation shifts, then the scaling exponent remains robust but the resource-use efficiency varies. Conversely, if the slope changes significantly, then the empirical trait relationships are shifted. In this study we show that, although the direction of trait change in response to environmental changes varies enormously between species and generally does not follow the LES, the slopes of interspecific trait relationships hold stable.

Results
Response of leaf traits to global environmental changes. The response ratios of leaf traits to global environmental changes were normally distributed (Fig. 2a). On average, experimental warming had no significant effect on A m and N m , except for a slightly increase in SLA values by 6.9% (Fig. 2b). In contrast, all trait values varied significantly when the plants experience drought condition. Experimental drought significantly decreased A m (−38.3%) and SLA (−8.7%), respectively, but increased N m (6.5%). The elevated atmospheric CO 2 concentration roughly induced 16.1% and 12.6% reductions in N m and SLA but increased A m by 12.6%. Nitrogen addition significantly increased N m (34.0%) and A m (12.8%), respectively. However, the commonly measured morphological trait SLA did not shift significantly under nitrogen addition.
Plant growth environment (field or environmentally controlled) and treatment conditions (treatment strength and duration) were important factors affecting trait response 35 . Our results showed that field experiments presented larger effect on leaf photosynthesis, but lower impact on leaf nitrogen and SLA than environmentally controlled experiments ( Supplementary  Fig. 2). Furthermore, the responses of some leaf traits were affected by treatment strength of warming, drought, and nitrogen addition as well as the experimental duration of elevated CO 2 (Supplementary Table 2, Supplementary Fig. 3). The percentage changes in A m was negatively correlated with warming strength and temperature conditions, and the warming effects varied from positive to negative with the increase of the temperature ( Supplementary Fig. 3a). The effect of drought on A m increased linearly with enhanced drought strength (R 2 = 0.21, P < 0.01). Similarly, the positive responses of leaf nitrogen to nitrogen addition increased linearly with applied dose (R 2 = 0.14, P < 0.01). Our results also revealed that response of SLA to eCO 2 logarithmically decreased with increasing experimental duration.
The direction of changes in pairwise traits under global environmental changes. Importantly, we also tested the consistency between direction of trait plasticity and the interspecific trait relationships. We plotted the direction of trait plasticity induced by environmental changes for each species into trait-trait space, and found that the changes in trait-pairs were divergent between species. In addition, we found that only the plasticity of N m -SLA relationship under eCO 2 and the plasticity of A m -SLA under nitrogen addition tended to be follow the LES direction (Fig. 3i, j). Most of the trait plasticity was equally abundant in all directions of the trait-trait space under warming and drought ( Fig. 3a-f). Furthermore, the directions of trait plasticity under eCO 2 and nitrogen addition were mainly asymmetric distributed on one side of the LES (Fig. 3g-l).  Table 3). The elevation for A m -SLA correlation was slightly altered by drought ( Fig. 4d and Supplementary Table 3). The relatively lower elevation revealed a reduction in photosynthetic rate at a given unit of SLA under the drought treatment. Similarly, the increased CO 2 treatment altered the A m -N m and A m -SLA relationships via significant enhancement on their elevations (Fig. 4g, h and Supplementary Table 3). The increased elevations respectively indicated higher photosynthetic rate for a given leaf nitrogen concentration and higher photosynthetic rate for the equivalent leaf area. Meanwhile, the nitrogen addition significantly lowered the elevations of A m -N m correlation, but lifted the elevation of SLA-N m relationship (Fig. 4k, l and Supplementary Table 3). Overall, the elevation of A m -N m correlation was most sensitive to global environmental changes, with the elevation changes were −0.19 and 0.15 under nitrogen addition and elevated CO 2 , respectively, (Fig. 5). The elevation of N m -SLA relationship was relatively robust, only significantly altered by nitrogen addition. Note that the relationship change of SLA-N m under warming was not detected because of their nonsignificant correlation (Supplementary Table 3).
Recent studies have emphasized the equivalent importance of mass-and area-based interspecific trait relationships 36,37 . Here, we also tested the robustness of area-based interspecific trait relationships under environmental changes. First, we reconfirmed the disparate patterns of mass-and area-based trait relationships 3,37 ( Supplementary Fig. 4). The area-based A-SLA and A-N were weakly related, whereas stronger correlations for mass-based traits were presented. In addition, contrasting relationships were found between the area-and mass-based N-SLA. Then we showed that the slopes of area-based N a -SLA relationships remained unchanged (all P > 0.05, Supplementary Table 4). In addition, we found different elevation changes for mass-and area-based N-SLA relationship under elevated CO 2 ( Supplementary Fig. 4o, r). The robust N m -SLA correlation under eCO 2 resulted from the proportional changes of N m and SLA along the original axis, while larger eCO 2 -dependent decrease in SLA than N a was responsible for the significant reduction in elevation of N a -SLA.
It has been confirmed that species from different functional groups and growth environment (field and environmentally controlled) may have different eco-physiological constraints when experiencing global environmental changes 38,39 . We found that interspecific trait relationships have consistent slopes between angiosperm and gymnosperm woody species as well as between dicotyledons and monocotyledons were consistent. Only the slopes of the N m -SLA relationship were different between C 3 and C 4 herbs (Supplementary Table 5  woody but not gymnosperm woody species ( Supplementary  Fig. 6a). Dicotyledons and monocotyledons showed similar trait adjustment strategies under eCO 2 and nitrogen addition ( Supplementary Fig. 6b). We also detected a lower sensitivity of the photosynthetic nitrogen use efficiency in C 4 than C 3 plants to eCO 2 ( Supplementary Fig. 6c). Overall, different response among functional groups have implications for how these plant groups are likely to perform in future climate change. In addition, we found that the slopes of trait relationships from field and environmentally controlled experiments were both unchanged (all P > 0.05, Supplementary Table 9), and their elevations changed in the same direction ( Supplementary Fig. 7a). The impact of treatment strength on the elevation shift of trait relationship was only detected in A m -SLA under drought (Supplementary Table 10 and Supplementary Fig. 7b).

Discussion
This study shows the high plasticity of leaf economic traits under global environment changes. Plants generally have great plasticity in leaf characteristics to optimize their function under the prevailing changes in environmental conditions 40 . Globally, climate warming showed zero-sum impact on leaf traits (Fig. 2a). Our results corroborate the dependence of warming effects on temperature conditions 41 , and show that warming accelerate leaf photosynthesis in cold environments, but negatively affect leaf photosynthesis in warmer climates 35,42 (Supplementary Fig. 3a). The negative effect of warming may result from warming-induced water deficit or the excessive temperature that beyond its optimum points 43 . In addition, both SLA and N m of plants are sensitive to elevated CO 2 . The opposite responses of leaf photosynthesis and other leaf traits strongly suggest that leaf traits move in the opposite direction of the common recognized LES under the increased CO 2 concentration 44,45 , as expressed by higher leaf photosynthesis and higher carbon investment in leaves. The reduced N m under elevated CO 2 has been attributed to the dilution effect due to rapid growth 46 or redistribution of ecosystem N stocks 47 . Moreover, our results show that plants tend to alter N m rather than SLA under higher nitrogen availability. SLA has been found to be a reliable indicator of plant resource-use strategy 48 . However, our findings corroborate a recent study 49 that leaf nutrient, but not SLA, is an appropriate indicator of plant response to increased nutrients inputs. Considering the extremely high variability of leaf photosynthesis due to stomata closure under drought, A m is not a reliable indicator of functional response under drought stress. As suggested by some previous studies 13,30 , the quantification of the leaf traits plasticity could also provide observational constraints on modeled vegetation and ecosystem function.
The unchanged slopes support the robustness of the reported trait relationships in disparate environmental conditions. We further propose a conceptual framework to show why the interspecific LES maintains despite the divergent trait plasticity (Fig. 3m). There are three response patterns of pairwise traits under environmental changes: (1) the direction of trait plasticity follows the LES; (2) the direction of trait plasticity is contrary to the LES with asymmetric responses; and (3) the direction of trait plasticity is contrary to the LES with symmetric responses. As shown by Fig. 3,  environmental change. The observed elevation changes are mainly driven by the asymmetric responses which are contrary to the direction of LES (Fig. 3d, g, h, k, l). The across-species LES could shift when the trait plasticity is symmetric between upward and downward directions of the LES.
The interspecific trait relationships are mostly consistent across experimental types (field and environmentally controlled) and functional groups (Supplementary Note 1). However, the elevation of trait relationship varies in some cases. For example, the elevation changes are higher in the environmentally controlled  Table 9), which might result from the faster growth rates of plants than that in the field 39 . The elevations of trait relationship show larger variations in angiosperms than gymnosperms (Supplementary Fig. 6a). This result suggests a more conserved trait adjustment strategy in gymnosperms than angiosperms, which has also been found in the root traits 50 . Due to the higher nitrogen use efficiency in C 4 than C 3 herbs 51 , it is expected that the elevation of A m -N m relationship is greater in C 4 than C 3 herbs ( Supplementary  Fig. 6c). However, eCO 2 treatment significantly increased the elevation of A m -N m relationship for C 3 but not C 4 herbs. This finding suggests that the A m -N m relationship could be more stable in C 4 than C 3 herbs under the higher atmospheric CO 2 concentration. It should be noted that the sample size of our database is not large enough to compare all trait relationships between the control and treatment groups for all functional groups. However, the findings of some differences in the elevation of trait relationship between functional groups imply the need to explore trait variations across phylogenetically distant species. In fact, some recent global analyses have shown the importance of phylogenetic factors in explaining plant response to global changes 52,53 .
The elevation of A m -SLA correlation represents the photosynthetic efficiency per-unit leaf mass, which is affected by the investment in photosynthetic mass per-unit leaf area 54,55 . Similarly, the elevation of A m -N m correlation indicates photosynthetic N-use efficiency, shaped by nitrogen allocation to Rubisco or Rubisco-use efficiency 56 . Plants are expected to increase investment in structural mass and therefore confer higher physical strength and greater resistance to drought stress 57 , which could decrease photosynthetic efficiency per-unit leaf area (Fig. 4e). Without evidence of significant increase in photosynthetic mass per-unit area or Rubisco nitrogen fraction 47 ( Fig. 4h) could be attributed to the increased Rubisco-use efficiency. The reduced photosynthetic N-use efficiency under nitrogen addition has been reported due to the decreased fraction of nitrogen to Rubisco 59,60 (Fig. 4j). In addition, the elevation changes in the N m -SLA correlation result from the variation of leaf nitrogen per-unit area, which is only significantly enhanced by nitrogen addition (Fig. 4l). Further researches are still needed for a deeper understanding of the physiological mechanism underlying the elevation change under future environmental changes.
The robust slopes of leaf trait relationship across species implie that the traits coordination could be used to predict ecological consequence of global environmental changes. In fact, plant traits are important parameters in regulating vegetation processes in the framework of Earth system models 11,61 . However, an increasing body of evidence has indicated the insufficient realism of current models for their ignorance of leaf traits variability and plasticity 28,62 . The widely demonstrated leaf traits plasticity can alter leaf structure and function, thus, plant productivity and land-atmosphere fluxes. Recent advances have incorporated leaf trait plasticity into models for a more realistic presentation of the responses of terrestrial ecosystem to climate change 25,30 . Currently, connecting trait relationships across plants, microorganisms and animals remains a big challenge for ecology and biogeochemical modeling 63 . Our study suggests that scaling trait relationships from intra-to interspecies and even up to the community or ecosystem level 64 is an important next step to implement trait-based approaches into modeling future dynamics of the Earth system.
In summary, the various patterns of the leaf trait relationships along climate gradient have long been limiting the predictive utility of such empirical correlations 15,34 . Such a weak predictive ability could largely stem from the divergent directions of the shifts in leaf traits under environmental changes among species. It is interesting that the divergent changes in pairwise traits have no impact on the slope of interspecific trait relationships at the global scale. Overall, this study indicates that the direction of changes in pairwise traits in response to global change varies enormously between species, and the within-species changes generally do not follow the direction of interspecific trait relationships. However, the slopes of interspecific trait relationships across species are stable under environmental changes, which indicate the fundamental mechanism of the bivariate relationship does not change. These findings underscore the importance of identifying the key ecological processes which link the changes of different traits within and between species.

Methods
Data compilation. We searched for peer-reviewed journal articles using ISI Web of Science and Google Scholar in March 2018 with no restriction on publication year. We focused on manipulative studies which included key global environmental changes: (1) warming, (2) eCO 2 , (3) drought, and (4) nitrogen addition. Drought here refers to rainfall reduction (field experiment) or irrigation reduction (environmentally controlled experiment). More than 4000 published articles reported changes of leaf traits (SLA, leaf nitrogen content (N), and/or net photosynthetic rate (A)) under experimental manipulations. To minimize publication bias caused by different data sources, we considered only articles with at least two of the three leaf traits in this study. After systematic screening, a total of 404 published articles were included in our database (Supplementary Table 1 and Supplementary Fig. 8).
Overall, the experimental approaches included greenhouse, growth chamber, pot, garden, and field habitat. In this study, the experiments conducted in garden and field habitat were defined as field experiments. The experiments conducted in greenhouse, growth chamber, and pot were classified as environmentally controlled experiments, in which the disturbances of the other variables were minimized. As a result, this study compiled a trait plasticity database from 102 field experiments and 108 environmentally controlled experiments ( Fig. 1a and Source Data). The consistency of trait relationships between field and environmentally controlled experiments has been tested before using the whole dataset in global analyses (Supplementary Note 1). The experimental condition of field experiment was characterized by MAT and mean annual precipitation, while the environmentally controlled experiment was described by experimental temperature (T) and relative humidity (%).
All original data were extracted from the text, tables, figures, and appendices of the publications. When data were graphically presented, GetData Graph Digitizer v2.26 was used to obtain numeric data (http://getdata-graph-digitizer.com/). Areaand mass-based measures of any traits were converted through the following relationship, where T m and T a are the values of trait T expressed on per-unit mass and per-unit area, respectively; SLA is the specific leaf area (in units of cm 2 g −1 ). Error propagation through the conversion is calculated using the following equations: Standard deviation (SD) of T m equal to, where SD T m , SD T a and SD SLA are the SD of trait T m , T a , and SLA, respectively.
Data modifications. Considering the statistical assumption of independence among observations, we used only one measurement of each species from the same study. For all the variables, if more than one observation were reported during the same experiment for the same species, a weighted average value was calculated by, with standard deviation where j is the number of observations, T i , SD i , and n i are the mean, SD, and sample size of the ith sampling data, respectively 65 .
For the potential nonindependence of the same species from different studies, we created another version of the dataset purged of intraspecific variation by replacing multiple observations per species with the mean of those values 3 . However, the two datasets showed same response patterns of trait relationships to different environmental changes (Supplementary Fig. 9 and Supplementary  Table 11). Therefore, considering huge experimental difference (field or environmentally controlled, different treatment strength, different genotypes, etc.) of the same species from different studies, we presented the results with the larger sampling size (i.e., keeping the same species from different studies) in the main text. We refer to any observation as a species observation in the subsequent analysis.
Statistical analysis. The effects of the global environmental changes on plant functional traits were quantified following the methods described by Hedges et al. 66 using Metawin 2.0. The natural logarithm-transformed response ratio (RR) was used to evaluate the environmental changes effects on leaf traits: where T t and T c are the experimental treatment mean and control mean, respectively. The variance (ν) is estimated by: where n t and n c represent the sample size, and SD t and SD c are the SD for the treatment and control variables, respectively. The weight of each RR is the reciprocal of its variance w ¼ 1 V À Á . Then the weighted response ratio (RR ++ ) is calculated as below (m is the number of groups and n is the number of comparison): with the standard error (SE) is calculated as Then the 95% confidence interval (95% CI) is RR þþ ± 1:96 S RR þþ À Á . The percentage changes of a variable were calculated as: The effects of environmental changes were evaluated as significant, if the 95% CI does not overlap zero. We also used Q-statistic 67 to test the heterogeneity of the effect sizes between different functional groups (angiosperm woody vs. gymnosperm woody, dicotyledons vs. monocotyledons, and C 3 herb vs. C 4 herb) and growth environment (field vs environmentally controlled conditions, low strength vs high strength, and treatment duration). If Q b is larger than a critical value, there would be significant difference between the categories. Statistical significance was tested at the P < 0.05 level.
SMA regression. Considering the concurrent errors in both axes, we used SMA regressions to quantify allometric relationships of pairwise traits under control and treatment conditions. The DOS-based SMATR package used for SMA regressions allows testing both for homogeneity among SMA slopes via a permutation test and for differences in SMA elevation via the SMA analog of standard ANCOVA 68 . When homogeneity was demonstrated (P > 0.05), a common slope was estimated. Elevation homogeneity comparisons were performed only when slopes were homogeneous. Where noted in the results, log 10 transformations were carried out on the original data and SMA regression then fitted between leaf functional traits.
Covariance error ellipse. Covariance error ellipses were created to summarize the distribution characteristics of leaf traits along the common axis: central tendency, dispersion, and directional trends. These measures define the axes of an ellipse encompassing the distribution of features. The covariance error ellipse was drawn with Matlab 69 .
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.