Main

The ‘niche’ theory is fundamental to ecology, because niches are both drivers and consequences of evolutionary processes1,2,3,4. The concept is easy to understand theoretically: each species tends to occupy a particular position along the gradients of all abiotic and biotic variables that define and determine fitness. This view is consistent with the niche concepts defined by Tracy and Christian1, Wright et al.5 and Swanson et al.6, among many others, where the niche is directly associated with a multivariate space. Accurate measurement of the exact parameters of a niche, however, is challenging due to the large number of variables that affect organisms within ecosystems. Several approaches have been developed in recent decades to more easily manage this complex and multivariate concept—for example, niche regeneration7 and functional niche8,9,10. A more general and easier method to define and measure species niches using field data, however, is needed.

The recently proposed biogeochemical niche (BN) hypothesis11,12,13 incorporates most, if not all, niche parameters using species-specific elemental composition and stoichiometry. The assumptions underlying this hypothesis are based on the idea that each species is a unique genetic pool of individuals and a product of long-term evolutionary processes, so that each species should have a specific morphological structure and functionality (from gene expression to physiological processes). Fundamental biological processes (for example, growth, secondary metabolism, reproduction and storage of bio-elements) have distinct rates in different species depending on selection pressures, so different species must differentially allocate elements to various traits of tissues and organs. Each species should thus tend to have its own elemental composition and stoichiometry (homeostatic component of BN). The changing circumstances during the lives of organisms, however, should also determine a necessary phenotypical plasticity to allow the individuals of each species to adapt its functionality and morphology during their lives (plasticity component). This ability differs in extent and quality among species. BN plasticity depends on the current genotypic and thus genetic variability of a population, and also on the phenotypical plasticity of individuals to respond to environmental shifts throughout their lives. The BN hypothesis allows us to detect plasticity at two levels: (1) within populations due to the intraspecific variability in elemental composition of a set of individuals of the same species living under the same environmental conditions at a specific time (for example, due to individual genotypic differences or different ontogenetic stages), and (2) at the individual level (phenotypic plasticity), by indicating how each individual of a population varies its elemental composition when environmental conditions shift13.

The BN hypothesis is useful for representing the ecological niche of each species in a hyperdimensional volume generated by different bio-elemental concentrations and stoichiometric relationships, which could be simply and practically tested by a combined chemical and mathematical approach using multivariate and phylogenetic analyses13. The position of each species in the hyperdimensional volume can shift with time and changing environmental conditions, depending on the degree of species-specific stoichiometric plasticity (plasticity component), but the BN hypothesis should also tend to maintain its own identity relative to the BNs of other species (homeostatic component)13. BN space at a specific time should therefore be a consequence of historical and current trends toward maximization of fitness in response to abiotic and biotic circumstances such as trophic relationships and water, light or nutrient availabilities and fluxes13,14,15. The various levels of plasticity among species can be detected because homeostatic species will occupy a smaller volume in multidimensional space, and plastic species will occupy a larger volume13. For example, comparing the movement/expansion/contraction of the BNs of two populations of different species or genotypes submitted to the same environmental shifts will thus provide information about their levels of BN plasticity13. The temporal shifts of the BN of a species or population can also be calculated, so we can follow the signatures of evolution on the BN. BNs also allow us to describe and quantify the expansion, contraction and extinction of niches and the appearance of new niches (for example, when new species colonize an ecosystem)13. The BN hypothesis, however, has been experimentally tested only at small spatial and phylogenetic scales11,12,13,14,15,16,17.

We aimed to test the BN hypothesis in forest trees across all forest biomes and continents (Extended Data Fig. 1). Our study is based on the elementome of leaves, the plant organ where most compounds (from those allocated to growth or reproduction to those allocated to energy metabolism, defence or storage) are synthesized and where photosynthesis, the most crucial plant function, occurs. Leaves thus constitute a key organ in plant functioning, and we can define the ‘species biogeochemical niche’ for each species in its environmental circumstances by analysis of changes in foliar elemental composition. We built a global data set to test the extent to which (1) shared ancestry, (2) abiotic factors (for example, climate, N deposition and soil traits) and (3) biotic factors (for example, composition of the community inhabited by a tree) affect the BN. We hypothesized that each species would have a different need and use for each bio-element to optimize function, and thus that shared ancestry would account for a large part of the foliar elemental composition in a wide set of data from a broad spectrum of species. Species, however, also evolve to some degree during fluctuating environmental conditions, so some of the variability in elemental composition should be due to the phenotypic component and thus to current abiotic and biotic conditions such as climate, soil type, atmospheric N deposition and competition, which should account for another part of the variability of elemental composition. We thus expected to demonstrate the suitability of using the elemental compositions of organisms, the elementome13, to define species-specific ‘niche differences’ in a tangible and measurable way, providing a valuable tool for establishing and identifying species niches.

Results

Phylogeny and BN size

The analysis of Pagel’s λ identified significant phylogenetic signals in foliar N, P, K, S, Ca and Mg concentrations, N/Ca, N/Mg, P/Ca, P/Mg, P/S, K/Ca and K/S ratios and the scores of the first six principal component analysis (PCA) axes (Supplementary Table 1, Figs. 1–4, Extended Data Figs. 2–5 and Supplementary Figs. 1–9), with Pagel’s λ values in several variables >0.5. These variables significantly tended to be more similar in the clades of more recently separated species than in those of more phylogenetically distant species, thus demonstrating that the divergence of the values of the variables among clades over time was largely and significantly driven by evolutionary processes. For example, foliar N concentrations were significantly more similar among Pinaceae species than between Pinaceae and Fagaceae species, with Pinaceae species generally having lower values than Fagaceae species (Fig. 1). The phylogenetic signals using Pagel’s λ were generally more similar in the subset of the database that also contained information for foliar C concentrations (7,479 datapoints representing 138 species) than in the general database without foliar C concentrations (Supplementary Tables 2 and 3 and Supplementary Figs. 10–14). Of the 33 variables studied, only the foliar C/K ratio did not have a significant λ. Furthermore, these Pagel’s λ values were high (>0.6 for most of the nutrient variables) (Supplementary Table 2).

Fig. 1: Phylogenetic diagram of foliar N concentration (percentage of dry weight) (25,112 datapoints) in the phylogenetic tree.
figure 1

This was obtained using the contMap function of the phytools package in R, representing the value structure used to estimate phylogenetic signals in variables. The contMap function estimates the ancestral characters at internal nodes using maximum likelihood and assuming Brownian motion as a model for trait evolution61, and then interpolates the ancestral condition along the branches of the tree62.

Full size image
Fig. 2: Phylogenetic diagram of foliar P concentration (percentage of dry weight) (25,112 datapoints) in the phylogenetic tree.
figure 2

This was obtained using the contMap function of the phytools package in R, representing the value structure used to estimate phylogenetic signals in variables. The contMap function estimates the ancestral characters at internal nodes using maximum likelihood and assuming Brownian motion as a model for trait evolution61, and then interpolates the ancestral condition along the branches of the tree62.

Full size image
Fig. 3: Phylogenetic diagram of foliar N/P ratio (25,112 datapoints) in the phylogenetic tree.
figure 3

This was obtained using the contMap function of the phytools package in R, representing the value structure used to estimate phylogenetic signals in variables. The contMap function estimates the ancestral characters at internal nodes using maximum likelihood and assuming Brownian motion as a model for trait evolution61, and then interpolates the ancestral condition along the branches of the tree62.

Full size image
Fig. 4: Phylogenetic diagram of PC1 scores (25,112 datapoints) in the phylogenetic tree.
figure 4

This was obtained using the contMap function of the phytools package in R, representing the value structure used to estimate phylogenetic signals in variables. The contMap function estimates the ancestral characters at internal nodes using maximum likelihood and assuming Brownian motion as a model for trait evolution61, and then interpolates the ancestral condition along the branches of the tree62.

Full size image

The range of values in a functional discriminant analysis (FDA) (mainly along root 1, which explained 90.8% of the total variance) that represented the sizes of species-specific BNs was significantly larger for the species subjected to lower climatic stress (Quercus robur and Quercus petraea) than for Mediterranean species adapted to drought (Pinus halepensis, P. pinaster and Quercus ilex) (Fig. 5 and Supplementary Tables 4 and 5). These five species were the most abundant in our survey and were clearly separated in the multivariate space by significant distances (squared Mahalanobis distances; Supplementary Table 4) and in all foliar variables that significantly contributed to the separation of all species (Supplementary Table 5). We also detected a strong phylogenetic effect in the distribution of scores along the first three root axes of the FDA (Supplementary Table 6).

Fig. 5: Plot of the first two roots of functional discriminant analysis using P. pinaster, P. halepensis, Q. ilex, Q. petraea and Q. robur as dependent categorical grouping factors, and foliar N, P, K, S, Ca and Mg concentrations and pairwise ratios as continuous independent variables.
figure 5

The plot depicts the 95% confidence interval of the mean corresponding to each species and the variable loadings on roots 1 and 2.

Full size image

Abiotic factors

Climate, N deposition and soil type

Phylogeny explained significantly most of the variance in foliar elemental composition: 58.7–91.7% (mean 80.9%) of the variance in foliar concentration of the six bio-elements, 39.0–94.1% (mean 68.8%) of the variance in their pairwise ratios and 43.2–89.6% (mean 74.0%) of the variance in the scores of the first three axes of a PCA (Supplementary Table 7). Species significantly explained 1.4–14.6% (mean 5.58%) of the variance in foliar concentrations of the six bio-elements, 0.4–28.1% (mean 8.1%) of the variance in pairwise ratios and 0.6–13.8% (mean 5.07%) of the variance in the scores of the first three PC scores. Inheritance thus significantly explained 73.3–93.6% (mean 86.5%) of the variance in foliar concentration of the six bio-elements, 67.1–94.5% (mean 81.7%) of the variance in their pairwise ratios and 57–90.2% (mean 79.1%) of the variance in the scores of the first three PCs.

The climatic variables and N deposition, independently from their effect on current species distribution, significantly explained 0.43–4.2% (mean 2.26%) of the variance in foliar concentration of the six bio-elements, 0.27–7.4% (mean 5.05%) of the variance in their ratios and 0.35–11.0% (mean 4.23%) of the variance in the scores of the first three PCs (Supplementary Table 7). Mean annual precipitation (MAP) and mean annual temperature (MAT) were the most important climatic variables explaining the variance in elemental concentrations, with some relevant exceptions. For example, the variability in foliar N concentration was partially due to positive correlations with MAP, MAT, N deposition and mean annual solar radiation (Supplementary Table 6). Interestingly, MAT and mean solar radiation had contrasting relationships with MAP and N deposition on overall foliar elemental composition, as indicated by the PC1 scores (Supplementary Table 7).

Higher precipitation was significantly associated with lower foliar concentrations of metal elements (K, Ca and Mg) and P and was significantly and positively correlated with foliar N and S concentrations (Supplementary Table 7). MAT was correlated significantly and positively with foliar N and metal (K, Ca and Mg) concentrations and N/P ratio, and negatively with foliar P concentrations. N deposition was correlated positively with foliar N, P and S concentrations and negatively with foliar Ca, Mg and K concentrations.

Soil type explained a low but significant percentage of the variance in elemental composition and ratios (Supplementary Table 8), ranging between 0.1% for the foliar P/S ratio to 2.0% for foliar K concentrations. A PCA, however, indicated that trees growing in different soil types occupied significantly different areas of the two-dimensional plot of the first two PC axes (Fig. 6) and that this distribution was mainly explained by phylogeny (R2 = 0.72, 0.84 and 0.87 for PC1, PC2 and PC3 scores, respectively) (Supplementary Table 8). Trees growing on Inceptisols and Alfisols, typical of temperate forests, occupied a central position in this space. Trees growing in Alfisols, typical of wet and mesic temperate forests, had intermediate foliar N and K concentrations and the second highest foliar P concentration (Extended Data Fig. 6), with intermediate foliar N/P, N/K and P/K ratios (Extended Data Fig. 7). Trees growing in Spodosols, very common in boreal and alpine coniferous forests, had the highest P and lowest K foliar concentrations (Extended Data Fig. 6) and thus the highest foliar P/K ratios (Extended Data Fig. 7). Trees growing in Oxisols, typical of wet tropical forests, had the second highest foliar N concentrations and the lowest foliar P concentrations (Extended Data Fig. 6), the highest foliar N/P and N/K ratios and the lowest foliar P/K ratio (Extended Data Fig. 7).

Fig. 6: Plot of PCA cases and variables superimposed, defined by the first two components of the PCA, with foliar N, P, K, Ca, Mg and S concentrations as variables and with soil orders as cases.
figure 6

Different lowercase letters on the arrows indicate significant differences in the soil orders position along the two first components of the PCA (P < 0.05). The colours of the arrows highlight the soil order that represents. Error bars indicate the 95% confidence intervals.

Full size image

Biotic factors

Competition of foliar elementomes among coexisting species

Foliar N and P concentrations, the N/P ratio and the scores for the first two PCA axes for species frequently shifted significantly when the distributions of two species overlapped. These differences usually explained little of the variance in foliar variables, but they were statistically significant (Supplementary Tables 9–15). For example, foliar N and P concentrations, the N/P ratio and the scores for the first two PCA axes differed significantly when Pinus abies grew in areas with and without Q. robur (Extended Data Fig. 8a), although the percentage of explained variance was ≤3.5% (Supplementary Table 9). These shifts were reciprocal. Foliar P concentrations, N/P ratios and PC1 scores also differed significantly between subsets of Q. robur growing in areas with and without P. abies (Extended Data Fig. 8a and Supplementary Table 10), explaining a maximum of only 5% of total variance. Foliar P concentrations, N/P ratios and PC1 scores differed significantly between Abies alba coexisting or not with Q. petraea (Extended Data Fig. 8b and Supplementary Table 11), explaining a maximum of 5% of total variance, and foliar N and P concentrations, N/P ratios and PC1 scores differed significantly between Q. petraea coexisting or not with A. alba (Supplementary Table 12), explaining a maximum of 4% of total variance. Foliar N concentrations and PC1 scores differed significantly between Fagus sylvatica growing in areas with and without Pinus sylvestris (Extended Data Fig. 8c and Supplementary Table 13), explaining a maximum of 2% of total variance. Foliar N and P concentrations, N/P ratios and PC1 scores differed significantly between P. sylvestris growing in areas with and without F. sylvatica (Extended Data Fig. 8c and Supplementary Table 14), explaining a maximum of 1% of total variance. The percentage of variance explained for some variables was much higher in some cases—for example, 17% for the foliar N concentration of P. sylvestris growing with or without Q. robur (Supplementary Table 14) and 35% for the foliar P concentration of Larix decidua growing with or without Q. robur (Supplementary Table 15).

Discussion

BN size and phylogeny

Use of the foliar concentrations of several bio-elements significantly separated the species in the hypervolume generated by the corresponding multivariate analyses, as previously reported only in smaller studies11,12,13,14,15,16,17,18. The FDA plot (showing 95.8% of total explained variance) clearly separated the BNs of species, with Pinus species having positive values on root 1 and Quercus species having more negative values, also consistent with the positive link between phylogenetic distance and species-specific BN identity among species.

The elements N, P, K, S, Ca and Mg contribute differentially to plant metabolic and physiological functions and to cells, tissues and organs. We therefore expected that different species, as evolutionary products, would have different optimal elementomes11,12,13,14. The results of our analysis on the global set of 227 of the most representative tree species worldwide (163, 58 and 6 from tropical, temperate and boreal biomes, respectively) strongly support this hypothesis. The results also indicated that foliar BNs of the species became more similar as their phylogenetic distance decreased. These results are fully consistent with Kerkhoff et al.19, who also found a consistent and significant phylogenetic signal in N and P concentrations in plant organs in a set of 1,287 plant species. A small number of published studies of BNs, all including fewer species, has also reported significant organ or body stoichiometric dependence on taxonomy and/or phylogeny of plant and animal species, although not all studies detected links between species phylogeny and N/P ratios13,14,15,20. The great majority of the studies nevertheless found significant relationships between species' elemental composition and taxonomic and phylogenetic distance13,14,15,20. Similar results have also been obtained in ionomic studies21.

Phylogeny and species, as proxies of overall genomic difference, however, did not explain 100% of variability in the elementome. Some of the phylogenetic lines of distant clades may have been exposed to similar environmental conditions that would have driven parallel selection of the characteristics that determine elemental concentrations, consequently eliciting convergence to more similar elementomes than would be expected from their phylogenetic distance. In other words, species that are phylogenetically distant (for example, those that have developed on different continents but under current similar environmental conditions) may occupy a similar BN. For example, a change toward a warmer climate can increase the speed of evolution of several characters differently in different species5,22,23. Several other factors, such as species migration, changes in species interactions (for example, with herbivorous or parasitic species) and climatic convergence, can increase the speed of evolutionary convergence among species in different clades24,25,26. Distant clades could thus evolve under new, more similar, environmental conditions, favouring a trend toward convergence in functionality and thus in elementome. The results nevertheless indicated that evolutionary processes have significantly contributed to differences in foliar elementomes that originated during species diversification, directly explaining 57–94.5% (average 85.7%) of the variance in foliar concentrations and ratios. Anacker and Strauss27 also reported that niche differences among species increased with phylogenetic distance, again consistent with our results. Part of the inheritability factor that differs from phylogeny explained an average of about 7% of variance in the variables studied, perhaps due to the recent divergent evolution of more proximal taxonomical species adapting to distinct and divergent environmental shifts in their respective distribution areas.

Climate and N deposition

Several studies have reported trends in foliar N and P concentrations and N/P ratios in trees growing along climatic and latitudinal gradients28,29,30,31, but their results have not always fully agreed. Most studies have observed a general trend toward decreasing foliar and litter P concentrations and increasing N/P ratios as latitude decreased and MAT, MAP and length of growing season increased28,29,30,31,32,33,34. Not all studies, however, have detected clear patterns of N/P ratios among or within climatic areas33. Townsend et al.34 and Lovelock et al.35, for example, found no relationship between N/P ratio and either latitude or MAP in tropical areas. Even when a significant relationship was detected in these studies, climatic variables explained only a small fraction of the variation in foliar elements—for example, 16–25 and 5–35% reported by Yuan and Chen31 and Reich and Oleksyn28, respectively, depending on the variable. These studies used linear models that considered only climatic variables and N deposition without phylogeny and species as random factors. The variance explained by climatic variables in our study decreased in our Bayesian analyses when we added species as a random factor (0.3–11.0%, with a mean of 4.3% among all variables) (Supplementary Table 6).

Our results thus indicate significant relationships of foliar nutrient composition with N deposition and basic climatic traits such as MAP and MAT at a global scale. A decrease in foliar P concentrations and an increase in foliar N concentrations and N/P ratios with increasing MAT are consistent with the higher N and lower P concentrations in plants frequently observed toward equatorial latitudes. We also identified a general and significant relationship of foliar concentrations of the main macronutrients with N deposition. Lower foliar metal concentrations with higher levels of N deposition are totally consistent with the higher leaching of soil bases associated with N deposition and with competition among bases for plant absorption with ammonium36. Higher foliar N and S concentrations with increased N deposition are also due to the consequent higher availability of soil N and S37,38. Interestingly, we also identified a globally positive correlation between higher levels of N deposition and higher foliar P concentrations when local studies reported all types of results, from increases to decreases in foliar P concentrations. Increases in P concentrations have been associated with higher capacities of plants and microbes to mobilize and take up more P due to higher N availability39,40. Lower foliar P concentrations under higher N loads, however, have been associated with stronger P limitation41,42,43. Our results thus indicated that N deposition in forests at the global scale significantly tends to increase foliar P concentrations but also N/P ratios, thus generally trending toward more P-limited forests.

Soil type

Soil type explained a modest but significant amount of variance (0.1–2%) in tree foliar composition and stoichiometry. In fact, soil type and its capacity to supply some of the most important bio-elements to plants is partially due to historic and current climatic conditions. Species foliar elementomes were consistent with the traits of the various soil types (Fig. 6). Trees growing in Inceptisols and Alfisols, typical of temperate forests, occupied a central position in the PCA space, suggesting a more balanced elemental composition than trees growing in other soil types. Trees growing in Andisols and Vertisols, two soil types rich in readily weathered minerals such as Ca2+ and Mg2+ (refs. 44,45), typically had higher than average foliar Ca and Mg concentrations. Andisols are volcanic soils that are frequently rich in Fe-Mg silicates and in anorthite, a Ca feldspar. Ertisols are characterized by high concentrations of expandable clays such as vermiculite and montmorillonite that are also rich in Mg and Ca, respectively. Trees growing in Spodosols (typical of sandy soils) had the lowest concentrations of Ca, Mg and K, which could be linked to the high leaching of basic cations in these acidic soils and, consequently, the low content of exchangeable complexes and slow mineralization. Trees growing in Spodosols also had the highest foliar P concentrations (Supplementary Fig. 10) and thus the highest foliar P/K ratios. Trees growing in Oxisols (wet tropical forests) had the second highest foliar N concentrations, the lowest foliar P concentrations, the highest foliar N/P and N/K ratios and the lowest foliar P/K ratio (Supplementary Fig. 11). These results for Oxisols are consistent with recent observations of low foliar P concentrations and high foliar N/P ratios in wet tropical forests28,32. Our global study thus associates high foliar N/K and low foliar P/K ratios with wet tropical forests. Relationships between foliar BN and soil type along natural gradients have recently been observed46, but these relationships may not be as strong as expected and may not necessarily be universal. Ordoñez et al.29 observed that the concentrations of some elements and ratios were correlated between soil type and photosynthetic tissues but that others were not.

Competition among coexisting species

Foliar N and P concentrations, N/P ratio and scores for the first two PCA axes for species shifted significantly when the distributions of two species overlapped. These differences usually explained little of the variance in foliar elemental concentrations and ratios but were significant (Supplementary Tables 8–14). For example, foliar N and P concentrations, N/P ratio and scores for the first two PCA axes differed significantly when P. abies grew in areas with and without Q. robur, although the percentage of explained variance was ≤3.5% (Supplementary Table 8). These shifts were reciprocal and occurred for pairs of tree species that were the most dominant in Europe (Supplementary Tables 8–14). These results from field analyses are consistent with those of an experiment in seminatural grasslands where the target species shifted their elemental compositions depending on neighbouring species17.

Homeostasis versus plasticity

Intraspecific variability explained a significant amount (2–20%) of the total random variability of all elemental concentrations and ratios (Supplementary Tables 6 and 7). Species that have evolved in highly fluctuating environments are expected to have a greater capacity of functional and/or morphological shifts and thus require a more plastic stoichiometry than species that have evolved in a more stable environment13,39. Our results confirmed these expectations: the range of values in the FDA that represented the sizes of species-specific BNs was higher (P < 0.001 along roots 1 and 2) for species subjected to less climatic stress (Q. robur and Q. petraea) than for the Mediterranean species adapted to drought (P. halepensis, P. pinaster and Q. ilex) (Fig. 5 and Supplementary Tables 5 and 6). These results indicated a trade-off between adaptation to being competitive in a stable environment versus being successful in a more fluctuating environment. Different levels of environmental stress cause a continuum of strategies between homeostasis and plasticity. Species growing in more stressful environments and with poor resource availability have less BN plasticity than those growing in less stressful and richer environments47,48.

Conclusions

The results of this study provide clear support for the BN hypothesis13. First, each species had a different BN, with a significant trend of larger differences in BNs as phylogenetic distance and evolutionary time increased. Recent evolutionary convergence due, for example, to recent adaptation of distant clades to similar soil or climatic environments, however, indicated that the differences in BN among species could not be fully resolved by phylogenetic analyses alone. Second, environmental factors such as climate and soil type also explained an important part of the intraspecific variance in BN. These effects were moderate but significant and independent of taxonomy. Each species could be represented by its specific space in the hypervolume generated by multivariate analysis of its foliar elemental composition and stoichiometry (elementome), so its specific plasticity was observed in the shift of its space in response to environmental changes. Third, coexisting, competing species tended to have distinct BNs to minimize competitive pressure. Fourth, a trade-off between adaptations to being competitive in a stable environment versus being successful in a more fluctuating environment generated a continuum of strategies between homeostasis and plasticity.

Methods

Data acquisition

Foliar data

We gathered 23,962 datapoints of foliar N, P, K, Ca, Mg and S concentrations, expressed as percentage dry weight. These data corresponded to 227 tree species at a global scale, including all latitudes and ecosystems. We considered only those tree species with more than three locations. The data were obtained from 192 publications (Supplementary Table 1) and inventories such as the Catalan Forest Inventory49. We also gathered and used a subset of 7,479 datapoints with 138 species that included information on foliar C concentration, in addition to foliar N, P, K, Mg, Ca and S concentrations, for identification of possible differences in the analyses with or without C concentrations. All data were obtained from leaves using comparable and homologated analytical methods (see the ICP forests manual, Sampling and Analysis of Needles and Leaves, http://icp-forests.net/page/icp-forests-manual). The N, P, K, S, Ca and Mg pairwise ratios were calculated on a mass basis. Nutrient concentrations for the same species from different databases were analysed using mixed models, with database as a fixed factor and country as a random factor; no significant differences were found. All foliar samples were collected between 1990 and 2015. We used data derived from georeferenced plots only; Supplementary Fig. 1 shows the distribution of those plots.

Data for climate, soil and N and P deposition

Climatic and soil data were added to the foliar stoichiometric data using the raster package in R (v.2.6–7). These data were obtained from the WorldClim 2.0 database50, with a resolution of 1 km2 at the equator: minimum average temperature, maximum average temperature, average solar radiation, maximum wind speed, average wind speed, mean water vapour pressure, MAT, mean diurnal range, isothermality, temperature seasonality, maximum temperature of the warmest month, minimum temperature of the coldest month, annual temperature range, mean temperature of the wettest quarter, mean temperature of the driest quarter, mean temperature of the coldest quarter, MAP, mean precipitation of the wettest month, mean precipitation of the driest month, mean precipitation seasonality, mean precipitation of the wettest quarter, mean precipitation of the driest quarter, precipitation of the warmest quarter and precipitation of the coldest quarter. This climatic model was calculated for a long meteorological time series (1970–2000) based on interpolated values of climatic data provided by meteorological stations throughout the territory and adjusted to the observed topography. Five aridity indices were calculated using the climatic data51,52,53,54.

The data for the deposition of atmospheric N and P were obtained from Global Threats to Human Water Security and River Biodiversity55, with a resolution of 1 km2 at the equator. Soil taxonomies (order and suborder) were obtained from the USDA Global Soils Region Map (https://www.nrcs.usda.gov/wps/portal/nrcs/detail/soils/use/?cid=nrcs142p2_054013), which provides a resolution of 1 km2 at the equator.

Phylogenetic and statistical analyses

Phylogenetic signal

We prepared a phylogenetic tree containing the species in our database to test for phylogenetic signals using R statistical software56. We thereby obtained a phylogenetic tree containing a selection of species from PhytoPhylo, an available megaphylogeny of vascular plants57. We used the read.tree and drop.tip functions from the R ape package58 to load the PhytoPhylo tree and removed all species that were not in our database.

We used the phylosig function from the R phytools package59 to test for phylogenetic signals in the foliar elemental compositions of species and therefore to determine the extent to which foliar N, P, K, S, Ca and Mg concentrations, pairwise ratios and PCA scores had phylogenetic signals. The phylosig function calculates statistics of a phylogenetic signal (Pagel’s λ) and P values based on the variance in phylogenetically independent contrasts relative to tip shuffling randomization60. We chose to analyse phylogenetic signals in the data using Pagel’s λ assumption, based on a study by Münkemüller et al.61 comparing the advantages and disadvantages of various methods for estimation of phylogenetic signals. Pagel’s λ method can provide reliable measurements of effect size and discriminate between more complex models of trait evolution (such as polygenic organismic traits)62. Mean λ in Pagel’s method does not change as the number of species in a phylogeny increases, and is recommended for large phylogenies with >50 species (or taxa)63, unlike other methods.

We also used the contMap function of the phytools package to graphically reconstruct the values of traits with a phylogenetic signal across our phylogeny. We used the ape package50 to load the phylogenetic tree and select the species included in it. The contMap function estimates the ancestral characters at internal nodes using maximum likelihood and assuming Brownian motion as a model for trait evolution55, and then interpolates the ancestral condition along the branches of the tree64

BN size

Another interesting trait when comparing taxa is BN size. We thus conducted FDA to determine whether different but closely related species typical of different environments (from more to less climatic stress) tended to have different BN sizes. We compared five of the most important forest species in Europe: Q. petraea, Q. robur, Q. ilex, P. pinaster and P. halepensis. All five species were represented at 600–800 sites in our database across their distributions. FDA is a multivariate analysis that derives optimal separation between groups (in this case, the different sets of individuals of each species) by maximizing between-group variance and minimizing within-group variance of the set of independent continuous variables used in the analysis (in this case, foliar N, P, K, S, Ca and Mg concentrations and their pairwise ratios)65. We compared the range of scores in the first two roots of the FDA as a measure of the size of the ‘niche space’ of each species. The first two roots of the FDA explained 95.8% of variance in the foliar elemental composition of the various species. We then analysed the roles of phylogeny and species in the dispersion of canonical scores on the root axes of the FDA using Bayesian phylogenetic linear mixed models and the MCMCglmm package66 in R. Phylogeny and species were included as random factors: the phylogenetic term accounted for variability in the shared ancestry while the species term accounted for species-specific traits independent of shared ancestry.

Analysis of the relationships of foliar elemental composition with climatic variables and N deposition

We tested the effects of climate and N deposition on the foliar concentrations of bio-elements, their ratios and PC scores (from the PCA of all elemental foliar concentrations and their pairwise ratios) using Bayesian phylogenetic linear mixed models and the MCMCglmm package66 in R. We used MAT, MAP, annual radiation, mean annual vapour pressure deficit, range of diurnal temperatures and N deposition as fixed effects. Phylogeny and species were included as random factors: the phylogenetic term accounted for variability in shared ancestry while the species term accounted for species-specific traits independent of the shared ancestry. Both random factors together thus accounted for variance explained by heritability. We repeated these analyses using soil type rather than climate and N deposition as fixed effects.

Analysis of the relationships between foliar elemental composition and soil type

We analysed differences in foliar variables among the various soil types (taxonomic orders). A soil map was generated using the R packages raster and rgdal to obtain soil classification for each sample location. We chose soil type (orders of soil taxonomy) as the most accurate taxon at the pixel scale in the USDA Global Soils Region Map. No data were found for Gelisols or Aridisols. We tested the effect of soil order on the foliar concentration of bio-elements, their ratios and PC scores (from the PCA of all elemental foliar concentrations and their pairwise ratios) using Bayesian phylogenetic linear mixed models and the MCMCglmm package66 in R. We used soil orders as fixed effects. Phylogeny and species were included as random factors: the phylogenetic and species terms were introduced as random factors accounting for variance explained by heritability, as described previously.

We performed PCAs for foliar N, P, K, Ca, Mg and S concentrations and N/P ratios to further explore the relationships between trees growing under different soil types and their overall elemental compositions. We then analysed the scores of the PC1 and PC2 axes to detect differences in overall foliar elemental composition depending on the order of the soil in which they grew, using Bayesian phylogenetic linear mixed models and the MCMCglmm package66 in R, in which the first three PCA axes were the response variables and soil order was the fixed predictor. Phylogeny and species were included as random factors, as in the previous analysis of FDA scores. Model parameters (soil types) with non-overlapping 95% credible intervals were considered to differ significantly.

Analysis of differences in species foliar elemental composition and stoichiometry between populations growing in different communities with different species compositions

We used the map of species distribution in the European Information System on Forest Genetic database, http://portal.eufgis.org/data/. We compared the foliar N, P and K concentrations of pairwise species with comparable co-occurring and non-co-occurring surfaces with an overlapping distribution of 25–75%. We established sets of individuals in the overlapping area between the compared species and both areas where only one of the species was present. The data were analysed in R using the packages raster (v.3.4.3), rgeos (v.3.4.4), maptools (v.3.4.3), maps (v.3.4.3), rworldmap (v.3.4.4), ggmap (v.3.4.3) and rworldxtra (v.3.4.4). We used these tools to select species with large distributions and many datapoints in our database (600–800) and combined species in several possible pairs. We compared the two portions of global distributions for each species of each pair that overlapped or not with the distribution of the other species. We used analysis of variance to compare N, P and K concentrations and PC1 and PC2 scores (from the PCA of the six bio-elements and their pairwise ratios) for each species inside and outside the overlapping zone (with or without competition between the two species, respectively).

Reporting Summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.