Leaf carbon and nitrogen stoichiometric variation along environmental gradients

Leaf stoichiometric traits are central to ecosystem function and biogeochemical cycling, yet no accepted theory predicts their variation along environmental gradients. Using data in the China Plant Trait Database version 2, we aimed to characterize variation in leaf carbon and nitrogen per unit mass (C_mass, N_mass) and their ratio and to test an eco-evolutionary optimality model for N_mass. Community-mean trait values were related to climate variables by multiple linear regression. Climatic optima and tolerances of major genera were estimated; Pagel's λ was used to quantify phylogenetic controls, and Bayesian phylogenetic linear mixed models to assess the contributions of climate, species identity, and phylogeny. Optimality-based predictions of community-mean N_mass were compared to observed values. All traits showed strong phylogenetic signals. Climate explained only 18 % of C:N ratio variation among species but 45 % among communities, highlighting the role of taxonomic replacement in mediating community-level responses. Geographic distributions of deciduous taxa were separated primarily by moisture and evergreens by temperature. C_mass increased with irradiance but decreased with moisture and temperature. N_mass declined with all three variables. C:N ratio variations were dominated by N_mass. The coefficients relating N_mass to the ratio of maximum carboxylation capacity at 25 ^∘C (V_cmax25) and leaf mass per area (M_a) were influenced by leaf area index. The optimality model captured 68 % and 53 % of variation between communities for V_cmax25 and M_a, respectively, and 21 % for N_mass. We conclude that stoichiometric variations along climate gradients are achieved largely by environmental selection among species and clades with different intraspecific trait values. Variations in leaf C:N ratio are mainly determined by N_mass, and optimality-based modelling shows useful predictive ability for community-mean N_mass. These findings should help to improve the representation of C:N coupling in ecosystem models.

1 Introduction

Nitrogen (N) has long been recognized as a key nutrient that influences photosynthesis, plant biomass, and carbon (C) allocation and therefore the terrestrial C cycle (Fernández-Martínez et al., 2014; Terrer et al., 2019). Many land surface models (LSMs) have recently incorporated representations of coupled C and N cycling, the intention being to increase the realism of model predictions of C cycling under climate change (Wiltshire et al., 2021). The leaf C:N ratio plays an essential role in this coupling; however it is often assigned a constant value per plant functional type (PFT), due to the lack of data and/or theory that would predict more realistic, continuous stoichiometric variation along environmental gradients (Meyerholt et al., 2020). One aspect of model uncertainty could be reduced if such variation were better understood and quantified (Boonman et al., 2020; Niu et al., 2023).

Many studies on leaf C and N along climate gradients have been carried out, but there is still no consensus on the major controls of such leaf traits at individual and community levels, hindering our understanding of trait–environment relationships (Anderegg, 2023). There is evidence for leaf stoichiometry being affected by many factors, including species identity, phylogeny, climate, and soil properties (Elser et al., 2010; Ma et al., 2018; Tang et al., 2018; Yang et al., 2016). The demonstrated roles of species identity and phylogeny indicate that leaf C and N contents (C_mass, N_mass) and their ratios are phylogenetically conservative (Sardans and Penuelas, 2014; Vallicrosa et al., 2021; Zhang et al., 2012). Differences among life forms and vegetation types have also been widely noted (Ma et al., 2018; Tang et al., 2018). On the other hand, few studies have examined the differences between evergreen and deciduous leaves, which are expected to diverge as they represent alternative life-history strategies, expressed in different responses of leaf mass per area (M_a) to climate (Kikuzawa et al., 2013; Wang et al., 2023). The patterns of stoichiometric response to environment remain inconsistent across studies. A potential contributory problem is their reliance on annual average climate variables, such as mean annual temperature – which does not accurately reflect actual growing-season conditions, especially in regions with cold winters (Körner, 2021) – and mean annual precipitation, which is generally not a good metric for plant-available moisture because it does not take account of the large variations in potential evapotranspiration (driven by solar radiation and temperature) across the world.

Many land surface models (LSMs) treat leaf C:N ratios as fixed parameters for PFTs (Boonman et al., 2020; Zaehle et al., 2014); some allow the C:N ratio to vary, within a prescribed range, based on C and N allocation to different tissues (Ghimire et al., 2016; Meyerholt and Zaehle, 2015; Smith et al., 2014; Wang et al., 2010). But the fixed-PFT schemes fail to capture the observed range of leaf stoichiometry within each PFT, while the dynamic schemes have not been extensively tested against observations. Responses of N use efficiency and net primary production (NPP) to elevated CO₂ vary considerably among models and are not always realistic (Zaehle et al., 2014). Here eco-evolutionary optimality (EEO) principles may help by providing a route towards testable, general trait predictions (Caldararu et al., 2020; Dong et al., 2022; Harrison et al., 2021; Xu et al., 2021). Caldararu et al. (2020) applied an optimality-based approach (maximizing carbon export and growth) to improve leaf N_mass prediction, but this analysis did not consider the large and potentially confounding effect of M_a variation with environment (Wang et al., 2023). We infer that there is still a need to investigate the eco-evolutionary basis of leaf C:N ratio variations and to reconsider how they are treated in LSMs (Sistla and Schimel, 2012).

In this study, we applied EEO principles to predict N_mass from a trait-correlation perspective. We assumed that the metabolic and structural components of leaf N are proportional to carboxylation capacity (V_cmax25, at a reference temperature of 25 ^∘C) and M_a, respectively. The coordination hypothesis provides predictions of V_cmax25: it is assumed that the light- and RuBisCO-limited assimilation rates under daytime conditions tend to equality, thus minimizing both maintenance respiration and the metabolic component of leaf N (Chen et al., 1993). An EEO-based hypothesis for the leaf economics spectrum provides predictions of M_a: it is assumed that the average net carbon gain by a leaf during its life cycle is maximized. The sum of the N in metabolic and structural components then determines the optimal leaf N content (N_mass). We set out (1) to analyse the contributions of climate, species identity, and phylogeny to leaf C_mass, N_mass, and their ratio; (2) to characterize geographic patterns in these traits along environmental gradients; and (3) to test the extent to which variation in N_mass among communities could be captured by the EEO principles outlined above. Using a data set comprising 1705 samples at 79 sites throughout China, we quantified the phylogenetic signal in species' traits and fitted a Bayesian mixed-effects model to partition individual trait variation into effects of climate, species identity, and phylogeny. We examined trait relationships with bioclimate variables (which improve on annual mean quantities by accounting for seasonality and latitude) and gridded data on soil C:N ratios by multiple regression.

2 Materials and methods

2.1 Trait and environmental data

Our analyses are based on trait data in the China Plant Trait Database version 2 (CPTDv2, Wang et al., 2022). The CPTDv2 contains morphometric, chemical, and photosynthetic leaf trait data on 1529 species at 140 sites representing the different biomes in China as well as climate information for each site. In CPTDv2, a stratified sampling strategy was consistently used at each site to ensure that the dominant species in each canopy layer were sampled (detailed in Wang et al., 2018) and to avoid bias from different sampling strategies. A total of 25 trees, 5 shrubs, 5 lianas or vines, and 5 understorey species (grasses and/or forbs) were sampled at each site. When the number of trees was less than 25 at a site, all the tree species were sampled and additional samples from the other life forms were supplemented up to a maximum of 40 species. Thus, the species sampled at each site can be regarded as a representative sample of the plant community and average trait values at each site. We extracted leaf carbon concentration (C_mass, %), nitrogen concentration (N_mass, %), leaf mass per area (M_a, g biomass m⁻²), stable carbon isotope ratios (δ¹³C, ‰) and V_cmax25 (µmol C m⁻² s⁻¹) data from the CPTDv2. Although M_a data are available from 124 sites, C_mass, N_mass, and δ¹³C data are only available at 79 (C_mass, N_mass) and 74 sites (δ¹³C). However, these sites are well distributed across the temperature and aridity gradients (Supplement Fig. S1). Although V_cmax25 data are available only at 32 sites from southwestern and northeastern China and along an elevational transect in the Gongga Mountains, there are 960 measurements from these sites.

Xu et al. (2021) and Wang et al. (2018) provided full details of the species sampled and trait measurements made at each site. M_a was estimated from measurements of leaf area and dry weight following standard protocols (Cornelissen et al., 2003). Leaf area was taken as the projected area of a randomly selected leaf, or leaflet for compound leaves, using a LiDE 220 Scanner (Canon Inc., Huntington, NY, USA). The dry weight was measured after oven-drying at 75 ^∘C for 72 h to constant weight. The average of three measurements made on leaves from different individuals was taken as the M_a value of one species at each site. We used a portable infrared gas analyser system (LI-6400; LI-COR Inc., Lincoln, NE, USA) to make leaf-gas-exchange measurements in the field. Terminal branches from the outer canopy were collected and re-cut under water immediately prior to measurement. The relative humidity and chamber block temperature were set close to that of the ambient environment at the time of measurement with a constant airflow rate (500 µmol s⁻¹). V_cmax was calculated from the light-saturated rate of CO₂ fixation at ambient CO₂ using the one-point method (De Kauwe et al., 2016) and adjusted to a standard temperature of 25 ^∘C with the Arrhenius equation (Bernacchi et al., 2001). Due to the time-consuming measurement of leaf-gas exchange, the photosynthetic traits of one sample were measured for each species. For each species at a site, leaf C content, N content, and δ¹³C were measured using pooled samples of leaves from at least three individuals of the same species with an isotope ratio mass spectrometer (Thermo Fisher Scientific Inc., Carlsbad, CA, USA). Carbon isotope ratios were used to calculate isotopic discrimination (Δ) and then to estimate the ratio of leaf-internal to ambient CO₂ partial pressure (χ) using the method of Cornwell et al. (2018) with a standard formula using the recommended values of a^′ and b^′ of 4.4 ‰ and 27 ‰, respectively (Cernusak et al., 2013; Farquhar et al., 1989):

$$\begin{array}{}\text{(1)}& \mathit{\chi }={\displaystyle \frac{\mathrm{\Delta }-a^{\prime }}{b^{\prime }-a^{\prime }}}.\end{array}$$

The bioclimate variables available for each site include an annual plant-available moisture index (α_p, an estimate of the ratio of annual actual evapotranspiration to potential evapotranspiration), the mean temperature of the coldest month (MTCO, ^∘C), the mean temperature during the thermal growing season, defined as the period with temperatures above 0 ^∘C (mGDD0, ^∘C), and the leaf-area-index-(LAI-)weighted photosynthetic photon flux density (I_abs, mol m⁻² s⁻¹) during the thermal growing season. The climate variables were interpolated to each site from 1814 weather stations in China using ANUSPLIN (Hutchinson and Xu, 2013). The LAI during the sampled month and year for each site from both data sets was extracted from the MODIS LAI product (MCD15A3H: https://modis.gsfc.nasa.gov/, last access: 3 February 2022) to provide a measure of canopy cover. We used the C:N ratio in topsoil (0–30 cm) extracted from the gridded soil data set of Shangguan et al. (2013) as an inverse index of soil fertility.

2.2 Trait prediction

The maximum capacity of carboxylation (V_cmax) was predicted using an EEO model based on the coordination hypothesis (Eq. 2), which states that plants coordinate RuBisCO-limited and light-limited photosynthesis rates to be equal under daytime conditions so that the available light is used without incurring futile maintenance costs (Prentice et al., 2014; Wang et al., 2017):

$$\begin{array}{}\text{(2)}& {\displaystyle }{V}_{\mathrm{cmax}}\approx {\displaystyle \frac{{\mathit{\phi }}_{\mathrm{0}}{I}_{\mathrm{abs}}(c_{\mathrm{a}}\mathit{\chi }+K)}{c_{\mathrm{a}}\mathit{\chi }+\mathrm{2}{\mathrm{\Gamma }}^{\ast }}},\text{(3)}& {\displaystyle }{\mathit{\phi }}_{\mathrm{0}}={\displaystyle \frac{\mathrm{0.352}+\mathrm{0.021}T-\mathrm{0.00034}{T}^{\mathrm{2}}}{\mathrm{8}}},\end{array}$$

where φ₀ is the intrinsic quantum efficiency of photosynthesis (µmol C µmol⁻¹ photon), which can be estimated for C₃ plants using Eq. (3) (Bernacchi et al., 2003); c_a is the ambient partial pressure of CO₂ (Pa); χ is the ratio of leaf-internal to ambient CO₂ partial pressure (Pa Pa⁻¹); K is the effective Michaelis–Menten coefficient of RuBisCO (Pa); Γ^∗ is the photorespiratory compensation point (Pa); and T is temperature (^∘C). We used mGDD0 as the temperature input.

M_a was predicted using an eco-evolutionary optimality model that predicts the relationship between M_a and leaf longevity, based on the assumption that leaves maximize net carbon gain during their life cycle (Wang et al., 2023). The predicted environmental effects on M_a differ between evergreen and deciduous species:

$$\begin{array}{}\text{(4)}& {\displaystyle }\begin{array}{rl}\ln \left(M_{\mathrm{a},\mathrm{de}}\right)& =\ln \left(f\right)+\ln \left(I_{\mathrm{abs}}\right)-\mathrm{0.052}T\\ & -\mathrm{0.27}\ln \left({\mathit{\alpha }}_{\mathrm{p}}\right)+\mathrm{2.65},\end{array}\text{(5)}& {\displaystyle }\begin{array}{rl}\ln \left(M_{\mathrm{a},\mathrm{ev}}\right)& =\mathrm{0.25}\ln \left(f\right)+\mathrm{0.5}\ln \left(I_{\mathrm{abs}}\right)-\mathrm{0.013}T\\ & -\mathrm{0.51}\ln \left({\mathit{\alpha }}_{\mathrm{p}}\right)+\mathrm{3.53},\end{array}\end{array}$$

where M_a,de and M_a,ev are the predicted M_a for deciduous and evergreen species, respectively, and f is the ratio of thermal growing-season length (days) to the number of days in the year.

Dong et al. (2017) proposed a model for N_area as the sum of components proportional to M_a and V_cmax25. A simple manipulation of this model gives:

$$\begin{array}{}\text{(6)}& N_{\mathrm{mass}}=a+{\displaystyle \frac{b{V}_{\mathrm{cmax}\mathrm{25}}}{M_{\mathrm{a}}}},\end{array}$$

where a (g N g biomass⁻¹) and b (g N s µmol C⁻¹) are empirical coefficients fitted across all species. To test whether nitrogen allocation varied within the canopy (Charles-Edwards et al., 1987), the random effect of binned LAI on the intercept (a) and slope (b) was tested using a mixed-effects model. The fitted values of a and b can be found in Supplement Table S1.

The C:N ratio was estimated as the ratio of predicted C_mass from Eq. (7) and N_mass from Eq. (6) (Fig. 6). We also calculated C:N ratios using observed C_mass and predicted N_mass (Fig. 7a) to check whether C_mass values influence the prediction of C:N ratios. Due to the lack of an existing theoretical basis to predict C_mass, we fitted the following linear regression using all observed trait data and three climate variables:

$$\begin{array}{}\text{(7)}& \ln \left(C_{\mathrm{mass}}\right)=\mathrm{3.06}\ln \left(I_{\mathrm{abs}}\right)-\mathrm{0.18}T-\mathrm{0.48}\ln \left({\mathit{\alpha }}_{\mathrm{p}}\right).\end{array}$$

2.3 Data analysis

Statistical analyses were carried out in R4.1.1 (R Core Team, 2021). The relative importance of C_mass and N_mass in controlling the leaf C:N ratio was evaluated using the relaimpo package (Groemping, 2006). Within- and between-site variability in traits was measured by the standard deviation (SD). Bioclimatic effects on leaf stoichiometry at the community level (i.e. with unweighted community-mean values as the data points) were examined using standard multiple (fixed-effects) linear regression (lm), and partial effects of each climate variable were visualized using visreg (Breheny and Burchett, 2017). To account for the potential effect of spatial processes on trait variation, we performed multiple regression on distance matrices (MRM) to separate the spatial and bioclimatic effects using the ecodist package (Lichstein, 2006). Phylogenetic analyses were carried out on all species (including 561 genera in 175 families and 57 orders). The phylogenetic trees were constructed and coloured with species-averaged trait values using the S.PhyloMaker and ggtree packages (Qian and Jin, 2016; Yu et al., 2017). Phylogenetic signal was calculated for each trait, using Pagel's λ, which measures the extent to which related species tend to have similar trait values. Pagel's λ varies from 0 to 1, indicating a low to high phylogenetic signal. It was calculated using the phytools package (Münkemüller et al., 2012; Revell, 2012). The significant values obtained indicate that values of these traits tend to be conserved within lineages. A Bayesian phylogenetic linear mixed model was applied at the species level (i.e. each occurrence of each species was treated as a data point), with species identity and phylogeny as random effects, using the MCMCglmm package (Hadfield, 2010). The model was repeated using three different phylogenetic hypotheses to account for the uncertainty in phylogenetic trees generated in these scenarios. Marginal (climate effects alone as fixed effects, without random effects) and conditional r² (with species and phylogeny as random effects) were compared (Nakagawa et al., 2017; Nakagawa and Schielzeth, 2013). The phylogenetic comparative method, which implicitly attributes overlapping effects of phylogeny and climate entirely to phylogeny (Westoby et al., 1995), was used to estimate trait variations explained by climate alone using the ape package (Paradis et al., 2004).

The temperature and moisture optima and tolerances of frequently occurring genera were calculated as follows (Meng et al., 2015). Each bioclimatic variable was binned and the mean abundance was calculated for the sites within each bin. The frequency distributions of abundance for each species and bioclimatic variable were obtained by selecting widths of the bins. The optimum was estimated as the average of the bioclimatic variable in the bins where a species was present, weighted by its mean abundance in the bins. Similarly, the tolerance was estimated as the abundance-weighted standard deviation of the bioclimatic variable. Frequently occurring species were defined as those that occurred more than 25 times for deciduous and more than 10 times for evergreen species.

We used trait-gradient analysis to access the intraspecific variation along the environmental gradient (Ackerly and Cornwell, 2007). The individual trait values were plotted against community-mean trait values which were strongly determined by external filters at a large scale, such as abiotic factors. By definition, the slope of the relationship between individual and community-mean trait values is unity. To characterize the plasticity (here including both genetic and phenotypic variation) of one species, the regression of individual trait values among one species against community-mean values represents its within-species ability to shift along the environmental gradient. The steep slope indicates plastic species and a large contribution of intraspecific trait variation to the overall trait gradient. Species that were measured at more than five sites were selected to estimate the slope.

Figure 1 Phylogenetic tree and signal for (a) leaf carbon content (%), (b) nitrogen content (%), and (c) carbon-to-nitrogen ratio. The red-to-blue colours of phylogenetic trees indicate low to high trait values. The phylogenetic signal is indicated by a statistic metric of Pagel's λ with its significance level (${}^{***}$ p<0.001). Orders with more than five species in the data set are labelled.

3 Results

3.1 Roles of phylogeny and species in stoichiometric variation

In general, related species tended to have similar stoichiometric traits; Pagel's λ was significant for all traits (Fig. 1). The mixed model produced higher conditional r² values when species and phylogeny were included as random effects, with species contributing 21 %–35 % and phylogeny contributing 16 %–18 % on average (Table 1). No significant relationships were found between leaf stoichiometric traits and the soil C:N ratio (Supplement Fig. S2).

Table 1 Regression models for each trait. The multiple linear regression (MLR) model was fitted at the community level, using only climate predictors. The Bayesian phylogenetic linear mixed model (BPLMM) was fitted at the species level. The marginal r² includes climate effects only; the conditional r² also includes species identity and phylogeny as random effects. The phylogenetic comparative method (PCM), also fitted at the species level, quantifies the variation attributed to climate alone, after factoring out effects of phylogenetic relatedness. The standard deviations of r² in BPLMM come from three different scenarios of phylogeny.

These findings are consistent with species turnover (taxonomic replacement) being a principal mechanism accounting for the observed trends in stoichiometry along environmental gradients. Figure 2 illustrates the turnover of major woody genera along the climatic gradients. The deciduous genera covered a wide range in moisture (α_p from 0.2 to 1) and light (I_abs from 7 to 16 mol m⁻² s⁻¹) but showed limited temperature tolerance ranges (Fig. 2a). Conversely, evergreen genera occupied a wide range on the temperature axis (from 6 to 21 ^∘C) but occurred only in wetter areas with α_p>0.6 (Fig. 2b). Growing-season (mGDD0) and coldest-month (MTCO) temperatures were positively correlated (not shown). The distribution of deciduous genera along the MTCO axis was similar to their distribution along the mGDD0 axis (Fig. 2c), whereas evergreen genera were more separated on the mGDD0 gradient than by MTCO – with the exception of Pinus, which showed a wide cold-tolerance range from around −30 to 0 ^∘C (Fig. 2d).

Figure 2 Optima and tolerances of major genera in climate space. α_p is a moisture index, mGDD0 is the mean temperature during the thermal growing season, MTCO is the mean temperature of the coldest month, and I_abs is the leaf-area-index-weighted photosynthetic photon flux density. Colours of circles represent the values of the leaf C:N ratio and the sizes of circles represent I_abs. The grey triangles are sampling sites. Also shown are abbreviated names of genera: Ac, Acer; Al, Allium; Ar, Artemisia; Cg, Caragana; Cx, Carex; Ls, Lespedeza; Ln, Lonicera; Pt, Potentilla; Qc, Quercus; Rb, Rubus; Sl, Salsola; Th, Thalictrum; Cm, Camellia; Cs, Castanopsis; Cy, Cyclobalanopsis; Fc, Ficus; Il, Ilex; Lt, Lithocarpus; Pn, Pinus; Rh, Rhododendron; Sm, Smilax; Sy, Symplocos.

Trait-gradient analysis showed that in addition to species turnover, intraspecific trait variation played a role in determining trait shifts at a regional scale (Fig. 3). The intraspecific slopes for C_mass, N_mass, and their ratio were calculated for 19, 19, and 42 species, respectively. Only 9, 8, and 16 of these species showed significant slopes. The intraspecific slopes for N_mass and C:N ratio ranged from 0.7 to 2.1 and 0.6 to 1.9, respectively. The slopes for C_mass ranged from 0.8 to 1.4 except for one species (Asparagus dauricus) that had a negative slope.

3.2 Leaf stoichiometric trait responses to climate

At the community level, climatic variables explained 10 %, 13 %, and 45 % of variation in C_mass, N_mass, and C:N ratio, respectively. At the species level, climatic variables explained 8 %, 3 %, and 18 % according to the mixed model. Smaller amounts of variation (2 %, 2 %, and 10 %) were captured by climate according to the phylogenetic comparative method (Table 1). MRM analysis also showed that trait variations were strongly explained by climatic factors but not significantly related to geographic distance – indicating that the purely spatial effect on trait values was weak (Supplement Table S2).

Stoichiometric trait responses to climate were generally similar in deciduous and evergreen species (Fig. 4). C_mass was significantly positively related to light and negatively related to moisture and growing-season temperature in both deciduous and evergreen species (Fig. 4a–c). N_mass significantly decreased with increasing light and moisture in both deciduous and evergreen species. N_mass also decreased with temperature in deciduous species but showed no significant relationship with temperature in evergreen species (Fig. 4d–f). The response of the leaf C:N ratio to climate was a combination of the C_mass and N_mass responses but was dominated by climate effects on N_mass. The leaf C:N ratio was positively related to light and moisture in both deciduous and evergreen species. It was also positively related to temperature for deciduous species but marginally negatively related to temperature for evergreen species (Fig. 4g–i).

Figure 3 Distribution of site-mean trait values and within-species regression lines for leaf stoichiometric traits. The black dots are individuals, and black dashed line is the regression of all individuals with slope equal to 1 by definition. The red lines represent significant regressions within species that were sampled at more than five sites.

Figure 4 Empirical partial relationships between leaf traits and climate. I_abs is the leaf-area-index-weighted photosynthetic photon flux density, α_p is a moisture index, and mGDD0 is the mean temperature during the thermal growing season. Yellow dots: deciduous species; green dots: evergreen species. The N_mass scale is inverted, so that the slopes of the regression lines in panels (a)–(c) and (d)–(f) should add up to the slopes in panels (g)–(i). The significance of each regression is listed in the panels. ${}^{***}$ p<0.001; ${}^{**}$ p<0.01; ^* p<0.05.

3.3 Eco-evolutionary optimality models for leaf traits

C_mass was relatively constant at different values of the leaf C:N ratio (grey lines, Supplement Fig. S3), while N_mass showed much greater variability. This pattern held for both deciduous and evergreen species. The analysis of relative importance showed that N_mass explains on average 90 % of variation in the leaf C:N ratio.

Leaf N_mass was positively related to its theoretical predictor ($V_{\mathrm{cmax}\mathrm{25}}/M_{\mathrm{a}}$) (Fig. 5). We found a significant LAI effect on the slope (b) and intercept (a) of this relationship, with the slope increasing and intercept decreasing towards greater LAI. The r² of Eq. (6) was improved from 0.14 to 0.21 at the species level after including the LAI effect. The optimality models captured 68 % and 53 % of the community-level variation in V_cmax25 and M_a, respectively (Fig. 6a, b). V_cmax25 was somewhat underestimated at most sites, with the largest bias when observed V_cmax25 was at alpine sites above 4000 m. M_a was distributed evenly near the 1:1 line, with the largest bias occurring at a semi-arid site with very high observed M_a. The optimality model, with LAI effect included, explained 21 % of N_mass variation using predicted values of V_cmax25 and M_a (Fig. 6c). The predicted leaf C:N ratios fell within the observed range in each PFT and outperformed fixed values prescribed in LSMs for most PFTs (Fig. 7). The prediction of the leaf C:N ratio using constant C_mass (45.6 %) was similar to that using observed C_mass.

Figure 5 The relationship between N_mass and $V_{\mathrm{cmax}\mathrm{25}}/M_{\mathrm{a}}$ along the leaf area index (LAI) gradient. Colour saturation represents different levels of LAI. Lines are separate regressions for sites within each LAI bin.

Figure 6 Optimality-based predictions versus observations of leaf traits at the site level. Grey lines are ordinary least-squares regressions. The black dashed line is the 1:1 line.

4 Discussion

We have demonstrated that across-site variations in leaf stoichiometric traits along climate gradients are driven mainly by species turnover and that an optimality-based model can predict 30 % of N_mass variation, highlighting the potential of applying EEO principles to leaf stoichiometry. Predicted leaf C:N ratios are within the range of observations. These findings provide a potential avenue for improving the representation of leaf stoichiometry in LSMs.

4.1 Climate effects mediated by compositional shifts and intraspecific variation

For C_mass and N_mass separately, the contribution of climate variables in the multiple regression was modest (10 %–13 %), while the mixed model attributed larger fractions of variation to species identity and phylogeny than to climate (Table 1). For the C:N ratio, climate explained 45 % of variation in the multiple regression, while the mixed model attributed similar fractions (19 %–20 %) of variation to species identity, phylogeny, and climate (Table 1). The phylogenetic comparative method (Table 1) attributed only 2 % of variation in C_mass and N_mass and 10 % of variation in the C:N ratio to climate alone. The trait-gradient analysis also indicated an important role of intraspecific variation in trait shift at a regional scale. These results are consistent both with strong phylogenetic control of leaf stoichiometry and with strong patterns of variation (especially for C:N ratio) in community-mean values determined to a substantial degree by environmental selection among species and clades characterized by different trait values (Liu et al., 2022), and they emphasize the contribution of within-species variation, which has often been neglected in community ecology (Violle et al., 2012).

Figure 7 Comparison of observed and predicted C:N ratios with target values in LSMs. The blue boxes represent predicted C:N ratios using observed C_mass (a) and using constant mean C_mass (b). The green boxes represent observed C:N ratios. Red crosses show target values in LSMs (adopted from CLM4).

Within-site variations – unconstrained by macroclimate – were usually larger than between-site variations (Supplement Fig. S4). This large within-site variability might explain why no significant effect of climate on leaf stoichiometry was detected in some previous regional studies (Yang et al., 2016; Zhang et al., 2017; Zhao et al., 2018). Zhang et al. (2019) showed a weak phylogenetic signal for the leaf C:N ratio evaluated by Blomberg's K. We used Pagel's λ due to its better performance and reliability with a large number of species (Münkemüller et al., 2012). The significant phylogenetic signals for leaf stoichiometric traits confirmed that species with similar evolutionary history tend to have similar leaf stoichiometry, indicating that leaf stoichiometric traits of extant species at a site may not remain adaptive under a changing environment (He et al., 2010; Li et al., 2021; Yang et al., 2016). It has been suggested that a high phenotypic plasticity of the leaf C:N ratio would be associated with a high mortality risk, supporting the idea that tight regulation of leaf stoichiometry within species helps to ensure plant survival (Luong et al., 2021). Within-site diversity may help communities to maintain their function in the face of climate variability and extremes.

Plant species may occupy different “biogeochemical niches” to ensure the full use of available resources and avoid competition (Sardans and Penuelas, 2014; Sardans et al., 2021). At the community level, climate variables captured more of the observed leaf stoichiometric variations, due to the averaging of data from co-occurring species and intraspecific variability (Vallicrosa et al., 2021). Systematic variation in community-mean leaf stoichiometric traits along climate gradients can be achieved through progressive species replacement at a macroclimatic scale and intraspecific trait variability at a regional scale (Liu et al., 2019; Yang et al., 2016).

The distributions of common deciduous genera were shown to be more sensitive to moisture, while the distributions of evergreen genera were mainly driven by temperature. This distinction may be related to the different adaptation strategies represented by differences in leaf longevity (LL). Kikuzawa et al. (2013) indicated that temperature is the best predictor of LL for evergreen species, while a consideration of an additional moisture factor was expected to improve the explanatory power of climate for LL in deciduous species. According to Kikuzawa's optimality model, LL of evergreen species is higher at low temperatures, in order to compensate for low total carbon gain during the short growing season. For deciduous species, however, LL should not exceed the length of the growing season – which can be affected by moisture as well as temperature in semi-arid and arid areas. Thus, our study suggests that climate shapes leaf stoichiometric variation at a macroclimatic scale via environmental selection among taxa and emphasizes the neglected role of phenology in biogeochemical cycles (He et al., 2006; Vallicrosa et al., 2021; Xiong et al., 2021).

Although some studies have shown an important role of soil fertility in determining plant stoichiometry, published studies have shown inconsistent results (Fang et al., 2019; Fyllas et al., 2009; He et al., 2010; Ordoñez et al., 2009; Xiong et al., 2021). Soil fertility as indexed by the soil C:N ratio had no significant effect on leaf stoichiometry in our analysis, indicating a decoupling of soil and leaf stoichiometry (Delgado-Baquerizo et al., 2017; Elser et al., 2010). Plant–soil interactions may affect whole-plant stoichiometry nonetheless, through effects on C allocation to different tissues. Allocation of N to leaves shows stronger homeostasis than in other tissues, possibly as a consequence of the need to maintain the crucial functions of photosynthesis and leaf respiration; the stoichiometry of other tissues may adjust to soil conditions in order to support leaf-level function (Chen et al., 2013; Delgado-Baquerizo et al., 2017; Zhang et al., 2017). Uncertainty in our soil fertility data may was inevitably introduced due to our reliance on a gridded soil map (Shangguan et al., 2013). More studies including in situ soil measurements are needed to more comprehensively investigate the effect of soil properties on plant stoichiometry.

4.2 Trait responses reflect plant strategies

Leaf stoichiometry integrates traits that reflect different plant functions, resulting in a potentially complex response to climate. In contrast with many previous studies, we have considered (and found significant effects of) functionally significant bioclimatic variables, including light, in leaf-level stoichiometry. Our analyses indicate general relationships that are quite similar between evergreen and deciduous plants.

Higher leaf C_mass was observed in cold and dry areas with high radiation (Fig. 4). Chen et al. (2021) found that leaf C_mass is positively related to vein density, which relates to the efficiency of water transport. At strong light, plants tend to have a higher photosynthetic rate, requiring more water for transpiration to maintain open stomata – which could be achieved by high carbon investment in venation (Sack and Scoffoni, 2013). In dry areas, high vein density is a common adaptation to drought, allowing plants to respond quickly to available water for carbon fixation and to keep leaves cool in the face of high air temperature (Scoffoni et al., 2011; Yao et al., 2021). Meanwhile, plants may accumulate nonstructural carbohydrates to adjust osmotic potential and avoid leaf desiccation (Bartlett et al., 2014). The leaf C_mass response to temperature as observed here, however, is the opposite of some previous reports (Ma et al., 2018; Xing et al., 2021). Global analysis showed an overall positive response to temperature with a range from −10 to 30 ^∘C, whereas leaf C_mass decreased when the mean annual temperature was lower than 20 ^∘C (Ma et al., 2018). M_a, which is positively related to leaf C_mass (Xing et al., 2021), is generally negatively correlated with temperature (Wright et al., 2004). A higher starch concentration is observed at low temperatures, due to conditions that allow photosynthesis but not growth (Hoch and Körner, 2012). This suggests that the leaf C_mass response to temperature may not be monotonic, owing to different functions dominating at the extremes.

The components of leaf N variation adapt to the climate in different ways (Dong et al., 2017; Peng et al., 2020; Xu et al., 2021). N_mass is also constrained by the trade-offs inherent in the leaf economics spectrum. Leaves with high M_a (and LL) have low photosynthetic rates per unit mass and low nutrient contents by mass (Wright et al., 2004). Thus, climate drives N_mass variation both directly and indirectly. Moisture has a negative effect on N_mass (Yang et al., 2016; Zhang et al., 2019; Zhao et al., 2018). It has been reported that N-containing compounds (such as amides) accumulate in plants in order to adjust osmotic pressure under drought (Raggi, 1994). In addition, N-rich leaf defence compounds increase towards more arid climates, at the expense of C-based defences such as spines and thorns (de Oliveira et al., 2020; Ghimire et al., 2017; Meloni et al., 2012). In the existing optimality model N_area should be positively related to radiation, since light has a positive effect on M_a and V_cmax25 (Smith et al., 2019; Wang et al., 2023). However, when N is expressed on a mass basis (N_mass), light has a negative effect, suggesting a lower sensitivity of V_cmax25 to light than M_a. The negative effect of temperature on N_mass has been observed across vegetation types (Han et al., 2005; He et al., 2008; Tang et al., 2018; Weih and Karlsson, 2001). This is consistent with the hypothesis that more nutrients are required to compensate for low enzyme activity at low temperatures (Reich and Oleksyn, 2004). Although the leaf C:N ratio response to climate is a combination of the responses of both C and N, it is dominated by the variation in N. Thus, understanding of N_mass variation should help elucidate variation in C:N ratios (Reich, 2005). Positive effects of temperature and moisture on leaf C:N ratios have also been observed in previous studies, implying higher N use efficiency in hot and wet areas (Fang et al., 2019; Zhang et al., 2019).

4.3 Leaf nitrogen content predicted by optimality models

C:N ratios couple C and N cycling, thus influencing the estimation of carbon assimilation and plant growth in LSMs (Wang et al., 2010; Zaehle et al., 2014). Fixed leaf C:N ratios assigned to PFTs, as for example in CLM4 and ED2.1, may result in inaccurate representations of this coupling (Bonan and Doney, 2018; Lawrence et al., 2011; Medvigy et al., 2009). Although model outputs such as ecosystem responses to elevated CO₂ are more consistent with observations in models where flexible C:N ratios are allowed (Lawrence et al., 2019; Meyerholt and Zaehle, 2015), large differences between models persist (Du et al., 2018).

Meyerholt and Zaehle (2015) highlighted the potential of optimality theory to improve the representation of N cycling in LSMs. Caldararu et al. (2020) showed that models that implement dynamic leaf stoichiometry schemes based on EEO principles can perform better than those with fixed-PFT schemes. Here we have shown that leaf N_mass covariation with V_cmax25 and M_a as predicted by EEO principles can provide further insights. Since N_mass is the key to determining the leaf C:N ratio, given the relative constancy of C_mass (Reich, 2005), we focused on the predictability of N_mass. The variation in leaf N per unit area (N_area) can be represented as the sum of two components, proportional to leaf mass per area (M_a) and the maximum capacity of carboxylation at 25 ^∘C (V_cmax25) (Dong et al., 2017), and now both M_a and V_cmax25 can be predicted from EEO principles (Smith et al., 2019; Wang et al., 2023; Xu et al., 2021). Community-level variations in M_a, V_cmax25, and N_area can indeed be largely captured (r²=0.53, 0.68, and 0.62, respectively) using climate variables as predictors. We also showed a tendency for the relationship between N_mass and the ratio $V_{\mathrm{cmax}\mathrm{25}}/M_{\mathrm{a}}$ to become steeper with increasing LAI. This finding is consistent with N redistribution within the canopy, as an acclimation to light conditions that maximizes total carbon gain (Hirose and Werger, 1987; Niinemets et al., 2015). The strong vertical light gradient in high-LAI canopies implies a large advantage for optimized N distribution, in contrast with more open canopies (Field, 1983). The same model framework as N_area, with this additional LAI effect included, showed good predictive skill for N_mass (and better than that of Boonman et al. (2020), obtained using an ensemble modelling approach) based on climate. However, our predicted N_mass was constrained within a narrow range, despite the well-captured variations in M_a and V_cmax25. The predicted N_mass in tropical forest with high LAI was systematically underestimated due to the low intercept (Supplement Table S1). We recognize that our method to predict N_mass may overlook additional functions of N in leaves, such as chemical defences, perhaps causing greater variation than predicted. This requires further investigation. Some species in this study apparently adjusted their leaf stoichiometry along major environmental gradients, possibly via genetic adaptation over multigenerational timescales. Due to the lack of intraspecific data within communities, we could not assess the degree of variation among conspecific plants in the same environment. Intraspecific variation within communities may, however, increase functional diversity and promote species coexistence (Westerband et al., 2021) and potentially provide a buffer against climatic variation and change (Ahrens et al., 2021). Further studies are needed to better understand intraspecific trait variation (Moran et al., 2016) in order to assign appropriate timescales for the dynamic responses of traits to environmental changes in Earth system models.

Predicted C:N ratios, whether using observed or constant C_mass, lie within the range of observed data, supporting the dominant role of N_mass in driving leaf C:N ratios (Fig. 6). The target (PFT-specific) values used in several LSMs such as CLM4, ORCHIDEE, and YIBs (Fig. 7) are based on data sets nearly 20 years old and fail to represent continuous trait variations that can now be inferred from much larger data sets. Our EEO-based approach thus suggests a way forward to improve the dynamic representation of leaf stoichiometry in LSMs.

5 Conclusions

This study shows that the leaf C:N ratio is mainly driven by mass-based leaf nitrogen content which can be estimated via the sum of metabolic and structural components of leaf nitrogen using eco-evolutionary optimality-based models. This provides another perspective to improve dynamic representation of stoichiometry in Earth system models. The variations in leaf stoichiometric traits at the individual level are mainly controlled by species identity and phylogeny; thus, the shift in leaf stoichiometry variations at the community level along climate gradient is achieved via species turnover and intraspecific variability. This allows the prediction of community-mean values of leaf stoichiometric traits using EEO-based models. We show that the coefficient representing nitrogen allocation to metabolic and structural components is related to leaf area index, which highlights the importance of nitrogen allocation in its prediction. The unexplained variation in leaf nitrogen content may attribute to other unclear physiological processes, which requires further effort to improve the prediction of the leaf C:N ratio.