Environmental versus phylogenetic controls on leaf nitrogen and phosphorous concentrations in vascular plants

Domain: earth_science
Published: 2024-06-24
Editorial outcome: accept

Abstract
Global patterns of leaf nitrogen (N) and phosphorus (P) stoichiometry have been interpreted as reflecting phenotypic plasticity in response to the environment, or as an overriding effect of the distribution of species growing in their biogeochemical niches. Here, we balance these contrasting views. We compile a global dataset of 36,413 paired observations of leaf N and P concentrations, taxonomy and 45 environmental covariates, covering 7,549 sites and 3,700 species, to investigate how species identity and environmental variables control variations in mass-based leaf N and P concentrations, and the N:P ratio. We find within-species variation contributes around half of the total variation, with 29%, 31%, and 22% of leaf N, P, and N:P variation, respectively, explained by environmental variables. Within-species plasticity along environmental gradients varies across species and is highest for leaf N:P and lowest for leaf N. We identified effects of environmental variables on within-species variation using random forest models, whereas effects were largely missed by widely used linear mixed-effect models. Our analysis demonstrates a substantial influence of the environment in driving plastic responses of leaf N, P, and N:P within species, which challenges reports of a fixed biogeochemical niche and the overriding importance of species distributions in shaping global patterns of leaf N and P.

Introduction
As fundamental elements for vascular plants on Earth, nitrogen (N) and phosphorus (P) are important in controlling photosynthesis, growth, and ecological functions of plants1–4. The relationship between leaf N and P concentrations is a key plant characteristic reflecting leaf economics5–7 and has been used to interpret ecosystem nutrient limitation8–11, carbon (C)- and N-cycle interactions under global change12–14, or to link macroecology and biogeography with trait-based functional ecology15–17. Understanding the variation in leaf N and P stoichiometric patterns and their underlying controls is therefore crucial for predicting responses of terrestrial ecology and biogeochemical cycles to environmental change.

Variations in leaf N and P stoichiometry along geographic gradients are pervasive18–28. For example, leaf N and P concentrations generally increase, but N:P ratios decrease from the equator to the cooler and drier mid-latitudes22. Such patterns have been interpreted in the light of different hypotheses regarding dominant controls.

Previous hypotheses have emphasized either the continental-scale variations of soil age and fertility29; the interactive effects of leaf N and temperature on biochemical reactions, including photosynthesis30; the effect of soil temperatures and chemistry on nutrient mineralization and availability to plants22; or the role of species distribution21 in shaping global patterns of leaf N, P and N:P. Overall, these contrasting viewpoints provide contentious interpretations for global leaf N and P stoichiometry variation and take conflicting viewpoints regarding the importance of phylogenetic vs environmental factors in shaping these stoichiometric patterns.

Recently, studies21,31 found an overriding effect of phylogeny and a vanishingly small influence of environmental variables on leaf nutrient concentrations based on linear mixed-effect models. Results were interpreted to support the importance of species distribution in shaping global-scale leaf nutrient patterns and the so-called “Biogeochemical Niche Hypothesis”21. This hypothesis posits that each species is characterized by a fixed leaf stoichiometry (low within-species variability) that matches environmental constraints and that geographical patterns in leaf nutrient distributions arise through environmental filtering of species occurrences. However, this interpretation of an overwhelming phylogenetic control on leaf N and P stoichiometry contrasts with previously documented influences of climatic and edaphic variables on leaf N and P stoichiometry7,19,32.

The conflicting attributions of observed variation in leaf N and P stoichiometry to phylogenetic vs environmental variables are related to an inherent methodological challenge. The separation of these variables is usually undermined by their lack of independence. The distribution of plant species is largely driven by the abiotic environment33,34. Yet, species do occur over a certain range of environmental conditions. To what extent the environment drives phenotypic plasticity or genetic adaptation in leaf N and P stoichiometry also within species remains challenging to detect but is informative for testing the Biogeochemical Niche Hypothesis. Linear mixed-effect models (LMMs) have been widely employed for separating phylogenetic and environmental effects on leaf traits21,35, motivated by their suitability to model structured data and their ability to control for phylogenetic effects and species identity as random terms, implicitly assuming that they are unrelated to the environment and given precedence over the latter in model fitting. More recently, tree-based statistical learning methods, for example, random forest models (RF), have been shown to be suitable for modeling leaf N and P36,37. These models, too, provide a natural way to simultaneously account for environmental (continuous) and phylogenetic (categorical) information. However, the implications of methodological choices for separating environmental vs phylogenetic variables so far have not been explicitly considered.

In view of these conflicting reports and methodological challenges, the question arises to what extent large-scale leaf N and P stoichiometric patterns are a reflection of different species (with their relatively fixed leaf nutrient stoichiometry) occurring at different sites along environmental gradients, and to what extent plasticity and genetic adaptation, driven by the environment, drive variation within species and contribute to large-scale patterns of leaf nutrient stoichiometry (Fig. 1).

To address this question, we compiled a global dataset of 36,413 paired leaf N and P concentrations per unit leaf mass of vascular plants, complemented by a comprehensive set of climatic, edaphic, and other environmental variables extracted from global datasets. Using these data, we first compared the power of linear regression models (LMs), LMMs, and RF models in explaining different components of variation in the data (variation across sites, within species, and across species). To enable comparability and facilitate the interpretation of our results in the context of the published literature, we specified LMMs to reflect the methodological choices of previous studies21,31. Then, we performed a trait gradient analysis (TGA)38–40 to address the second question: does spatial variation in leaf N and P stoichiometry arise predominantly from species distribution and their respective stoichiometry (biogeochemical niche), or is there substantial within-species variation? We tested whether, according to the Biogeochemical Niche Hypothesis21, the slopes of species-level trait gradient regressions were flat, as shown in Fig. 1. Alternatively, stoichiometric plasticity will yield positive slopes in the TGA. Perfect stoichiometric plasticity is indicated by slopes distributed around 1.

Our results show that variations within species are similarly strong as variations between species and that they are clearly influenced by the environment. These findings fill the gap of distinguishing and quantifying the role of species identity vs environmental controls on leaf nutrient stoichiometry and indicate that the previous Biogeochemical Niche Hypothesis, with its interpretation of an overriding and almost exclusive effect of phylogeny on leaf N, P, and N:P, should be revised.

Results
Variables selection and effects
We started by identifying the most important environmental variables for explaining variation in leaf N, P, and N:P. Reduced predictor sets, specific for leaf N, P, and N:P, respectively, enabled an improved model performance compared with models that included all 45 predictors (Fig. 2a–c) and were used for all subsequent analyses. In LMM models (Fig. 2d–f), N-deposition (ndep) had the strongest effect on leaf N and leaf N:P variation within species (both positive). The temperature of the coldest month (tmonthmin) had the strongest effect on leaf P (negative). Atmospheric CO2 concentrations (co2), ndep, and tmonthmin were among the most important predictors for leaf N, P, and N:P. Soil variables were only selected among the most important variables for leaf N (aluminum saturation of the soil solution, ALSA) and for leaf P (soil texture, measured by the water holding capacity class, AWC_CLASS).

Contrasting model performances
The random forest (RF) models, fitted to site-level aggregated data (mean across all observations by site), with the selected subset of environmental variables as predictors, achieved an R2 of 0.46, 0.34, and 0.34 for leaf N, P, and N:P, respectively, in contrast to an R2 of 0.17, 0.19, and 0.19 in LMs (Fig. 3a, d, g). R2marg, measuring the proportion of variation explained by fixed (environmental) variables in LMMs, fitted to the full data, was only 0.04, 0.05, and 0.09 for leaf N, P, and N:P, respectively. In contrast, the proportions of variation explained by species identity were 0.68, 0.63, and 0.44 (intraclass correlation coefficient, ICC), respectively, for leaf N, P, and N:P ratio (Fig. 3c, f, i).

When RF models were fitted to the full data, environmental variables explained a larger proportion of the variation than they did in LMMs, namely 0.13, 0.26, and 0.16 vs 0.04, 0.05, and 0.09, respectively. A similar contrast in the predictive power of environmental variables in RF and LMs is seen with models fitted to the modified data that contained only within-species variation (Fig. 3b, e, h). Here, RF models achieved an R2 of 0.29, 0.31, and 0.22 for leaf N, P, and N:P, while LMs achieved an R2 of 0.01, 0.03, and 0.07, respectively.

The power of environmental variables vs species identity in LMMs is subject to methodological aspects of the model-fitting procedure. By design, random factors are fitted with priority and “absorb” potentially shared effects with environmental variables. This is reflected also when comparing LMMs with LMs that are specified with distinct orders of the predictors (Supplementary Table 3). Subject to that order, species identity, and environmental variables contribute different sums of squares in LMs (Supplementary Table 3). Based on LMs, the shared effect of species identity and sites was dominant for leaf N, P, and N:P, explaining more than double of the variance explained by their separate effects (Supplementary Fig. 1).

Trait gradient analysis
The trait gradient analysis showed patterns similar to RF for leaf N, P, and N:P (Fig. 4). For most species, there was considerable variation of leaf N, P, and N:P within species across sites and this variation paralleled the variation in mean leaf N, P, and N:P of multiple species recorded at respective sites. The most frequent slope, i.e., the mode of the density distribution of slopes (Fig. 4d), is at 0.98 for leaf N, 0.80 for leaf P, and 1.01 for leaf N:P. A peak of the distribution of slope values close to unity indicates predominant plasticity or within-species genetic variation (henceforth referred to as ‘intra-specific variation’) of the respective trait along environmental gradients. However, for leaf N, a substantial number of species exhibited intra-specific variation with slopes < 0.5. For leaf P, the most common degree of intra-specific variation appears somewhat smaller than the most common degree of intra-specific variation in leaf N and N:P. For leaf N, the distribution of slopes appears to be broader than for leaf P and N:P.

Ranges of site-level mean leaf N, P, and N:P, normalized by their respective overall mean, along which species occurred, tended to be the smallest for leaf N:P (median: 0.27), followed by leaf P (median: 0.31), and the largest for leaf N (median: 0.36, Fig. 4e). We found no correlation between normalized ranges and slopes, neither for N, nor for P, and no correlation between slopes for N and slopes for P (not shown). However, ranges in N and in P were positively correlated (Pearson’s r = 0.54, Supplementary Fig. 3).

We found clear differences in the degree of intra-specific variation across species (Fig. 4f, g). For example, Pinus sylvestris and Quercus petraea appear to have relatively static leaf N, P, and N:P. In contrast, Betula pendula and Picea abies appear to be plastic (or exhibit intra-specific genetic variation) in leaf N, P, and N:P. Others species differ in their intra-specific variation between leaf N and leaf P. For example, P. sitchensis and Q. ilex have plastic (or genetically variable) leaf N, but relatively static leaf P. A substantial portion of species is “super-plastic”, with slopes > 1. This may reflect a strategy to overexpress a response to the environment through enhanced sensitivity to environmental variations than what would be expected from patterns in the site-level mean data.

Other metrics of the importance of species identity
As indicated by the species variation decomposition, roughly half of the variations in leaf N and P arose within species, while the remainder was linked to variations across species (Fig. 4f). Among the three traits investigated, leaf N appeared to be most strongly linked to species identity, whereby interspecific variations explained 60% of overall variations (R2across = 0.60), while 40% of variations arose within species. The importance of species identity was weaker for leaf P (R2across = 0.53) and leaf N:P (R2across = 0.42).

Methods
Leaf N and P data
The dataset used here is extended from Tian et al.82—a compilation of matched leaf N and leaf P concentrations (both in units of mg g−1 based on dry leaf biomass) and mass-based leaf N:P measurements from the literature and large-scale field investigations in China18,46,47,83,84, respectively, from the TRY Plant Trait Database85, and from the International Co-operative Programme on Assessment and Monitoring of Air Pollution Effects on Forests, ICP Forests21,86. N:P values higher than 70 g N (g P)−1 (0.04% of the data) were removed from the dataset. All sites had only one sampling date. In total, the dataset contains 36,413 individual records, collected from 7549 distinct sites and 3625 distinct species in 1383 genera and 203 families. Sites cover a wide climatic gradient and all major biomes (Supplementary Fig. 2). Species names are homogenized following Flora of China (http://frps.eflora.cn/), Useful Tropical Plants (http://tropical.theferns.info/), Australian Native Plants (https://www.anbg.gov.au/index.html), Angiosperm Phylogeny Website (http://www.mobot.org/MOBOT/research/APweb/)87, and The Plant List (www.theplantlist.org; accessed March 2021).

Environmental data
All plots from which data were obtained were georeferenced (WGS84 standard). Using information on longitude, latitude, and elevation, we complemented our dataset by extracting data for a set of environmental variables from maps with global coverage. The geographical position of sites was verified and improved using original publications and Google Earth (https://earth.google.com/web/). To complement missing records of elevation, we extracted information for respective site locations from ETOPO188. All data extraction was done using the ingestR package89. In total, we used 45 environmental variables, representing climatic, edaphic (soil-related), and other environmental variables (Supplementary Tables 1–2).

We computed 12 climatic variables (Supplementary Table 1) based on WATCH-WFDEI climate data90, covering the years from 1979 to 2012, and down-scaled these climatic variables based on high-resolution WorldClim climatology91. In our dataset, sampling dates of leaf N and P spanned from 1935 to 2015. We used climatic time series here to account for long-term changes in climate in combination with known leaf N and P sampling dates. For sampling dates before 1979, we considered the mean climate for 1979–1988. Additionally, the annual mean ratio of actual over potential evapotranspiration (alpha), mean annual total evapotranspiration (AET), and the aridity index of annual mean precipitation over potential evapotranspiration (AI) were estimated using the SPLASH land water balance model and potential evapotranspiration based on the Priestly-Taylor Equation92. We additionally included model-based estimates of the maximum rate of carboxylation normalized to 25°C (Vcmax25), the electron transport for ribulose-1,5-bisphosphate (RuBP) regeneration (Jmax25), and a multi-day average stomatal conductance (gs). These were obtained from point-scale simulations of the P-model93, using its implementation in the rsofun R package94, predicting leaf-level acclimation of photosynthetic traits (Vcmax25, Jmax25, gs) based on optimizing the trade-off between carbon gain and water loss95. Due to the intrinsic96 and widely observed97 link between Vcmax25, Rubisco, and leaf N and P contents, we used these estimates here as predictors for leaf N, P, and N:P, thus attempting to account for simultaneous effects of multiple climatic variables (air temperature, VPD, elevation, irradiance, CO2) on leaf nutrient stoichiometry. As an alternative aridity-related predictor, we included an estimate of the 80-year maximum cumulative water deficit (CWDX80) from ref. 98. In total, we used 19 climate-related variables (Supplementary Table 1).

Soil properties related to texture and fertility were extracted from different digital global soil maps. Specifically, several properties related to soil structure, texture, and ion exchange capacity, as listed in Supplementary Table 2, were extracted from the Harmonized World Soil Database v1.2 (HWSD, https://www.fao.org/soils-portal/soil-survey/soil-maps-and-databases/harmonized-world-soil-database-v12/)99. Aluminum saturation (ALSA), organic carbon content (ORGC), total nitrogen content (TOTN), and soil C:N ratios were extracted from the harmonized dataset of derived soil properties for the world (WISE30sec, https://data.isric.org)100. Soil phosphorus concentration using the Bray method (PBR), total phosphorus concentration (TP), and potassium (TK) concentration were extracted from The Global Soil Dataset for Earth System Modeling (GSDE, http://globalchange.bnu.edu.cn/research)101. To keep the soil properties comparable, the standard for soil layers was set around 30 cm (i.e., 0–30 cm for HWSD, 0–40 cm for WISE30sec, and 0–28.9 cm for GSDE).

Using information on site-specific sampling dates, we complemented the dataset with atmospheric CO2 measurements from Mauna Loa Observatory102, averaged over respective years, uniformly for all sites (assuming globally well-mixed concentrations). Also using information of sampling dates, estimates of dry and moist atmospheric reactive N-deposition were extracted from outputs of global simulations of atmospheric chemistry over the historical period by ref. 103. Attempting to account for the small-scale redistribution of soil nutrients along the hillslope and local climatic effects mediated by the landscape position of the measurement plot, we included the Compound Topography Index (CTI)104 as an additional predictor.

Statistical models
We started our analyses by selecting a subset of the most important predictors for leaf N and P concentrations and leaf N:P using a recursive feature elimination. Relying on a limited set of predictors reduces collinearity among them and the potential for overfitting in subsequent model-based analyses. Starting from models that included all 45 available predictors, we iteratively removed the single predictor that led to the smallest decrease in the R2 of observed versus modeled values, determined from a five-fold cross-validation. We retained a final set of predictors considering the R2 determined on a five-fold cross-validation. The feature elimination was based on random forest (RF) models and was performed on data aggregated to site-level community-weighted mean leaf N and P concentrations, and leaf N:P. The algorithm’s hyperparameter mtry, specifying the number of predictors considered at each split in individual decision trees, was set to mtry=(K–1)/3, where K is the total number of predictors retained at the respective step of the feature elimination. The minimum node size (min.node.size, controlling the depth of decision trees) was set to 5 at all steps of the feature elimination and was chosen based on a prior hyperparameter search considering the root mean square error determined from a five-fold cross-validation. RF models were fitted using the ranger105 and caret106 libraries in R.

With the selected subset of predictors for each of the three target variables, we then performed a model-fitting comparison to cross-compare different modeling approaches (linear mixed models, LMM, and linear regression models, LM, vs RF) and the role of data aggregation. Models were fitted with the subset of selected predictors (i) to the full, species-level data, (ii) to data aggregated to the site-level considering community-weighted means of all variables, and (iii) to modified data that contained only within-species variations. For the latter, we modified leaf N and P concentrations, and N:P values by subtracting species-mean leaf N, P, and N:P values from all values of the respective species and considered these modified values as targets for modeling. RF models were fitted to all three types of data. For RF models fitted to the full data (i), information about species, genus, and family identity was one-hot encoded. Hyperparameters min.node.size and mtry were tuned for each model separately, considering the mean root mean square error determined across five cross-validation folds. Reported R2 values were quantified as the mean across five cross-validation folds. We quantified the proportion of variation in the full data explained by environmental variables alone by taking the difference between RF models fitted with environmental variables and species, genus, and family identity as predictors, and RF models with only species, genus, and family identity as predictors.

LMMs were fitted to the full data with the subset of environmental predictor variables as fixed effects, considering species identity as a grouping variable for random intercepts, and using the nlme library107 in R. Before fitting LMMs, data were Yeo-Johnson-transformed. Partitioning of the variation explained by environmental variables (fixed variables) vs species identity (grouping factor for random offsets) was done using the R package performance108. We quantified the marginal R2 (R2marg) to estimate the proportion of variation explained by the fixed effects — representing environmental variables. The intraclass correlation coefficient (ICC) was calculated to quantify the proportion of variation explained by species identity. For data aggregated to the site-level and the modified data that contained only within-species variation, we fitted linear regression (ordinary least squares) models (LM), using the same subset of selected predictors. Note that LMs were fitted to those data because across-species variations, accounted for by random factors in LMMs, are removed by design at the site-level aggregated data and in the modified data that contains only within-species variations.

We additionally performed an analysis of variance (ANOVA) of species identity and environmental variables, investigating its dependency on including species as a random factor in LMMs and on the order of specifying species as a “fixed” factor in LMs. This was done by comparing ANOVA tables resulting from fitting species identity as fixed effect before or after environmental variables, using the lm function in base-R. Additionally, we separated individual and shared effects of site and species based on linear models where the respective predictors were fitted in different orders (Supplementary Fig. 1).

Trait gradient analysis
To investigate variations of leaf N and P stoichiometry within vs across species, a trait gradient analysis was performed, following refs. 38–40. We fitted ordinary least squares (OLS) regressions to species-specific relationships between individual, species-level data points of leaf N (P, N:P) and the site-mean leaf N (P, N:P) of the site belonging to the respective data record. That is, a data point xij recorded from species i at site j was regressed against xj, the site-mean leaf N (P, N:P) for site j. Note that site-mean values are aggregated from data of all species sampled at the respective site. The regression slope (bi) of the regression of xij vs xj indicates the degree of intra-specific plasticity across sites with different environmental conditions, characterized by the site-level mean leaf N, P, and N:P, respectively40. A regression slope of zero can be interpreted as a non-plastic behavior (Fig. 1). A regression slope of one indicates an average and a regression slope > 1 an above-average plastic response of leaf nutrient concentrations and stoichiometry within a species across sites with different community-mean concentrations and stoichiometry. We also quantified the species-specific range along the horizontal axis (Ri), given by values xj. We defined ranges as the difference between 1% and 99% quantiles to reduce the influence of outliers. Species-specific range values were normalized by the overall range of site-level means xj to make ranges comparable across the three different traits investigated here. Normalized ranges are interpreted here as an indicator of the range of site conditions under which the respective species occurs. For the TGA, we considered only data from species that were recorded in at least five sites and used only data from sites where at least five different species were sampled. This filtering retained 3385 data points from 372 species. For fitting linear regression models, all data were log-transformed before performing the trait gradient analysis to improve the normality of the residuals. Normalized ranges were quantified based on original, not log-transformed data.

We also conducted a “species variation decomposition”, where we quantified and compared the coefficient of determination from comparing observed (unmodified) leaf N, P, and N:P values with modified values. Modifications followed two alternative assumptions, reflecting hypotheses in Fig. 1. The first assumes that variations arise only across species (no within-species variations), and values xi,k of species i and observation k are replaced by the respective species’ mean, x(i). The second assumes that variations arise only within species and each record is “normalized” such that the resulting species mean is equal to the global mean x(i).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The present study appears to be an important work. However, it is uncertain whether it\nrepresents a substantial advancement in the field and if the methods employed sufficiently\nsupport the arguments made in this paper. Specific concerns are as follows. Please consider\naddressing these concerns to enhance the clarity and depth of your research findings in the\npaper.\nSpecific comments:\nLines 44-48: What is the biological significance of comparing Random Forest models, trait\ngradient analysis, and linear models? Do these comparisons reflect specific biological\nmeanings or mechanisms?\nLines 90-92: This statement is difficult to comprehend. What is meant by \"cannot be\nindependent\"? In the natural world, biological organisms and their environments are\ninherently adaptive processes, making complete independence challenging. Does it refer to\nthe lack of independence between species identity and historical or current environmental\nfactors? Linear mixed effect models can also analyze the interaction between species\nidentity and environment.\nLine 93: What about the interaction between environmental controls and species identity?\nIt might not be an absolute dichotomy; their interactions could vary with changing\nenvironmental conditions or plant functional groups.\nLines 103-105: Why compare these three methods? Are there no better alternatives? Do the\ndifferences in these statistical methods hold any biological significance? Why compare\nvariation within species and between species? Does within-species variation change with\nplant developmental stages, age, or individual size?\nFigure 2: Are there issues of collinearity among the predictors? How does collinearity affect\nthe results of the three methods? How many predictors are used in the models? If there are\nnumerous predictors, could the models' abilities to handle complex data differ, potentially\naffecting the results?\nLines 241-242: How can the comparability of variation within and between species be\nensured?\nLines 249-250: What distinguishes \"variation\" from \"plasticity\"? What is the significance of\nwithin-species variations?\nLines 305-309: Linear mixed effect models can account for nonlinear relationships and\ninteractions, although Random Forest models might be more effective. However, how do\nwe interpret the meanings of Random Forest model parameters?", "reviewer_2": "The aims of this study are to disentangle the effects of community composition versus intra-\nspecific variation on leaf N and P concentrations, and to quantify the impacts of\nenvironmental variation. This is important for many reasons and purposes, such as\nprojections of future nutrient cycling in ecosystems and for vegetation model construction,\nspecifically the definition of plastic and species-fixed components of leaf N and P. Given the\ncentral role of nutrient limitation for vegetation dynamics this is an important research\nquestion with wide-ranging potential impact.\nIn contrast to many previous studies the authors use non-linear random forest models (RF) for\nanalysis, which is shown to greatly improve the models’ predictive ability compared to\ntraditional linear models (LM). The analysis convincingly (but diplomatically) debunks the\n“Biogeochemical Niche Hypothesis”, i.e. that inter-specific variation was much more\nimportant than intra-specific. Finally, the authors describe very well the limitations of their\nanalysis, e.g. that it doesn’t separate the effects of individual acclimation and intra-specific\ngenetic variation, which would be nice to do for modeling purposes. The text is very nicely\nwritten, and the figures are very accessible.\nThe only thing that I found slightly strange is the inclusion of modelled stomatal\nconductance and Vcmax as environmental predictors. These are plant physiological\nparameters. But because they are modelled here based on environmental factors, I guess\nthey could be seen as derived aggregated environmental factors. Maybe it would be worth\nmentioning which environmental variables were used for their calculation? A grain of\nuncertainty for me is that I am not at all an expert in the statistical models used, so I cannot\njudge whether the RF model used is the best approach or if even more appropriate\nstatistical methods are available. But as far as I can see, the methods used are very effective\nand quite innovative.\nIn summary, it looks like a brilliant paper.\nOskar Franklin", "reviewer_3": "Review of MS for Nature Comms\nThis ms reports on leaf nutrient contents (N, P N:P) along environmental gradients, asking\nwhat the relative contributions are of environmentally-driven variation within species,\ncompared to species turnover. The question is important because leaf nutrients are\nimportant components of the Leaf Economic Spectrum that reflects a fundamental plant\ntrade-off between acquisitive strategies for carbon assimilation and conservative strategies\nfor long lasting leaf tissue. These traits are key to biomass productivity and are reflective of\nenvironment and changes in it.\nAs with the classic “nature vs nurture” argument, the background here is that opposing\nviewpoints may propose either (1) that species have fixed traits and that variation among\nsites is due to species turnover along environmental gradients. Contrastingly, the hypothesis\nis (2) that species are highly plastic and that observed leaf nutrients reflect environmental\ndrivers within species. Arguments and evidence have been presented for both in the\nliterature. Here the authors attempt to resolve or measure the balance of evidence for\neach—their relative contributions—for leaf N & P. I like the idea, but I am concerned by\nsome of the methodology, and I found that the presentation needs improvement.\nI believe that the ms finds that both potential explanations have some support. Yet I found it\nquite difficult to glean an overall sense. This seems to be because: (a) methodology matters\nto the quantification; and (b) because in different places, the text variously emphasises one\nor other explanation, which may be a matter of expression; and (c) lastly because at least\nsome of the motivation for the ms appears to be to test /falsify a claim that species have\nfixed chemistry, with negligible effect of environment within-species. So, I found myself\nconfused, as other readers might be too.\nTo explain, the Title supports environment (by omitting species turnover). The results are\nreported in the Abstract as supporting an interpretation that Environment is of\ngreat/overriding importance. (there is a single line buried in the middle of the abstract\nstating similar effects, this is swamped by 7 successive lines emphasising environment).\nThen if I look at Fig 3 c,f,I, I have the impression species identity dominates. I hope that this\nindicates the problem.\nThe ms analyses a large, global dataset of leaf nutrients. Over 3,600 species across at >7500\nsites and retaining at least 5 Site measurements. This is a strong data set covering\nclimate space well.\nI found the flow of the manuscript confused me as to the focus. The abstract, in particular,\nshould be rewritten, because it appears that the topic of the paper was methodological\ncomparison of two vastly different analytical methodologies. That is not the case that the\nauthors had in mind I believe. I felt that the analytical comparison was subsidiary to the\nmain question or environmentally driven plasticity vs species turnover. But I may be wrong.\nA considerable part of the manuscript focusses on methodological concerns in the\nattribution of variance explanation between species and environment. Specifically, the\nauthors draw attention to the case of random effects in mixed effects modelling. It is true\nthat the random effects soak up variation, some of which may be explained by fixed effects\n(covariates of interest). So, in a way of speaking, the random effects—species in this case—\nhave priority. To work around this the authors use LM with species as a fixed effect either\nadded before or after environment, to get a sense of the explained variation.\nTrait-Gradient Analysis is an appealing methodology. But it has a flaw. It seeks to regress a\ntrait measured on species i at site j on the community mean trait at site j. That species i\nmust be part of the community at site j and contribute to the community mean trait.\nTherefore, the predictor and response are dependent. This dependence would be greater\nfor more abundant species, for more widespread species and for less diverse communities.\nThis makes the slope estimates unreliable, to an unknown degree.\nThe ms has potential to resolve some of the conflict in the literature about the relative\nimportance of explanations of Leaf N&P. I hope that my comments here help in achieving\nthat potential.\nSpecific comments\nTitle – this suggests (by omission) that variation across species is negligible.\nL39 plasticity should be mentioned.\nL43-44 emphasise this, it gets lost in following 2 sentences. Potentially you could append\nwith “thus challenging claims that species turnover is paramount”\nL44-48 these 2 sentences read as though they are contrasting patterns with respect to your\nresearch question (i.e. within spp environmentally-driven variation or species turnover).\nBut it isn’t, as far as I can tell. Suggest break into 2 sentences and introduce RF result with\nrespect to species identity.\nL48-51 yet this sentence appears to claim predominance of environment, contradicting\nearlier statement that both species and environment matter equally.\nL51 ff. where is this in the results?\nL76-79 this sentence does not seem to match the statement before about drivers.\nL103-105 unclear how this methodological comparison relates to the question, rather than\nbeing a critique of a previous (paper using a) method.\nL120 ff. this section is written without regard for the ecological question at hand. Highlight\nthe ecology, not the stats.\nL136 can you specify that the recursive feature elimination applies in RF models? If indeed it\ndoes.\nL140, and what are brown bars?\nL141 as far as I am aware, t-values are not measure of effect size or magnitude, they are\nmeasures of signal/noise.\nIf predictor variables have all been centred and scaled by a SD, i.e., ‘standardised’, one can\nuse the coefficient as a measure of effect size. Or in a nonlinear model one can simply use\nthe predicted range (i.e. predicted max – predicted min) across a particular covariate.\nL185-188 when you specify the species, if the species have restricted distributions, the\nname encodes something about the environment. it is not truly independent of\nenvironment.\nL204-207 what does this mean? the order of calculation is unclear to me, and thus what it is\nmeant to represent. Is it how much the trait varies among species at a site as a proportion\nof the site mean? Those kernels in Fig 4 e don’t seem accurate, presumably a range of 0 is\npractically impossible. It means that either you have single species at a site, or the trait\nvalues are identical for the species at a site. Or am I missing something here?\nFig 3 c,f,i suggest species identity dominates. Yet Fig 4 f suggests practical equivalence of the\nrelative contributions (40/50/60% across species). Where does Fig 4 f come from?\nFig 4 g,h,i these estimates should have uncertainty bars.\nL238-245 this is really good. It is clear about the questions in conclusions, but I found that\nclarity hard to find elsewhere.\nSupp Fig 1. Where does this result come from?\nL486-> which results specifically Fig 3? Supp Fig 1?\nL515-516 this is the dependency that I mentioned in general comments.\nL535-542 this is unclear to me. Is this the method used to produce Supp Fig 1 or Fig 4f??\nAlso, the actual calculation here is underspecified. What is x_i,k? I think i the species,\nj is the site, but what is k?"}, "author_response": "Authors' response to reviewer’s comments on the manuscript “Widespread\nenvironment-driven plasticity in species’ leaf nitrogen and phosphorus concentrations”\nby Tian et al.\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\n[Comment]: The present study appears to be an important work. However, it is uncertain\nwhether it represents a substantial advancement in the field and if the methods employed\nsufficiently support the arguments made in this paper. Specific concerns are as follows. Please\nconsider addressing these concerns to enhance the clarity and depth of your research findings\nin the paper.\n[Reply]: We thank the reviewer for his or her appreciation of the importance of our work. We\nare convinced that our paper constitutes a substantial advancement in the field for the following\nreasons.\n Leaf N, P, and N:P are key traits for photosynthesis and biogeochemical cycling and our\nunderstanding of large-scale controls are directly relevant for informing global vegetation\nmodels. However, published findings appear to contradict each other. Recent publications\nhave reported an almost exclusive control of leaf N, P, and N:P by phylogeny with a very\nminor role of environmental controls (Sardans et al. 2021; Vallicrosa et al. 2022). Yet,\nenvironmental variables are emphasized as important drivers for global patterns in another\ngroup of studies and used for producing globally upscaled maps of foliar nutrient\nconcentrations (e.g., Han et al. 2005; Reich et al. 2010; Zhang et al. 2012; Maire et al.\n2015; Asner et al. 2016; Tian et al. 2018; Tang et al. 2018; Dechant et al. 2023 Preprint).\n The methodological choices underlying these two groups of studies and implications for\nour understanding of the global ecology of foliar nutrient concentrations have never been\naddressed explicitly and their implicit conflict has not been discussed, let alone resolved in\nprevious literature. Our results demonstrate that methodological choices underpin the\nconflicting conclusions of large-scale controls on leaf N, P, and N:P presented in previously\npublished, high-profile studies (e.g., Sardans et al. 2021; Vallicrosa et al. 2022).\n We thus reconcile an implicit conflict and show a way forward, demonstrating the inherent\nlimitation of linear mixed-effect models (LMM) for separating within vs across species\nvariations in traits and demonstrating the power of more flexible, tree-based methods\n(Random Forest, RF) for this modeling task. We quantify the relative contribution of\nspecies identity and environmental factors in driving variations in the data for leaf N, P,\nand N:P. These estimates contrast with the results by Sardans et al. 2021, published in\nNature EE, and thus suggest that our understanding of dominating controls on these foliar\ntraits must be revised.\n In addition to the RF analysis, we apply traits-gradient analysis for the first time to quantify\nthe plasticity of global leaf N, P, and N:P within species and compare it among\nenvironmental settings and species. These are all novel insights that directly inform the\nfundamental question we address: what is the relative importance of species identity vs\nenvironment in controlling global variation in leaf N, P, and N:P?\nThe methods employed in our study support our arguments very directly in that the types of\nstatistical models (LMM and RF) used in our study reflect the two main methods that have been\nused in the published literature on leaf N, P, and N:P modeling (e.g., Campetella et al. 2011;\nBergmann et al. 2017; Mundim et al. 2021; Fajardo et al. 2022; see our reply in response to\n‘Lines 44-48: What is the biological significance of comparing Random Forest models’ below).\nBy contrasting their results, we are able to directly discuss an implicit conflict that was inherent\nin the published literature (using environment for global spatial modeling vs the claim that the\nenvironment is unimportant for leaf N, P, and N:P). (See also our reply to the reviewer comment\non Lines 44-48).\nSpecific comments:\n[Comment]: Lines 44-48: What is the biological significance of comparing Random Forest\nmodels, trait gradient analysis, and linear models? Do these comparisons reflect specific\nbiological meanings or mechanisms?\n[Reply]: We revised text in several places to clarify the relevance of these methods in the\ncontext of deciphering the power of environment vs. phylogenetic controls on leaf N and P\nconcentrations of terrestrial plants. The revised text in the abstract now reads:\nThe substantial influence of environmental variables on within-species variations is\nidentified using Random Forest models. In contrast, widely used linear-mixed effect\nmodels missed these effects almost completely. Our analysis demonstrates a substantial\ninfluence of the environment in driving plastic responses of leaf N, P, and N:P within\nspecies. This challenges previous reports of a fixed biogeochemical niche and an\noverriding importance of species distributions in shaping global patterns of leaf N and P\nconcentrations. (l. 46-52)\nTwo aspects are relevant for why we compared Random Forest models (RF) and linear (mixed-\neffect) models. First, these are the most widely used model types in the related literature. Linear\nmixed-effect modeling is widely established in the literature on plant traits (Mitchell et al. 2014;\nKnapp et al. 2016). In more recent years, the use of RF has gained traction (Rahman et al. 2021;\nSoltanikazemi et al. 2022) — presumably for their ease of use (minimal data pre-processing\nand hyperparameter tuning required) and strong performance (Cutler et al. 2007; Liu, Y 2014).\nSecond, both model types can deal with data that contain both categorical (e.g., species identity)\nand continuous variables (e.g., environmental variables). In particular, linear mixed-effect\nmodels are a popular choice for dealing with structured data where structures introduce\nvariations that remain unexplained by the fixed factors (e.g., variations ascribed to species\nidentity) (Duursma et al. 2016; West et al. 2022). In that sense, yes, the model types relate to\nbiological aspects of the data and modeling task at hand.\nWe added text in the Introduction to better explain this background (l. 88-105):\nThe conflicting attributions of observed variation in leaf N and P stoichiometry to\nphylogenetic vs environmental variables are related to an inherent methodological\nchallenge. The separation of these variables is usually undermined by their lack of\nindependence. The distribution of plant species is largely driven by the abiotic\nenvironment35,36. Yet, species do occur over a certain range of environmental conditions.\nTo what extent the environment drives phenotypic plasticity or genetic adaptation in leaf\nN and P stoichiometry also within species remains challenging to detect but is\ninformative for testing the Biogeochemical Niche Hypothesis. Linear mixed-effect\nmodels have been widely employed for separating phylogenetic and environmental\neffects on leaf traits24,37, motivated by their suitability to model structured data and their\nability to control for phylogenetic effects and species identity as random terms, implicitly\nassuming that they are unrelated to the environment and given precedence over the latter\nin model fitting. More recently, tree-based statistical learning methods, for example\nRandom Forest, have been shown to be suitable for modeling leaf N and P38,39. These\nmodels, too, provide a natural way to simultaneously account for environmental\n(continuous) and phylogenetic (categorical) information. However, implications of\nmethodological choices for separating environmental vs phylogenetic variables so far\nhave not been explicitly considered.\n[Comment]: Lines 90-92: This statement is difficult to comprehend. What is meant by \"cannot\nbe independent\"? In the natural world, biological organisms and their environments are\ninherently adaptive processes, making complete independence challenging. Does it refer to the\nlack of independence between species identity and historical or current environmental factors?\nLinear mixed effect models can also analyze the interaction between species identity and\nenvironment.\n[Reply]: Following your suggestion, we reformulated text here. As mentioned above (response\nto comment “Lines 44-48 …”), relevant text now reads:\nThe separation of these variables is usually undermined by their lack of independence.\nThe distribution of plant species is largely driven by the abiotic environment35,36. Yet,\nspecies do occur over a certain range of environmental conditions and to what extent\nenvironment-driven plasticity or genetic adaptation drives variations in leaf N and P\nstoichiometry also within species remains challenging to detect but is informative for\ntesting the Biogeochemical Niche Hypothesis. (l. 90-95)\n[Comment]: Line 93: What about the interaction between environmental controls and species\nidentity? It might not be an absolute dichotomy; their interactions could vary with changing\nenvironmental conditions or plant functional groups.\n[Reply]: Thank you for this important point. We agree that there might not be an absolute\ndichotomy between environmental controls and species identity. We also agree that in addition\nto the non-orthogonal (“separation of these variables is usually undermined by their lack of\nindependence”) main effects of phylogeny and environment there could be phylogeny-by-\nenvironment interactions, i.e. different reaction norms of species to environmental variables.\nIndeed, that’s what we find using the trait gradient analysis. We revised text here. It now reads\n(l. 106-111):\nIn view of these conflicting reports and methodological challenges, the question arises\nto what extent large-scale leaf N and P stoichiometric patterns are a reflection of different\nspecies (with their relatively fixed leaf nutrient stoichiometry) occurring at different sites\nalong environmental gradients, and to what extent plasticity and genetic adaption, driven\nby the environment, drive variation within species and contribute to large-scale patterns\nof leaf nutrient stoichiometry (Fig. 1).\n[Comment]: Lines 103-105: Why compare these three methods? Are there no better\nalternatives? Do the differences in these statistical methods hold any biological significance?\n[Reply]: We replied to this point above (see point Lines 44-48 …). As described in the\nIntroduction section, previous contrasting viewpoints regarding the importance of phylogenetic\nvs environmental variables in shaping stoichiometric patterns in leaf N and P concentrations\nand N:P ratios are based on conflicting results from the linear regression model (LM) and linear\nmixed-effect models (LMMs). Compared with traditional LM and LMMs, RF have advantages\nincluding better performance for large datasets, ability to model complex interactions among\npredictor variables, learn non-linear effects of multiple environmental variables more\neffectively, and so on. Hence, we compared the power of the LM, LMMs and RF in explaining\ndifferent components of variations in leaf N and P stoichiometry including variations across\nsites, within species, and across species. Our results highlight the importance of choosing\nsuitable statistical models for revealing within-species variations and for properly reflecting the\ninfluence of environmental variables on interspecific and intraspecific variations of leaf N and\nP stoichiometry. We believe that differences in these statistical methods hold important\nbiological significance in quantifying roles of phylogeny and environment on plant nutrient\nstoichiometry which implies diverse plant species’ adaptation strategies during long-term\nevolution processes and their responses (especially trait plasticity) to current environmental\nchanges. These biological significances are now discussed on lines 359-375.\n[Comment]: Why compare variation within species and between species?\n[Reply]: We compared variation within species and between species to answer the key\nscientific question here. We revised text in the Introduction to explain how the identification of\nvariation within species is informative for our research questions.\nText on l. 81-84 reads:\n[The “Biogeochemical Niche Hypothesis”] posits that each species is characterized by a\nfixed leaf stoichiometry (low within-species variability) [...]\nText on l. 93-95 reads:\nTo what extent the environment drives phenotypic plasticity or genetic adaptation drives\nin leaf N and P stoichiometry also within species remains challenging to detect but is\ninformative for testing the Biogeochemical Niche Hypothesis.\n[Comment]: Does within-species variation change with plant developmental stages, age, or\nindividual size?\n[Reply]: To obtain empirical evidence for a clear effect of developmental stage, age, or size on\nwithin-species variation of leaf N and P concentrations for global terrestrial plants would need\nhuge projects in the future, which is out of the scope of our current dataset. Here, we had to\nimplicitly assume that our data have no strong sampling biases with respect to plant\ndevelopmental stage, age, or size. Because we did not find these variables coded together with\nthe variables used in our analysis, we could unfortunately not test for such biases nor remove\npotential variation caused by plant developmental stage, age, or size (which might have reduced\nresiduals and increased the strength of observed effects). The implicit assumption that original\nstudies sampled plants of standardized developmental stage, age, and size across phylogeny\nand environment are typical for published large-scale analyses of leaf N, P, and N:P (e.g., Han\net al. 2005; Reich et al. 2010; Zhang et al. 2012; Maire et al. 2015; Asner et al. 2016; Tian et\nal. 2018; Sardans et al. 2021; Vallicrosa et al. 2022).\nTo improve the global data set quality to the greatest extent, we included only those records of\npaired N and P concentrations of green and ripe leaves with detailed location information and\nexcluded all records without site information or with unpaired N-P records. For the data\ncollected via field investigation, we confirmed that all plant samples were harvested during the\ngrowing season and the duplicated records were carefully removed. Field sampling was\nprimarily completed from July to August (details reported in Tian et al. 2019). Considering the\nconsolidated sampling methods, large number of records and wide spatial distribution, we\nbelieve that our database here provides a useful and credible tool for the study of leaf nutrient\nstoichiometry.\n[Comment]: Figure 2: Are there issues of collinearity among the predictors? How does\ncollinearity affect the results of the three methods? How many predictors are used in the models?\nIf there are numerous predictors, could the models' abilities to handle complex data differ,\npotentially affecting the results?\n[Reply]: We designed our analysis with a strong view to mitigating the influence of collinearity.\nThis is extensively reported (see Results, first subsection ‘Variable selection and effects’ and\nFig. 2). The purpose of Fig. 2 is indeed to show that a relatively small subset of predictors is\nsufficient for modeling and that additional predictors do not add information — pointing to\ncollinearity among them. A computationally intensive recursive feature elimination with a\ncross-validation at each step was performed. The final selection of seven predictors for leaf N\nand nine predictors for leaf P and N:P corresponds to a complementary set of environmental\nvariables, including climatic and edaphic, CO and N-deposition.\nThe identified important\n2\npredictors are briefly discussed in the context of reports in the literature on l. 384-392:\nGlobal mean atmospheric CO was among the most important selected variables for all\n2\nthree target variables and the single most important predictor for leaf N:P. The direction\nof the effect on leaf N identified here (decline with increasing CO ) is consistent with\n2\nobservations from Free Air CO Enrichment (FACE) experiments71,74,80,81 and with the\n2\nmechanistic understanding of the influence of CO on N demand82.\nAmong the strongest effects identified here was a positive effect of N-deposition on leaf\nN (consistent with refs. 52-56) and N:P (consistent with refs. 10,57-59), a negative effect\nof the temperature of the coldest month on leaf P, and a positive effect of the same on\nleaf N:P (consistent with refs. 26,30,60). These apparently robust trait–environment\nrelationships have given rise to spatial upscaling25,61, mapping leaf traits with global\ncoverage.\nWe are therefore confident that our approach effectively mitigates collinearity and associated\nrisks of model overfitting, and that the selected factors represent a set of robustly identified\ndrivers of leaf N, P, and N:P variations.\n[Comment]: Lines 241-242: How can the comparability of variation within and between\nspecies be ensured?\n[Reply]: In quantitative terms, the separation of variation within and across species is based on\ncomparing the variance of species means with the pooled variance within species as described\non l. 574-581 (“species variation decomposition”). The physical units (g N (P) g DM-1)2 of the\ntwo variances within and across species are therefore identical and thus the coefficients of\ndetermination calculated from the two types of modifications (the first replacing individual\nvalues with the corresponding species means and the second replacing individual values with\nthe deviations from their species mean) can be compared. In terms of ecological meaning,\ncomparing variations across vs. within species is interpretable with respect to fundamental\nquestions related to the leaf economics spectrum, trait–environment relationships, or the role\nof phylogeny and other plant classifications — as demonstrated by a rich literature (e.g.,\nAnderegg et al. 2018; Umaña et al. 2019; Fyllas et al. 2020). We added a cross-reference to the\nMethods section in the caption of Fig. 4 where results from the species variation decomposition\nare displayed (panel f).\n[Comment]: Lines 249-250: What distinguishes \"variation\" from \"plasticity\"? What is the\nsignificance of within-species variations?\n[Reply]: The term plasticity, or phenotypic plasticity, is commonly used to refer to the\nvariability of a trait in response to environmental stimuli or the capacity of an organism to alter\nits phenotype in response to varying environmental conditions (Valladares et al. 2014; Adu et\nal. 2022). It is a biological phenomenon. In contrast, variation relates purely to the data.\nThe relationship between within-species variations and environmental variables are interpreted here\nas an indicator of phenotypic plasticity or genetic adaptation of populations within a species\n(Stotz et al. 2021).\nTo better explain the significance of within-species variation for our research question, we\nmodified text in several places.\nOn l. 81-84:\n[The “Biogeochemical Niche Hypothesis”] posits that each species is characterized by a\nfixed leaf stoichiometry (low within-species variability), [...]\nOn l. 93-95:\nTo what extent the environment drives phenotypic plasticity or genetic adaptation in leaf\nN and P stoichiometry also within species remains challenging to detect but is\ninformative for testing the Biogeochemical Niche Hypothesis.\n[Comment]: Lines 305-309: Linear mixed effect models can account for nonlinear\nrelationships and interactions, although Random Forest models might be more effective.\nHowever, how do we interpret the meanings of Random Forest model parameters?\n[Reply]: We thank the reviewer for pointing out that mixed-effect models may indeed be\nformulated to deal with non-linear and interactive responses. We revised text to be more precise\nin our critique, which pertains specifically to how LMMs were used in the published literature\non leaf N, P, and N:P (l. 310-317):\nThe large shared effect between environment and species identity (Supplementary Fig.\n1) cannot be decomposed without relying on targeted experimental designs and their\ninteractions and non-linearities are not considered in published LMM-based analyses20,24.\nOur LMM model specification imitated their methodological choices. By design, in their\nLMMs the shared effect gets attributed to species identity if used as random-effects term,\nrather than to environmental variables used as fixed-effects terms in LMMs.\nWe added a sentence on l. 325-327 pointing to the limitation regarding interpretability\nof fitted model parameters in Random Forest analysis.\nHowever, fitted model coefficients of RF are not always directly interpretable — in\ncontrast to coefficients of LMMs.\nFurther research will be useful for diagnosing functional relationships from RF models fitted\nusing machine learning. However, we consider this to be beyond the scope for the present paper\nwhich focuses on identifying effects based on their importance for out-of-sample predictions.\nReferences cited above:\nAdu, M. O. et al. Root system traits contribute to variability and plasticity in response to\nphosphorus fertilization in 2 field-grown Sorghum [Sorghum bicolor (L.) Moench]\ncultivars. Plant Phenomics 2022, 0002 (2022).\nAnderegg, L. D. et al. Within‐species patterns challenge our understanding of the leaf\neconomics spectrum. Ecol. Lett. 21, 734-744 (2018).\nAsner, G. P., Knapp, D. E., Anderson, C. B., Martin, R. E. & Vaughn, N. Large-scale climatic\nand geophysical controls on the leaf economics spectrum. PNAS 113, E4043-E4051\n(2016).\nBergmann, J., Ryo, M., Prati, D., Hempel, S. & Rillig, M. C. Root traits are more than analogues\nof leaf traits: the case for diaspore mass. New Phytol. 216, 1130-1139 (2017).\nCampetella, G. et al. Patterns of plant trait-environment relationships along a forest succession\nchronosequence. Agr. Ecosyst. Environ. 145, 38-48 (2011).\nCutler, D. R. et al. Random forests for classification in ecology. Ecology 88, 2783-2792 (2007).\nDechant, B. et al. Intercomparison of global foliar trait maps reveals fundamental differences\nand limitations of upscaling approaches. EarthArXiv 1-57 (2023).\nDuursma, R. & Powell, J. Mixed-effects models. Sydney: Hawkesbury Institute for the\nEnvironment (2016).\nFajardo, A., Piper, F. I. & García‐Cervigón, A. I. The intraspecific relationship between wood\ndensity, vessel diameter and other traits across environmental gradients. Funct. Ecol. 36,\n1585-1598 (2022).\nFyllas, N. M. et al. Functional trait variation among and within species and plant functional\ntypes in mountainous mediterranean forests. Front. Plant Sci. 11, 212 (2020).\nHan, W., Fang, J., Guo, D. & Zhang, Y. Leaf nitrogen and phosphorus stoichiometry across\n753 terrestrial plant species in China. New Phytol. 168, 377-385 (2005).\nKnapp, S., Stadler, J., Harpke, A. & Klotz, S. Dispersal traits as indicators of vegetation\ndynamics in long-term old-field succession. Ecol. Indic. 65, 44-54 (2016).\nLiu, Y. Random forest algorithm in big data environment. Computer modelling & new\ntechnologies 18, 147-151 (2014).\nMaire, V. et al. Global effects of soil and climate on leaf photosynthetic traits and rates. Glob.\nEcol. Biogeogr. 24, 706-717 (2015).\nMitchell, R. M., & Bakker, J. D. Quantifying and comparing intraspecific functional trait\nvariability: a case study with Hypochaeris radicata. Funct. Ecol. 28, 258-269 (2014).\nMundim, F. M., Vieira‐Neto, E. H., Alborn, H., & Bruna, E. M. Disentangling the influence of\nwater limitation and simultaneous above and belowground herbivory on plant tolerance\nand resistance to stress. J. Ecol. 109, 2729-2739 (2021).\nRahman, M. et al. Disentangling the role of competition, light interception, and functional traits\nin tree growth rate variation in South Asian tropical moist forests. Forest Ecol. Manag.\n483, 118908 (2021).\nReich, P. B. et al. Evidence of a general 2/3-power law of scaling leaf nitrogen to phosphorus\namong major plant groups and biomes. Proc. Natl Acad. Sci. USA 277, 877-883 (2010).\nSardans, J. et al. Empirical support for the biogeochemical niche hypothesis in forest trees. Nat.\nEcol. Evol. 5, 184-194 (2021).\nSoltanikazemi, M., Minaei, S., Shafizadeh-Moghadam, H. & Mahdavian, A. Field-scale\nestimation of sugarcane leaf nitrogen content using vegetation indices and spectral bands\nof Sentinel-2: Application of random forest and support vector regression. Comput\nElectron Agr. 200, 107130 (2022).\nStotz, G. C., Salgado-Luarte, C., Escobedo, V. M., Valladares, F. & Gianoli, E. Global trends\nin phenotypic plasticity of plants. Ecol. Lett. 24, 2267-2281 (2021)\nTang, Z. et al. Patterns of plant carbon, nitrogen, and phosphorus concentration in relation to\nproductivity in China's terrestrial ecosystems. Proc. Natl. Acad. Sci. USA 115, 4033-4038\n(2018).\nTian, D. et al. A global database of paired leaf nitrogen and phosphorus concentrations of\nterrestrial plants. Ecology 9, e02812 (2019).\nTian, D. et al. Global leaf nitrogen and phosphorus stoichiometry and their scaling exponent.\nNatl. Sci. Rev. 5, 738-739 (2018).\nUmaña, M. N. & Swenson, N. G. Does trait variation within broadly distributed species mirror\npatterns across species? A case study in Puerto Rico. Ecology. 100, e02745 (2019).\nValladares, F. et al. The effects of phenotypic plasticity and local adaptation on forecasts of\nspecies range shifts under climate change. Ecol. Lett. 17, 1351-1364 (2014).\nVallicrosa, H. et al. Global maps and factors driving forest foliar elemental composition: the\nimportance of evolutionary history. New Phytol. 233, 169-181 (2022).\nWest, B. T., Welch, K. B. & Galecki, A. T. Linear mixed models: a practical guide using\nstatistical software. Crc Press (2022).\nZhang, S. B., Zhang, J. L., Slik, J. & Cao, K. F. Leaf element concentrations of terrestrial plants\nacross China are influenced by taxonomy and the environment. Glob. Ecol. Biogeogr. 21,\n809-818 (2012).\nReviewer #2 (Remarks to the Author):\n[Comment]: The aims of this study are to disentangle the effects of community composition\nversus intra-specific variation on leaf N and P concentrations, and to quantify the impacts of\nenvironmental variation. This is important for many reasons and purposes, such as projections\nof future nutrient cycling in ecosystems and for vegetation model construction, specifically the\ndefinition of plastic and species-fixed components of leaf N and P. Given the central role of\nnutrient limitation for vegetation dynamics this is an important research question with wide-\nranging potential impact.\nIn contrast to many previous studies the authors use non-linear random forest models (RF) for\nanalysis, which is shown to greatly improve the models’ predictive ability compared to\ntraditional linear models (LM). The analysis convincingly (but diplomatically) debunks the\n“Biogeochemical Niche Hypothesis”, i.e. that inter-specific variation was much more important\nthan intra-specific. Finally, the authors describe very well the limitations of their analysis, e.g.\nthat it doesn’t separate the effects of individual acclimation and intra-specific genetic variation,\nwhich would be nice to do for modeling purposes. The text is very nicely written, and the figures\nare very accessible.\nThe only thing that I found slightly strange is the inclusion of modelled stomatal\nconductance and Vcmax as environmental predictors. These are plant physiological\nparameters. But because they are modelled here based on environmental factors, I guess\nthey could be seen as derived aggregated environmental factors. Maybe it would be worth\nmentioning which environmental variables were used for their calculation? A grain of\nuncertainty for me is that I am not at all an expert in the statistical models used, so I cannot\njudge whether the RF model used is the best approach or if even more appropriate\nstatistical methods are available. But as far as I can see, the methods used are very effective\nand quite innovative.\nIn summary, it looks like a brilliant paper.\nOskar Franklin\n[Reply]: We thank the reviewer for the appreciation of the value of our research.\nRegarding the derived environmental variables (V , J , g), we added justification for the\ncmax25 max25 s\nuse of these as reflections of climate input variables. The complemented text in the Methods\nsection (l. 459-463) reads:\nDue to the intrinsic97 and widely observed98 link between V , Rubisco, and leaf N\ncmax25\nand P contents, we used these estimates here as predictors for leaf N, P, and N:P,\nthus\nattempting to account for simultaneous effects of multiple climate variables (air\ntemperature, VPD, elevation, irradiance, CO ) on leaf nutrient stoichiometry.\nclimatic information relevant for modeling leaf N, P, and N:P. However, the feature\nelimination did not retain these variables and they were not used for all subsequent\nanalyses presented in Figs. 3 and 4. Therefore, we did not include a more extensive\ndiscussion of this result in the paper (also given the space constraints for the text).\nRegarding the choice of RF as a statistical modelling tool, we now provide more information\nfor justifying it (also in response to Reviewer #1). Text on l. 88-105 reads:\nThe conflicting attributions of observed variation in leaf N and P stoichiometry to\nphylogenetic vs environmental variables are related to an inherent methodological\nchallenge. The separation of these variables is usually undermined by their lack of\nindependence. The distribution of plant species is largely driven by the abiotic\nenvironment35,36. Yet, species do occur over a certain range of environmental conditions.\nTo what extent the environment drives phenotypic plasticity or genetic adaptation in leaf\nN and P stoichiometry also within species remains challenging to detect but is\ninformative for testing the Biogeochemical Niche Hypothesis. Linear mixed-effect\nmodels have been widely employed for separating phylogenetic and environmental\neffects on leaf traits24,37, motivated by their suitability to model structured data and their\nability to control for phylogenetic effects and species identity as random terms, implicitly\nassuming that they are unrelated to the environment and given precedence over the latter\nin model fitting. More recently, tree-based statistical learning methods, for example\nRandom Forest methods, have been shown to be suitable for modelling leaf N and P 38,39.\nThese models, too, provide a natural way to simultaneously account for environmental\n(continuous) and phylogenetic (categorical) information. However, implications of\nmethodological choices for separating environmental vs phylogenetic variables so far\nhave not been explicitly considered.\nReviewer #3 (Remarks to the Author):\nReview of MS for Nature Comms\n[Comment]: This ms reports on leaf nutrient contents (N, P N:P) along environmental gradients,\nasking what the relative contributions are of environmentally-driven variation within species,\ncompared to species turnover. The question is important because leaf nutrients are important\ncomponents of the Leaf Economic Spectrum that reflects a fundamental plant trade-off between\nacquisitive strategies for carbon assimilation and conservative strategies for long lasting leaf\ntissue. These traits are key to biomass productivity and are reflective of the environment and\nchanges in it.\nAs with the classic “nature vs nurture” argument, the background here is that opposing\nviewpoints may propose either (1) that species have fixed traits and that variation among sites\nis due to species turnover along environmental gradients. Contrastingly, the hypothesis is (2)\nthat species are highly plastic and that observed leaf nutrients reflect environmental drivers\nwithin species. Arguments and evidence have been presented for both in the literature. Here the\nauthors attempt to resolve or measure the balance of evidence for each—their relative\ncontributions—for leaf N & P. I like the idea, but I am concerned by some of the methodology,\nand I found that the presentation needs improvement.\n[Reply]: We appreciate the very valuable feedback that helped us to improve the manuscript.\nIndeed, our research can be framed as an attempt to balance evidence and quantify relative\ncontributions for contrasting interpretations. We have now adopted this “spot-on” formulation\nin the abstract (l. 39-40):\nHere, we balance these contrasting views via measuring their relative contributions to\nthe global leaf N/P stoichiometric pattern.\n[Comment]: I believe that the ms finds that both potential explanations have some support. Yet\nI found it quite difficult to glean an overall sense. This seems to be because: (a) methodology\nmatters to the quantification; and (b) because in different places, the text variously emphasises\none or other explanation, which may be a matter of expression; and (c) lastly because at least\nsome of the motivation for the ms appears to be to test /falsify a claim that species have fixed\nchemistry, with negligible effect of environment within-species. So, I found myself confused,\nas other readers might be too.\nTo explain, the Title supports environment (by omitting species turnover). The results are\nreported in the Abstract as supporting an interpretation that Environment is of\ngreat/overriding importance. (there is a single line buried in the middle of the abstract\nstating similar effects, this is swamped by 7 successive lines emphasising environment). Then if I look at Fig 3 c,f,i, I\nhave the impression species identity dominates. I hope that this indicates the problem.\n[Reply]: Indeed, our framing was initially guided strongly by the testing of the Biogeochemical\nNiche Hypothesis. To more accurately reflect our finding that within- and across-species\nvariations contribute about equally to overall variation in the data and the fact that we quantify\ntheir relative contributions and the influence of the environment on each, we revised text in\nseveral places. First, we changed the title to:\nEnvironmental versus phylogenetic controls on leaf nitrogen and phosphorous\nconcentrations of land plants\nSecond, we revised the abstract (l. 36-52):\nGlobal patterns of leaf nitrogen (N) and phosphorus (P) stoichiometry have been variably\ninterpreted as reflecting phenotypic plasticity in response to the environment, or an\noverriding effect of the distribution of species growing in their biogeochemical niche.\nHere, we balance these contrasting views via measuring their relative contributions to\nthe global leaf N/P stoichiometric pattern. We completed a comprehensive global data\nset and investigate how species identity and environmental variables control variation in\nleaf N, P, and N:P. We find that within-species variation contributes around half to the\ntotal variation in these three leaf variables, 29%, 31%, and 22% of which, respectively,\nare explained by environmental variables. Within-species plasticity along environmental\ngradients is highest for leaf N:P and lowest for leaf N and varies across species. The\nsubstantial influence of environmental factors on within-species variations was identified\nby Random Forest models. In contrast, widely used linear mixed-effect models missed\nthese effect almost completely. Our analysis demonstrates a substantial influence of the\nenvironment in driving plastic responses of leaf N, P, and N:P within species. This\nchallenges previous reports of a fixed biogeochemical niche and an overriding\nimportance of species distributions in shaping global leaf patterns of N and P\nconcentrations.\nThird, we complemented text in the Discussion to better explain that the environmental\nvariables influence both species distribution and leaf N, P, and N:P variations within species.\nThe text on l. 257-266 now reads:\nThe influence of the environment as a driver of leaf N, P, and N:P is demonstrated by\nthree results presented here. First, environmental variables explain 30–45% of variation\nin community-weighted means across sites (Fig. 3a, d, g). This reflects the environmental\nfiltering of species occurrence across environmental gradients. Second, besides species\ndistribution, the environment influences leaf N, P, and N:P directly, driving variation\nwithin species. Random Forest models explained around 20–30% of this variation, as\nshown in Figure 3b, e, and h. Third, a large proportion of species expressed strong\nplasticity in leaf N, P, and N:P, whereby within-species variation paralleled across-site\nvariation. This was expressed by the slope of species-specific regression lines in the trait\ngradient analysis (Fig. 4d).\n[Comment]: The ms analyses a large, global dataset of leaf nutrients. Over 3,600 species across\nat >7500 sites and retaining at least 5 Site measurements. This is a strong data set covering\nclimate space well.\nI found the flow of the manuscript confused me as to the focus. The abstract, in particular,\nshould be rewritten, because it appears that the topic of the paper was methodological\ncomparison of two vastly different analytical methodologies. That is not the case that the\nauthors had in mind I believe. I felt that the analytical comparison was subsidiary to the\nmain question of environmentally driven plasticity vs species turnover. But I may be wrong.\n[Reply]: The scientific question (controls on leaf N and P stoichiometry) and methodological\nchoices are intertwined. We now revised text across the whole manuscript to put the science\nquestion in the foreground and to provide a better explanation of how finding answers to the\nscience question is linked to methodological choices. Mentioning (and explaining) the link to\nthe methods is necessary for discussing our results in the context of the literature. With this in\nmind, we re-wrote the abstract and added a paragraph in the Introduction that explains the link\nto methodological choices (l. 88-90):\nThe conflicting attributions of observed variation in leaf N and P stoichiometry to\nphylogenetic vs environmental variables are related to an inherent methodological\nchallenge. [...] \n[Comment]: A considerable part of the manuscript focuses on methodological concerns in the\nattribution of variance explanation between species and environment. Specifically, the authors\ndraw attention to the case of random effects in mixed effects modeling. It is true that the random\neffects soak up variation, some of which may be explained by fixed effects (covariates of\ninterest). So, in a way of speaking, the random effects—species in this case—have priority. To\nwork around this the authors use LM with species as a fixed effect either added before or after\nenvironment, to get a sense of the explained variation.\nTrait-Gradient Analysis is an appealing methodology. But it has a flaw. It seeks to regress a\ntrait measured on species i at site j on the community mean trait at site j. That species i must be\npart of the community at site j and contribute to the community mean trait. Therefore, the\npredictor and response are dependent. This dependence would be greater for more abundant\nspecies, for more widespread species and for less diverse communities. This makes the slope\nestimates unreliable, to an unknown degree.\n[Reply]: We don’t consider the trait gradient analysis (TGA) as a statistical prediction method\nwhere the predictor and target variables, i.e., the variables plotted along the x and y-axis, must\nbe independent. Across all observations, a correlation between individual observations and the\nsite means emerges by design. The value of the TGA is in its information about the slope of\nspecies-specific regressions. A slope of zero would be indicative of a species-specific trait that\ndoes not vary within species. If indeed, this was an accurate model of the data-generating\nprocess, the TGA would identify it, as demonstrated below and in the supplementary material,\nuploaded along with our re-submission. A trait gradient analysis with synthetic data, generated\nby alternative processes that reflect the Biogeochemical Niche Hypothesis on the one hand and\nperfect stoichiometric plasticity on the other hand is shown in Fig. R1. A potential solution to\nthe problem mentioned by the reviewer would be the use of jackknifing, i.e. leave the species\nbeing tested out of the calculation of the site mean. However, then each species would have a\ndifferent x-axis for the calculation of its slope, making comparisons among slopes biased in a\ndifferent way. We did not use such an approach because TGA in the literature has been\nconceived and applied in the way we used it here (e.g., Ackerly & Cornwell 2007; Kooyman et\nal. 2010; Gallagher & Leishman 2012; Cheng et al. 2016; Ottaviani et al. 2018; Dong et al.\n2020; Zhou et al. 2022) and it would require a more thorough analysis to evaluate jackknifing\nin TGA.\nFig. R1: Trait gradient hypotheses demonstrated by applying the trait gradient method to synthetic data,\ngenerated with a data-generating process that followed the Biogeochemical Niche Hypothesis (left\ncolumn) or the perfect stoichiometric plasticity assumption (right column).\nThe reviewer suggests that the “dependence would be greater for more abundant species, for\nmore widespread species and for less diverse communities.” We tested whether the slopes show\na dependence on abundance (number of observations per species in the data used for the TGA)\nand spread (number of sites at which a given species is observed). Slopes are defined at the\nspecies level. It is unclear how a dependence on community diversity could be tested. We found\nno relationship between the slope and species abundance or spread (see Fig. R2).\nFig. R2: Relationships between slopes in the trait gradient analysis and species abundance (top row) and\nspread (bottom row).\n[Comment]: The ms has potential to resolve some of the conflict in the literature about the\nrelative importance of explanations of Leaf N&P. I hope that my comments here help in\nachieving that potential.\n[Reply]: Indeed, the comments by reviewer 3 were greatly helpful to improve our MS.\nThank\nyou very much.\nSpecific comments\n[Comment]: Title – this suggests (by omission) that variation across species is negligible.\n[Reply]: Thanks for your suggestion. We changed the title to:\nEnvironmental versus phylogenetic controls on leaf nitrogen and phosphorous\nconcentrations of land plants\n[Comment]: L39 plasticity should be mentioned.\n[Reply]: We revised this sentence. It now reads:\nGlobal patterns of leaf nitrogen (N) and phosphorus (P) stoichiometry have been variably\ninterpreted as reflecting phenotypic plasticity in response to the environment, or of an\noverriding effect by the global distribution of species growing in their biogeochemical\nniches [...] (l.36-39)\n[Comment]: L43-44 emphasise this, it gets lost in following 2 sentences. Potentially you could\nappend with “thus challenging claims that species turnover is paramount”\n[Reply]: Thanks. We rewrote the Abstract completely as shown above.\n[Comment]: L44-48 these 2 sentences read as though they are contrasting patterns with respect\nto your research question (i.e. within spp environmentally-driven variation or species turnover).\nBut it isn’t, as far as I can tell. Suggest break into 2 sentences and introduce RF result with\nrespect to species identity.\n[Reply]: Thanks. We rewrote the Abstract completely as shown above.\n[Comment]: L48-51 yet this sentence appears to claim predominance of environment,\ncontradicting earlier statement that both species and environment matter equally.\n[Reply]: We rewrote the Abstract completely — also with a view to better reflect that influences\nof species identity and environmental variables contribute about equally.\n[Comment]: L51 ff. where is this in the results?\n[Reply]: We removed this statement in the revised Abstract.\n[Comment]: L76-79 this sentence does not seem to match the statement before about drivers.\n[Reply]: We removed this sentence and now write:\nOverall, these contrasting viewpoints provide contentious interpretations for global leaf\nN and P stoichiometry variations and take conflicting viewpoints regarding the\nimportance of phylogenetic vs environmental controls in shaping these stoichiometric\npatterns.\n[Comment]: L103-105 unclear how this methodological comparison relates to the question,\nrather than being a critique of a previous (paper using a) method.\n[Reply]: We revised and added text in the Introduction to better explain the relevance of\nmethodological choices in the context of the science question at hand (see also response to\nReviewer 1). Most importantly, we added a paragraph in the Introduction to explain this (l. 88-\n105):\nThe conflicting attributions of observed variation in leaf N and P stoichiometry to\nphylogenetic vs environmental variables are related to an inherent methodological\nchallenge. The separation of these variables is usually undermined by their lack of\nindependence. The distribution of plant species is largely driven by the abiotic\nenvironment35,36. Yet, species do occur over a certain range of environmental conditions.\nTo what extent the environment drives phenotypic plasticity or genetic adaptation in leaf\nN and P stoichiometry also within species remains challenging to detect but is\ninformative for testing the Biogeochemical Niche Hypothesis. Linear mixed-effect\nmodels have been widely employed for separating phylogenetic and environmental\neffects on leaf traits24,37, motivated by their suitability to model structured data and their\nability to control for phylogenetic effects and species identity as random terms, implicitly\nassuming that they are unrelated to the environment and given precedence over the latter\nin model fitting. More recently, tree-based statistical learning methods, for example\nRandom Forest methods, have been shown to be suitable for modelling leaf N and P38,39.\nThese models, too, provide a natural way to simultaneously account for environmental\n(continuous) and phylogenetic (categorical) information. However, implications of\nmethodological choices for separating environmental vs phylogenetic variables so far\nhave not been explicitly considered.\n[Comment]: L120 ff. this section is written without regard for the ecological question at hand.\nHighlight the ecology, not the stats.\n[Reply]: We revised the text to better motivate the variable selection with regards to the\necological question at hand (l. 134-147):\nWe started by identifying the most important environmental variables separately for\nexplaining variations in leaf N, P, and N:P. Reduced predictor sets, specific for leaf N,\nP, and N:P, respectively, enabled an improved model performance compared with\nmodels that included all 45 predictors (Fig. 2a-c) and were used for all subsequent\nanalyses. In LMM models (Fig. 2d-f), N deposition (ndep) had the strongest effect on\nleaf N and leaf N:P variations within species (both positive). The temperature of the\ncoldest month (tmonthmin) had the strongest effect on leaf P (negative). Atmospheric\nCO2 concentrations (co2), ndep, and tmonthmin were among the most important\npredictors for leaf N, P, and N:P. Soil variables were only selected among the most\nimportant variables for leaf N (aluminum saturation of the soil solution, ALSA) and for\nleaf P (soil texture, measured by the water holding capacity class, AWC_CLASS).\n[Comment]: L136 can you specify that the recursive feature elimination applies in RF models?\nIf indeed it does.\n[Reply]: Yes, the recursive feature elimination was done based on RF. We now provide this\ninformation in the caption of Figure 2.:\n(a-c) Variable selection order determined by recursive feature elimination based on\nRandom Forest and 5-fold cross-validation.\n[Comment]: L140, and what are brown bars?\n[Reply]: The brown bars in Figs. 2 a-c indicate the R2of a model with one additional predictor,\nas indicated by the label along the y-axis, added. As there is no clear gain when including\nadditional predictors colored in brown, the final selection of variables is indicated by the green\nbars. In other words, the brown bars indicate the alternative environmental predictors that will\nnot increase the R2of a model. We added text in the caption:\nBrown bars indicate additional, next most important predictors which were not used for\nsubsequent analyses.\n[Comment]: L141 as far as I am aware, t-values are not measures of effect size or magnitude,\nthey are measures of signal/noise. If predictor variables have all been centered and scaled by a\nSD, i.e., ‘standardised’, one can use the coefficient as a measure of effect size. Or in a nonlinear\nmodel one can simply use the predicted range (i.e. predicted max – predicted min) across a\nparticular covariate.\n[Reply]: We thank the reviewer for this insight. Although commonly used as measures for the\nvariable importance, t-values cannot be interpreted as “effect magnitudes”. As suggested by the\nreviewer, we corrected this part of the analysis and describe it accordingly in the caption:\nEffect magnitudes of the selected variables, measured by the coefficients of normalized\nfixed effects in LMMs.\n[Comment]: L185-188 when you specify the species, if the species have restricted distributions,\nthe name encodes something about the environment. it is not truly independent of environment.\n[Reply]: We agree with you that plant species distribution encodes an influence of the\nenvironment and that they are not truly independent of the environment. We added text in the\nIntroduction to clarify this upfront. On l. 88-92, we write:\nThe conflicting attributions of observed variation in leaf N and P stoichiometry to\nphylogenetic vs environmental variables are related to an inherent methodological\nchallenge. The separation of these variables is undermined by their lack of independence.\nThe distribution of plant species is largely driven by the abiotic environment35,36.\n[Comment]: L204-207 what does this mean? the order of calculation is unclear to me, and thus\nwhat it is meant to represent. Is it how much the trait varies among species at a site as a\nproportion of the site mean? Those kernels in Fig 4 e don’t seem accurate, presumably a range\nof 0 is practically impossible. It means that either you have single species at a site, or the trait\nvalues are identical for the species at a site. Or am I missing something here?\n[Reply]: In this paragraph, we report results on the extent of the distribution of species analyzed\nin the trait gradient analysis. The quantification of the species distribution range is defined here\nas the range of values of site-mean leaf N, P, and N:P at which a given species occurs. Of course,\nthe species itself contributes to the site-mean. However, note that …\n“For the TGA, we considered only data from species that were recorded in at least five sites\nand used only data from sites where at least five different species were sampled.” (l.568-569).\nThe order of calculation is described in the Methods section (l. 563-573):\nWe defined ranges as the difference between 1% and 99% quantiles to reduce the influence of\noutliers. Species-specific range values were normalized by the overall range of site-level means\nx to make ranges comparable across the three different traits investigated here. Normalized\nj\nranges are interpreted here as an indicator of the range of site conditions under which the\nrespective species occurs. […] Normalized ranges were quantified based on the original, not\nlog-transformed data.\nThe kernel density plots are correct. They indicate a non-zero density around a range of zero.\nA clearer picture is provided by a histogram, indicating that there are observations in the lowest\nbin, but values themselves are not zero.\nN N:P P\n3\n2\n1\n0\n0.0 0.5 1.0 0.0 0.5 1.0 0.0 0.5 1.0\nNormalised range\nFig. R3: Histogram of normalized ranges. The kernel density lines as shown in Fig. 4 of the manuscript\nare shown on top of histogram bars.\nWe chose to show kernel density plots instead of histograms as they enable data for all three\ntraits to be plotted in the same panel, while enabling a clear visual distinction. Plotting\nhistograms for the three traits on top of each other obscures visual clarity.\n[Comment]: Fig 3 c,f,i suggest species identity dominates. Yet Fig 4 f suggests practical\nequivalence of the relative contributions (40/50/60% across species). Where does Fig 4 f come\nfrom?\n[Reply]: Fig. 4f is titled as Species variation decomposition. We added a reference to Methods\nin the caption of Fig. 4. The respective section in the Method is on l. 574-581. Fig. 4f shows\nytisneD N\nP\nN:P\nhow much of the variation arises within vs across species — irrespective of whether any of\nthese variations are associated with environmental predictors and independent of any modeling.\nIn contrast, Fig. 3 (c, f, i) shows how much of the variation environmental predictors explain\nand is based on modeling. We tried to make this distinction clear. For example, in the Abstract,\nwe write:\nWe find that within-species variations contribute around half to the total variation in\nthese three leaf variables, 29%, 31%, and 22% of which, respectively, are explained by\nenvironmental variables.\n[Comment]: Fig 4 g,h,i these estimates should have uncertainty bars.\n[Reply]: We added the 95% confidence intervals of the slope estimate as grey lines to the bar\nplots in the revision.\n[Comment]: L238-245 this is really good. It is clear about the questions in conclusions, but I\nfound that clarity hard to find elsewhere.\n[Reply]: Thank you. We tried to clarify that “Variation within species between sites is similarly\nlarge as variation among species” also in other places of the manuscript.\n[Comment]: Supp Fig 1. Where does this result come from?\n[Reply]: We added a cross-reference to clarify this (green text below). In the Methods section,\nexisting text reads:\nAdditionally, we separated individual and shared effects of site and species based on\nlinear models where the respective predictors were fitted in different orders (Supp. Fig.\n1).\nSupp. Fig. 1 is referenced on l. 546.\n[Comment]: L486-> which results specifically Fig 3? Supp Fig 1?\n[Reply]: We are not sure what the reviewer refers to. Existing text on l. 486 reads:\n“…variables alone by taking the difference between RF models fitted with\nenvironmental…”.\nIf we understand correctly, the reviewer asked the relationship between our results\nshowed in Fig. 3 and Supp. Fig. 1 in the manuscript. We rewrote this part ‘Contrasting\nmodel performances’ to clarify our result in the revision. The results showed specifically\nin Fig 3 locates in line 165-178 which reads:\nThe Random Forest (RF) models, fitted to site-level aggregated data (mean across all\nobservations by site), with the selected subset of environmental variables as predictors,\nachieved an R2 of 0.46, 0.34, and 0.34 for leaf N, P, and N:P, respectively, in contrast to\nan R2 of 0.17, 0.19, and 0.19 in LMs (Fig. 3 a,d,g). R2 , measuring the proportion of\nmarg\nvariation explained by fixed (environmental) variables in LMMs, fitted to the full data,\nwas only 0.04, 0.05, and 0.09 for leaf N, P, and N:P, respectively. In contrast, the\nproportions of variation explained by species identity were 0.68, 0.63, and 0.44 (ICC),\nrespectively, for leaf N, P, and N:P ratio (Fig. 3 c,f,i). When RF models were fitted to\nthe full data, environmental variables explained a larger proportion of the variation than\nthey did in LMMs, namely 0.13, 0.26, and 0.16 vs 0.04, 0.05 and 0.09, respectively. A\nsimilar contrast in the predictive power of environmental variables in RF and LMs is\nseen with models fitted to the modified data that contained only within-species variation\n(Fig. 3 b,e,h). Here, RF models achieved an R2 of 0.29, 0.31, and 0.22 for leaf N, P, and\nN:P, while LMs achieved an R2 of 0.01, 0.03, and 0.07, respectively.\nThe results showed specifically in Supp. Fig. 1 locates in line 185-187 which reads:\nBased on LMs, the shared effect of species identity and sites was dominant for leaf N, P,\nand N:P, explaining more than double of the variance explained by their separate effects\n(Supplementary Fig. 1).\n[Comment]: L515-516 this is the dependency that I mentioned in general comments.\n[Reply]: We replied to the point about the trait gradient analysis above and wrote that we found\nno relationship between the slope and species abundance or spread (see Fig. R2).\n[Comment]: L535-542 this is unclear to me. Is this the method used to produce Supp Fig 1 or\nFig 4f?? Also, the actual calculation here is underspecified. What is x_i,k? I think i the species,\nj is the site, but what is k?\n[Reply]: Yes, this is the method used to produce Fig 4f. We revised the wording in the revision\nas:\nWe also conducted a “species variation decomposition”, whereby we quantified and\ncompared the coefficients of determination from comparing observed (unmodified) leaf\nN, P, and N:P values with modified values. Modifications followed two alternative\nassumptions, reflecting hypotheses in Fig. 1. The first assumes that variations arise only\nacross species (no within-species variations), and values x of species i and observation\ni,k\nk are replaced by the respective species’ mean, . The second assumes that variations\narise only within species and each record is “nor\n𝑥\nm𝑖alized” such that the resulting species\nmean is equal to the global mean .\n𝑥\nReferences cited above:\nAckerly, D. D & Cornwell, W. K. A trait-based approach to community assembly: partitioning\nof species trait values into within- and among-community components. Ecol. Lett. 10:\n135-145.\nCheng, J. H., Chu, P. F., Chen, D. M., Bai, Y. F. Functional correlations between specific leaf\narea and specific root length along a regional environmental gradient in Inner Mongolia\ngrasslands. Funct. Ecol. 30: 985-997 (2016).\nDong, N., et al. Components of leaf-trait variation along environmental gradients. New Phytol.\n228: 82-94 (2020).\nGallagher, R. V. & Leishman, M. R. Contrasting patterns of trait-based community assembly\nin lianas and trees from temperate Australia. Oikos 121: 2026-2035 (2012).\nKooyman, R., Cornwell, W., Westoby, M. Plant functional traits in Australian subtropical rain\nforest: partitioning within-community from cross-landscape variation. J. Ecol. 98: 517-\n525 (2010).\nOttaviani, G., Tsakalos, J. L., Keppel, G., Mucina, L. Quantifying the effects of ecological\nconstraints on trait expression using novel trait-gradient analysis parameters. Ecol. Evol.\n8: 435-440 (2018).\nZhou, J. H., Cieraad, E., van Bodegom, P. M. Global analysis of trait–trait relationships within\nand between species. New Phytol. 233: 1643-1656 (2022).\nREVIEWERS' COMMENTS\nReviewer #1 (Remarks to the Author):\nThank you for the thorough revisions the author made to address my concerns. The quality\nof the article has significantly improved. I believe the paper has now reached the standard\nfor publication in Nature Communications.\nReviewer #3 (Remarks to the Author):\nThis work provides a balanced analysis of contrasting views for leaf chemical traits along\nenvironmental gradients concerning the relative roles of specie identity and turnover versus\nwithin species plasticity. I thank the authors for their consideration and response to my\nreview. I found the ms much improved. it should be an important paper for our field. I have\njust a couple of comments (not necessarily to be acted upon) and few minor suggestions.\nThere is no doubt that flexible Machine learning methods such as Random Forests and\nBoosted Regression Trees are very useful. Perhaps most useful when exploring partial\ndependence plots to (a) identify most important predictors, (b) learn about the functional\nform of relationships, (c) identifying which interactions are most important. This is\nimportant for our field to understand how and why the RF models explain so much more\nthan the LMMs and how do we improve LMMs, which being parametric offer advantages for\ninterpretation and generality over machine learning methods. Perhaps a comment in the\ndiscussion?\nThanks for the work to address my comments about TGA. I accept your conclusions. I am\nstill privately concerned about TGA. But I do not wish to impeded progress on publishing\nthis work. Using your own code, I found that shortening the range of species greatly reduced\nthe slope of Ti,k on Tbar_k. intermediate simulations with a mix of niche turnover and\nplasticity produced a range of slopes. And I am still trying to make sense of them and cannot\nrid myself of the sense of circularity. Indeed, your simulation of Perfect plasticity encodes\ndirect dependence of species trait observation on the site trait mean. Anyhow, this is not\nthe concern of the current paper. You do find quite a strong degree of plasticity; on that we\ncan agree. Good luck with your important work. I enjoyed thinking about it.\nMinor edits\nL40 suggest “compiled” for “completed”\nL250 suggest “compiled” for “completed”\nL397ff can this be rewritten, something like this:\n“The role of environment in both filtering species and influencing traits within species,\ncombined with…”\nReviewer #3 (Remarks on code availability):\nI installed and ran the supplementary code for the rejoinder. it was clear,\ncommented well. I did not run the main scripts. Looking at the organisation of the other scripts,\nthe good organisation and clear coding lead me to be fairly confident in them. running full\nanalyses would take considerable time. I noted that the data are not all available in the Github\nrepo indicated.\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\n[Comment]: Thank you for the thorough revisions the author made to address my concerns.\nThe quality of the article has significantly improved. I believe the paper has now reached the\nstandard for publication in Nature Communications.\n[Reply]: We thank the reviewer for his or her comments which give us opportunity to improve\nour manuscript.\nReviewer #3 (Remarks to the Author):\n[Comment]: This work provides a balanced analysis of contrasting views for leaf chemical\ntraits along environmental gradients concerning the relative roles of specie identity and turnover\nversus within species plasticity. I thank the authors for their consideration and response to my\nreview. I found the ms much improved. it should be an important paper for our field. I have just\na couple of comments (not necessarily to be acted upon) and few minor suggestions.\n[Reply]: We thank the reviewer for his or her appreciation of the importance of our work.\n[Comment]: There is no doubt that flexible Machine learning methods such as Random Forests\nand Boosted Regression Trees are very useful. Perhaps most useful when exploring partial\ndependence plots to (a) identify most important predictors, (b) learn about the functional form\nof relationships, (c) identifying which interactions are most important. This is important for our\nfield to understand how and why the RF models explain so much more than the LMMs and\nhow do we improve LMMs, which being parametric offer advantages for interpretation and\ngenerality over machine learning methods. Perhaps a comment in the discussion?\n[Reply]: Thank you for this important suggestion. We agree with you about the advantage of\nRF. Hence, we mentioned this point in the Discussion as follows:\nRandom Forest models learn interactive and non-linear effects of multiple environmental\nvariables more effectively, detect trait–environment relationships beyond those arising from\nspecies composition (Fig. 3b, e, h), and yield superior results compared with linear regression\nmodels in out-of-sample evaluations when fitted to leaf N and P data at hand (Fig. 3a, d, g).\nHowever, fitted model coefficients of RF are not always directly interpretable — in contrast to\ncoefficients of LMMs. (l. 328-334).\nIn Fig. 2, we also provided information about feature importance and effect magnitudes and\ndirection, although the feature importance order was established by feature elimination rather\nthan by analysing a single fitted model, and the effect magnitude and direction was from LMMs,\nnot the RF.\nSpecific comments:\n[Comment]: Thanks for the work to address my comments about TGA. I accept your\nconclusions. I am still privately concerned about TGA. But I do not wish to impeded progress\non publishing this work. Using your own code, I found that shortening the range of species\ngreatly reduced the slope of Ti,k on Tbar_k. intermediate simulations with a mix of niche\nturnover and plasticity produced a range of slopes. And I am still trying to make sense of them\nand cannot rid myself of the sense of circularity. Indeed, your simulation of Perfect plasticity\nencodes direct dependence of species trait observation on the site trait mean. Anyhow, this is\nnot the concern of the current paper. You do find quite a strong degree of plasticity; on that we\ncan agree. Good luck with your important work. I enjoyed thinking about it.\n[Reply]: We are also interested in conducting further exploration about TGA in the next work.\nThank you very much for your suggestion about this important point.\nMinor edits:\n[Comment]: L40 suggest “compiled” for “completed”\n[Reply]: Revised according to this suggestion.\n[Comment]: L250 suggest “compiled” for “completed”\n[Reply]: Revised according to this suggestion.\n[Comment]: L397ff can this be rewritten, something like this:\n“The role of environment in both filtering species and influencing traits within species,\ncombined with…”\n[Reply]: This sentence was revised as follows:\nThe role of environment in both filtering species and influencing traits within species, combined\nwith the widespread use of linear mixed effects models with limited function of distinguishing\nstatistical independence of phylogeny and environmental effects, predisposed published\nanalyses to miss strong effects of environmental variables (l. 423-427).\nReviewer #3 (Remarks on code availability):\n[Comment]: I installed and ran the supplementary code for the rejoinder. it was clear,\ncommented well. I did not run the main scripts. Looking at the organisation of the other scripts,\nthe good organisation and clear coding lead me to be fairly confident in them. running full\nanalyses would take considerable time. I noted that the data are not all available in the Github\nrepo indicated.\n[Reply]: We thank the reviewer for his or her appreciation of the codes. All the codes and data\nare available at Github (https://github.com/stineb/leafnp_data) and Zenodo\n(https://doi.org/10.5281/zenodo.11071944). Source data are available at Zenodo\n(https://zenodo.org/records/11105009). In addition, figures can be re-created directly using\nintermediate data objects stored at GitHub (https://github.com/stineb/leafnp_data) in\ndata/data_published_figures."}]

Decision letter