# Empirical support for the biogeochemical niche hypothesis in forest trees

Jordi Sardans, Helena Vallicrosa, Paolo Zuccarini, Gerard Farré-Armengol, Marcos Fernández-Martínez, Guille Peguero, Albert Gargallo-Garriga, Philippe Ciais, Ivan Janssens, Michael Obersteiner, et al.

► To cite this version:

Jordi Sardans, Helena Vallicrosa, Paolo Zuccarini, Gerard Farré-Armengol, Marcos Fernández-Martínez, et al.. Empirical support for the biogeochemical niche hypothesis in forest trees. Nature Ecology & Evolution, 2021, 5 (2), pp.184-194. 10.1038/s41559-020-01348-1. hal-04228261

HAL Id: hal-04228261

https://hal.science/hal-04228261v1

Submitted on 26 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

This item is the archived peer-reviewed author-version of:

Empirical support for the biogeochemical niche hypothesis in forest trees

Reference:

Sardans Jordi, Vallicrosa Helena, Zuccarini Paolo, Farre-Armengol Gerard, Fernandez-Martinez Marcos, Peguero Guillermo, Gargallo-Garriga Albert, Ciais Philippe, Janssens Ivan, Obersteiner Michael, .... - Empirical support for the biogeochemical niche hypothesis in forest trees

Nature Ecology & Evolution - ISSN 2397-334X - Berlin, Nature research, 5(2021), p. 184-194

Full text (Publisher's DOI): https://doi.org/10.1038/S41559-020-01348-1

To cite this reference: https://hdl.handle.net/10067/1749910151162165141

<figure><img src="image_3.png" /></figure>
1 Empirical support for the Biogeochemical Niche Hypothesis in forest trees

2

3 Sardans<sup>1,2</sup>, J., Vallicrosa<sup>1,2</sup>, H., Zuccarini<sup>1,2</sup>, P. Farré-Armengol<sup>1,2</sup> G., Fernández-Martínez<sup>3</sup>, M.,

4 Peguero<sup>1,2</sup>, G., Gargallo-Garriga<sup>1,2</sup>, A., Ciais<sup>4</sup> P., Janssens<sup>3</sup> I.A., Obersteiner<sup>5</sup> M., Richter<sup>6</sup>, A.,

5 Peñuelas<sup>1,2</sup>, J.

6

7 <sup>1</sup> CSIC, Global Ecology Unit CREAF-CEAB-UAB, Bellaterra, 08193 Catalonia, Spain.

8 <sup>2</sup> CREAF, Cerdanyola del Vallès, 08193 Catalonia, Spain.

9 <sup>3</sup> PLECO (Plants and Ecosystems), Department of Biology, University of Antwerp, B-2610 Wilrijk, Belgium.

10

11 <sup>4</sup>Laboratoire des Sciences du Climat et de l'Environnement, IPSL, 91191 Gif-sur-Yvette, France.

12 <sup>5</sup>International Institute for Applied Systems Analysis (IIASA), Ecosystems Services and Management,

13 Schloss platz 1, A-2361 Laxenburg, Austria.

14 <sup>6</sup> University of Vienna, Terrestrial Ecosystem Research, Department of Microbiology and Ecosystem

15 Science, 1090 Vienna, Austria

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34
The possibility of using the elemental compositions of species as a tool to identify species/genotype niche remains to be tested at a global scale. We investigated the relationships between the foliar elemental compositions (elementomes) of trees at a global scale with phylogeny, climate, N deposition and soil traits. We analyzed foliar N, P, K, Ca, Mg and S concentrations in 25544 trees of 227 species. Shared ancestry explained 60-94% of the total variance of foliar nutrient concentrations and ratios, whereas current climate, atmospheric N deposition and soil type together explained 1-7%, consistent with the biogeochemical niche (BN) hypothesis that expects each species to have a specific need and use of each bioelement. The remaining variance was explained by the avoidance of nutritional competition with other species and the natural variability within species. The BN hypothesis is thus able to quantify species-specific tree niches and their shifts in response to environmental changes.

The "niche" theory is fundamental to ecology, because niches are both drivers and consequences of evolutionary processes¹⁻⁴. 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 Christian (1986)¹, Wright et al. (2006)⁵ and Swanson et al. (2015)⁶, among many others, where the niche is directly associated with a multivariate space. Accurately measuring 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, e.g. niche regeneration⁷ and functional niche⁸⁻¹⁰. A more general and easier method for defining and measuring species niches using field data, however, is needed.

The recently proposed Biogeochemical Niche (BN) hypothesis¹¹⁻¹³ 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, a product of long-term evolutionary processes, so each species should have a specific morphological structure and functionality (from gene expression to physiological processes). Fundamental biological processes (e.g. 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
70 genetic variability of a population and also on the phenotypical plasticity of individuals to respond to
71 environmental shifts throughout their lives. The BN hypothesis allows us to detect plasticity at two levels:
72 within populations due to the intraspecific variability in elemental composition of a set of individuals of the
73 same species living under the same environmental conditions at a specific time (e.g. due to individual
74 genotypic differences or different ontogenetic stages), and at the individual level (phenotypic plasticity) by
75 indicating how each individual of a population varies its elemental composition when environmental
76 conditions shift¹³.

77 The BN hypothesis is useful for representing the ecological niche of each species in a hyper-
78 dimensional volume generated by different bio-elemental concentrations and stoichiometric relationships,
79 which could be simply and practically tested by a combined chemical and mathematical approach using
80 multivariate and phylogenetic analyses¹³. The position of each species in the hyper-dimensional volume
81 can shift with time and changing environmental conditions depending on the degree of species-specific
82 stoichiometric plasticity (plasticity component), but the BN should also tend to maintain its own identity
83 relative to the BNs of other species (homeostatic component)¹³. BN space at a specific time should
84 therefore be a consequence of historical and current trends toward maximizing fitness in response to abiotic
85 and biotic circumstances such as trophic relationships and water, light or nutrient availabilities and fluxes¹³-
86 ¹⁵. The various levels of plasticity among species can be detected because homeostatic species will occupy
87 a smaller volume in multi-dimensional space, and plastic species will occupy a larger volume¹³. For
88 example, comparing the movement/expansion/contraction of the BNs of two populations of different
89 species or genotypes submitted to the same environmental shifts will thus provide information about their
90 levels of BN plasticity¹³. The temporal shifts of the BN of a species or population can also be calculated, so
91 we can follow the signatures of evolution on the BN. BNs also allow us to describe and quantify the
92 expansion, contraction and extinction of niches and the appearance of new niches (e.g. when new species
93 colonize an ecosystem)¹³. The BN hypothesis, however, has been experimentally tested only at small
94 spatial and phylogenetic scales¹¹⁻¹⁷.

95 We aimed to test the BN hypothesis in forest trees across all forest biomes and continents (Fig. S1).
96 Our study was based on the elementome of leaves, the plant organ where most compounds (from those
97 allocated to growth or reproduction to those allocated to energy metabolism, defense or storage) are
98 synthesized and where photosynthesis, the most crucial plant function, occurs. Leaves thus constitute a key
99 organ in plant functioning, and we can define the ‘species biogeochemical niche’ for each species in its
100 environmental circumstances by analyzing changes in foliar elemental composition. We built a global data
101 set to test the extent to which (i) shared ancestry, (ii) abiotic factors (e.g. climate, N deposition and soil
102 traits) and (iii) biotic factors (e.g. composition of the community inhabited by a tree) affect the BN. We
103 hypothesized that each species would have a different need and use for each bio-element to optimize
104 function and thus that shared ancestry would account for a large part of the foliar elemental composition in
105 a wide set of data from a broad spectrum of species. Species, however, also evolve to some degree during
106 fluctuating environmental conditions, so some of the variability of 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 “elementome”¹³, to define species-specific “niche differences” in a tangible and measurable way, providing a new valuable tool for establishing and identifying species niches.

112

## 113 Results

114

### 115 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 PCA axes (Table S1, Figs. 1-4 and S2-S14), with Pagel's λ values in several variables >0.5.

These variables tended to be more similar in the clades of more recently separated species than in the clades of more phylogenetically distant species, thus demonstrating that the divergence of the values of the variables among the clades over time were largely and significantly driven by evolutionary processes. For example, foliar N concentrations were 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 (7479 datapoints representing 138 species) than in the general database without foliar C concentrations (Tables S2 and S3, Figs. S15-S19). 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 (Table S2).

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

138

### 139 Abiotic factors

Climate, N deposition and soil type.

141
Phylogeny explained most of the variance of the foliar elemental compositions. Phylogeny explained 58.7-91.7% (mean 80.9%) of the variance of the foliar concentrations 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 principal component analysis (PCA) (Table S7). Species explained 1.4-14.6% (mean 5.58%) of the variance of the foliar concentrations of the six bio-elements, 0.4-28.1% (mean 8.1%) of the variance in pair-wise ratios and 0.6-13.8% (mean 5.07%) of the variance in the scores of the first three PC scores. Inheritance thus explained 73.3-93.6% (mean 86.5%) of the variance in the foliar concentrations 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 PC.

The climatic variables and N deposition, independently from their effect on current species distribution, explained 0.43-4.2% (mean 2.26%) of the variance in the foliar concentrations 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 PC (Table S7). Mean annual precipitation (MAP) and mean annual temperature (MAT) were the most important climatic variables explaining the variances of the elemental concentrations, with some relevant exceptions. For example, the variability of foliar N concentration was partially due to positive correlations with MAP, MAT, N deposition and mean annual solar radiation (Table S6). 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 (Table S7).

Higher precipitation was associated with lower foliar concentrations of metal elements (K, Ca and Mg) and P and was positively correlated with foliar N and S concentrations (Table S7). MAT was correlated positively with foliar N and metal (K, Ca, 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 percentage of the variances in elemental compositions and ratios (Table S8), ranging between 0.1% in the foliar P:S ratio to 2.0% in foliar K concentrations. A PCA, however, indicated that trees growing in different soil types occupied significantly different areas of the 2D plot of the first two PC axes (Fig. 6) and that this distribution was mainly explained by phylogeny ($R^2$ = 0.72, 0.84 and 0.87 for the PC1, PC2 and PC3 scores, respectively) (Table S8). 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 (Fig. S20), with intermediate ratios of foliar N:P, N:K and P:K (Fig. S21). Trees growing in Spodosols, very common in boreal and alpine coniferous forests, had the highest P and lowest K foliar concentrations (Fig. S20) and thus the highest foliar P:K ratios (Fig. S21). Trees growing in Oxisols, typical of wet tropical forests, had the second highest foliar N concentrations and the lowest foliar P concentrations (Fig. S20), the highest foliar N:P and N:K ratios and the lowest foliar P:K ratio (Fig. S21).

## Biotic factors
179 Competition of foliar elementomes between coexisting species.

180 Foliar N and P concentrations, the N:P ratio and the scores for the first two PCA axes for species frequently
181 shifted significantly when the distributions of two species overlapped. These differences usually explained
182 little of the variance of the foliar variables, but they were statistically significant (Tables S9-S15). For
183 example, foliar N and P concentrations, the N:P ratio and the scores for the first two PCA axes differed
184 significantly when *P. abies* grew in areas with and without *Q. robur* (*Fig. S22a*), although the percentage
185 of explained variance was ≤3.5% (Table S9). These shifts were reciprocal. Foliar P concentrations, N:P
186 ratios and the PC1 scores also differed significantly between subsets of *Q. robur* growing in areas with and
187 without *P. abies* (*Fig. S22a*, Table S10), explaining a maximum of only 5% of the total variance. Foliar P
188 concentrations, N:P ratios and the PC1 scores differed significantly between *A. alba* coexisting or not with
189 *Q. petraea* (*Fig. S22b*, Table S11), explaining a maximum of 5% of the total variance, and foliar N and P
190 concentrations, N:P ratios and the PC1 scores differed significantly between *Q. petraea* coexisting or not
191 with *A. alba* (Table S12), explaining a maximum of 4% of the total variance. Foliar N concentrations and
192 the PC1 scores differed significantly between *Fagus sylvatica* growing in areas with and without *P.
193 sylvestris* (*Fig. S22c*, Table S13), explaining a maximum of 2% of the total variance. Foliar N and P
194 concentrations, N:P ratios and the PC1 scores differed significantly between *P. sylvestris* growing in areas
195 with and without *F. sylvatica* (*Fig. S22c*, Table S14), explaining a maximum of 1% of the total variance.
196 The percentage of variance explained for some of the variables was much higher in some cases, e.g. 17%
197 for the foliar N concentration of *P. sylvestris* growing with or without *Q. robur* (Table S14) and 35% for
198 the foliar P concentration of *Larix decidua* growing with or without *Q. robur* (Table S15).
199

200 Discussion

201 BN size and phylogeny

202 The use of the foliar concentrations of several bio-elements clearly separated the species in the
203 hypervolume generated by the corresponding multivariate analyses, as previously reported only in smaller
204 studies<sup>11-18</sup>. The FDA plot (showing 95.8% of the total explained variance) clearly separated the BNs of the
205 species, with *Pinus* species having positive values on Root 1 and *Quercus* species having more negative
206 values, also consistent with the positive link between phylogenetic distance and species-specific BN
207 identity among species.
208

208 N, P, K, S, Ca and Mg contribute differentially to plant metabolic and physiological functions and
209 to cells, tissues and organs. We therefore expected that different species, as evolutionary products, would
210 have different optimal elemental compositions (elementomes)<sup>11-14</sup>. The results of our analysis of the global
211 set of 227 of the most representative tree species worldwide (163, 58 and 6 from tropical, temperate and
212 boreal biomes, respectively) strongly support this hypothesis. The results also indicated that the foliar BNs
213 of the species became more similar as their phylogenetic distance decreased. These results are fully
214 consistent with Kerkhoff et al.<sup>19</sup>, who also found a consistent and significant phylogenetic signal in N and P
concentrations in plant organs in a set of 1287 plant species. The small number of published studies of BNs, all including fewer species, have 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 ratios²⁰. The great majority of the studies nevertheless found significant relationships between species elemental composition and taxonomic and phylogenetic distance²⁰. Similar results have also been obtained in ionomic studies²¹.

Phylogeny and species, as proxies of overall genomic difference, however, did not explain 100% of the 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 (e.g. 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 species⁵,²²,²³. Several other factors, such as species migration, changes of species interactions (e.g. with herbivorous or parasitic species) and climatic convergence can increase the speed of evolutionary convergence among species in different clades²⁴⁻²⁶. Distant clades could thus evolve under new, more similar, environmental conditions, favoring a trend toward convergence in functionality and thus in elementome. The results nevertheless indicated that evolutionary processes have significantly contributed to the differences in foliar elementomes that originated during species diversification, directly explaining 57-94.5% (averaging 85.7%) of the variance of their foliar concentrations and ratios. Anacker and Strauss (2016)²⁷ also recently reported that the 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 the variance of the variables studied, perhaps due to the recent divergent evolution of more proximal taxonomical species recently 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 gradients²⁸⁻³¹, 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 increased²⁸⁻³⁴. Not all studies, however, have detected clear patterns of N:P ratios among or within the climatic areas³³. Townsend et al. (2007)³⁴ and Lovelock et al. (2007)³⁵, for example, found no relationship between the 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, e.g. 16-25 and 5-35% reported by Yuan and Chen (2009)³¹ and Reich and Oleksyn (2004)²⁸, 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) (Table S6).

Our results thus indicated 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 frequently observed higher N and lower P concentrations in plants toward equatorial latitudes. We also identified a general and significant relationship of the foliar concentrations of the main macronutrients with N deposition. The lower foliar metal concentrations with higher levels of N deposition are completely consistent with the higher leaching of soil bases associated with N deposition and with the bases competition for plant absorption with ammonium³⁶. The higher foliar N and S concentrations with more N deposition are also due to the consequent higher availability of soil N and S³⁷,³⁸. Interestingly, we also identified a global 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. The 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 availability³⁹,⁴⁰. Lower foliar P concentrations under higher N loads, however, have been associated with a stronger P limitation⁴¹⁻⁴³. Our results thus indicated that N deposition in forests at the global scale 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 the variance (0.1-2%) of tree foliar composition and stoichiometry. In fact, soil type and its capacity to supply some of the most important bio-elements to plants are 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 easily weathered minerals such as Ca²⁺ and Mg²⁺⁴⁴,⁴⁵, 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 with 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 (Fig. S15) 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
288 foliar P:K ratio (Fig. S16). These results for Oxisols were consistent with recent observations of low foliar
289 P concentrations and high foliar N:P ratios in wet tropical forests²⁸,³². To the best of our knowledge, our
290 global study is the first to associate high foliar N:K and low foliar P:K ratios with wet tropical forests.
291 Relationships between foliar BN and soil type along natural gradients have recently been observed⁴⁶, but
292 these relationships may not be as strong as expected and may not necessarily be universal. Ordoñez et al.
293 (2009)²⁹ observed that the concentrations of some elements and ratios were correlated between soil type
294 and photosynthetic tissues but others were not.

295

296 Competition among coexisting species

297 The foliar N and P concentrations, the N:P ratio and the scores for the first two PCA axes for species
298 shifted significantly when the distributions of two species overlapped. These differences usually explained
299 little of the variance of the foliar elemental concentrations and ratios but were significant (Tables S8-S14).
300 For example, foliar N and P concentrations, the N:P ratio and the scores for the first two PCA axes differed
301 significantly when *P. abies* grew in areas with and without *Q. robur*, although the percentage of explained
302 variance was ≤3.5% (Table S8). These shifts were reciprocal and occurred for the pairs of the tree species,
303 which were the most dominant in Europe (Tables S8-S14). These results of field analyses are consistent
304 with the results of an experiment in seminatural grasslands where the target species shifted their elemental
305 compositions depending on the neighboring species¹⁷.
306

306

307 Homeostasis versus plasticity

308 Intraspecific variability explained a significant amount (2-20%) of the total random variability of all
309 elemental concentrations and ratios (Tables S6 and S7). Species that have evolved in highly fluctuating
310 environments are expected to have a larger capacity of functional and/or morphological shifts and thus
311 require a more plastic stoichiometry than species that have evolved in a more stable environment¹³,³⁹. Our
312 results confirmed these expectations: the range of values in the FDA that represented the sizes of species-
313 specific BNs was larger (*P*<0.001 along Roots 1 and 2) for the species subjected to less climatic stress (Q.
314 *robur* and *Q. petraea*) than for the Mediterranean species adapted to drought (*P. halepensis*, *P. pinaster* and
315 *Q. ilex*) (Fig. 5, Tables S5 and S6). These results indicated a trade-off between adaptation to being
316 competitive in a stable environment versus being successful in a more fluctuating environment. Different
317 levels of environmental stress cause a continuum of strategies between homeostasis and plasticity. Species
318 growing in more stressful environments, with poor resource availability, have less BN plasticity than
319 species growing in less stressful and richer environments⁴⁷,⁴⁸.

320

321 Conclusions

322 The results of this study provide clear support for the BN hypothesis¹³. First, each species had a different
323 BN, with a significant trend of larger differences in BNs as phylogenetic distance and evolutionary time
324 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 the 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.

335

## 336 METHODS

337

338 **Data acquisition**

339 **Foliar Data.** We gathered 25 544 datapoints of foliar N, P, K, Ca, Mg and S concentrations, expressed as percent dry weight. These data corresponded to 227 tree species at a global scale, including all latitudes and ecosystems. We only considered tree species with more than three locations. The data were obtained from 192 publications (Table S1) and inventories such as the Catalan Forest Inventory⁴⁹. We also gathered and used a subset of 7479 datapoints with 138 species that contained information of foliar C concentration, in addition to foliar N, P, K, Mg, Ca and S concentrations, for identifying possible differences in the analyses with or without C concentrations. All data had been obtained from leaves using comparable and homologated analytical methods (see pages 9-13, ICP forests manual Sampling and Analysis of Needles and Leaves, http://icp-forests.org/manual.htm). 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 analyzed using mixed models, with database as a fixed factor and country as a random factor. No significant differences were found. All foliar samples had been collected between 1990 and 2015. We only used data from georeferenced plots. **Fig. S1** shows the distribution of the plots.

352

353 **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 database⁵⁰, with a resolution of 1 km² at the equator: minimum average temperature, maximum average temperature, average solar radiation, maximum wind speed, average wind speed, mean water vapor pressure, mean annual temperature (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, mean annual precipitation (MAP), mean precipitation of the wettest month, mean precipitation of the driest month, mean precipitation seasonality, mean precipitation of the wettest quarter,

361
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 data⁵¹⁻⁵⁴.

The data for the deposition of atmospheric N and P were obtained from Global Threats to Human Water Security and River Biodiversity⁵⁵, with a resolution of 1 km² 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 km² 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 software⁵⁶. We thereby obtained a phylogenetic tree containing a selection of species from PhytoPhylo, an available megaphylogeny of vascular plants⁵⁷. We used the *read.tree* and *drop tip* functions from the R *ape* package⁵⁸ to load the PhytoPhylo tree and remove all species that were not in our database.

We used the *phylosig* function from the R *phytools* package⁵⁹ to test for phylogenetic signals in the foliar elemental compositions of the 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 of phylogenetically independent contrasts relative to tip shuffling randomization⁶⁰. We chose to analyze the phylogenetic signals in the data using Pagel's λ assumption, based on a study by Münkemüller et al. (2012)⁶¹ comparing the advantages and disadvantages of various methods for estimating phylogenetic signals. Pagel's λ method can provide reliable measurements of effect size and can discriminate between more complex models of trait evolution (such as polygenic organismic traits)⁶². 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)⁶³, unlike other methods.

We also used the *contMap* function of the *phytools* package to graphically reconstruct the values of the traits that had a phylogenetic signal across our phylogeny. We used the *ape* package⁵⁰ to load the phylogenetic tree and select the species in it. The *contMap* function estimates the ancestral characters at internal nodes using maximum likelihood and assuming Brownian motion as a model for trait evolution⁵⁵ and then interpolates the ancestral condition along the branches of the tree⁵⁴

**BN size.** “BN size” is another interesting trait when comparing taxa. We thus conducted a functional discriminant analysis (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: *Quercus petraea*, *Q. robur*, *Q. ilex*, *Pinus 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 the optimal separation between groups (here 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 (here the foliar N, P, K, S, Ca and Mg concentrations and their pairwise ratios)⁶⁵. 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 the variance in the foliar elemental composition of the various species. We then analyzed the roles of phylogeny and species in the dispersion of the canonical scores on the root axes of the FDA using Bayesian phylogenetic linear mixed models and the MCMCglmm package⁶⁶ in R. Phylogeny and species were included as random factors: the phylogenetic term accounted for variability in the shared ancestry, and the species term accounted for species-specific traits independent of the shared ancestry.

411

**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 the bioelements, their ratios and PC scores (of the PCA of all elemental foliar concentrations and their pair-wise ratios) using Bayesian phylogenetic linear mixed models and the MCMCglmm package⁶⁶ in R. We used MAT, MAP, annual radiation, mean annual vapor pressure deficit VPD, range of diurnal temperatures and N deposition as fixed effects. Phylogeny and species were included as random factors: the phylogenetic term accounted for the variability in the shared ancestry, and the species term accounted for species-specific traits independent of the shared ancestry. Both random factors together thus accounted for the variance explained by heritability. We repeated these analyses using soil type instead of climate and N deposition as fixed effects.

421

**Analysis of the relationships between foliar elemental composition and soil type.** We analysed the differences of the foliar variables among the various soil types (taxonomic orders). A soil map was generated using the R *raster* and *rgdal* packages to obtain the soil classifications 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 bioelements, their ratios and PC scores (of the PCA of all elemental foliar concentrations and their pair-wise ratios) using Bayesian phylogenetic linear mixed models and the MCMCglmm package⁶⁶ 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 the variance explained by heritability, as described previously.

We performed PCAs of 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 analyzed 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 package⁶⁶ 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 the FDA scores. Model parameters (soil types) with non-overlapping 95%
credible intervals were considered to differ significantly.

440

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 (EUFGIS) 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 between 25 and  
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 analyzed 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 choose 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  
the global distributions for each species of each pair that overlapped or not with the distribution of the other  
species. An ANOVA compared the 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).

455

456 Acknowledgements

The authors would like to acknowledge the financial support from the European Research Council Synergy grant ERC-SyG-2013-610028 IMBALANCE-P, the Spanish Government grant PID2019-110521GB-I00 and the Catalan Government grant SGR 2017-1005. M.F.-M. is a postdoctoral fellow of the Research Foundation – Flanders (FWO).

461

462 References

1. Tracy, C.R. & Christian, K.A. Ecological relations among space, time, and thermal niche axes. *Ecology*, **67**, 609-615 (1986).

465 2.Peterson, A. T., Soberon, J. & Sanchez-Cordero, V. Conservatism of ecological niches in evolutionary time. *Science* **285**, 1265-1267 (1999).

467 3.Hirzel, A. H., Hausser, J., Chessel, D. & Perrin, N. Ecological-niche factor analysis: How to compute habitat-suitability maps without absence data? *Ecology* **83**, 2017-2036 (2002).

469 4.Sexton, J. P., Montiel, J., Shay, J. E., Stephens, M. R. & Slatyer, R. A. Evolution of ecological niche breadth. *An. Rev. Ecol. Evol. Syst.* **48**, 183-206 (2017).

471 5.Wright, J. W., Davies, K. F., Lau, J. A., McCall, A. C. & McKay, J. K. Experimental verification of ecological niche modelling in a heterogeneous environment. *Ecology* **87**, 2433-2439 (2006).

473 6. Swanson, H. K. et al. A new probabilistic method for quantifying n-dimensional ecological niches and niche overlap. *Ecology* **96**, 318-324 (2015).

475 7.Grubb, P. J. The maintenance of species-richness in plant communities: the importance of the regeneration niche. *1977. Biol. Rev.* **52**, 107-145 (1977).

477 8.Herrel, A., Spithoven, L., Van Damme, V. & De Vree F. Sexual dimorphism of head size in *Gallotia galloti*: testing the divergence hypothesis by functional analyses. *Funct. Ecol.* **13**, 289-297 (1999).

478
479 9.Mouillot, D. et al. Niche overlap estimates based on quantitative functional traits: a new family of non-parametric indices. *Oecologia* **145**, 345-353 (2005).

480 10.Kraft, N., Valencia, R. & Ackerly D.D. Functional traits and niche based tree community assembly in an amazonian forest. *Science* **322**, 580-582 (2008).

481 11.Peñuelas, J., Sardans, J., Ogaya, R.,& Estiarte, M. Nutrient stoichiometric relations and biogeochemical niche in coexisting plant species: effect of simulated climate change. *Pol. J. Ecol.* **56**, 613-622 (2008).

485 12.Peñuelas J. et al. Faster returns on "leaf economics" and different biogeochemical niche in invasive compared with native plant species. *Glob. Change Biol.* **16**, 2171-2185 (2010).

487 13.Peñuelas, J. et al. The bioelements, the elementome and the "biogeochemical niche". *Ecology* **100**, e02652 (2019).

488 14.Sardans, J. et al. Factors influencing the foliar elemental composition and stoichiometry in forest trees in Spain. *Persp. Plant Ecol. Evol. Syst.* **18**, 52-69 (2016).

490 15.Sardans, J. et al. Foliar elemental composition of European forest tree species associated with evolutionary traits and present environmental and competitive conditions. *Glob. Ecol. Biogeogr.* **24**, 240-255 (2015).

492 16.Urbina, I. et al. Shifts in the elemental composition of plants during a very severe drought. *Environ. Exp. Bot.* **111**, 63-73 (2015).

494 17.Urbina, I. et al. Plant community composition affects the species biogeochemical niche. *Ecosphere* **8**, e01801 (2017).

496 18. White, P.J. et al. Testing distinctness of shoot ionomes of angiosperm families using the Rothamsted Park grass continuous hay experiment. *New Phytologist* **196**, 101-109 (2012).

498 19.Kerkhoff, A. J., Fagan, W. F., Elser, J. J., Enquist, B. J. Phylogenetic and growth form variation in the scaling of nitrogen and phosphorus in the seed plants. *The Am. Nat.* **168**, E103-E122 (2006).

500 20. Sardans, J. Janssens, I. A., Ciais P., Obersteiner, M., Peñuelas J. Recent advances, approaches and future research in ecological stoichiometry: towards a novel holistic theory in ecology. *Biogeochemistry*. Submitted.

502 21. Neugebauer, K. et al. Variation in the Angiosperm Ionome. *Physiol. Plant.* **163**, 306-322 (2018).

503 22.Gillman, L. N., Keeling, D. J., Gardner, R.C. & Wright, S.D. Faster evolution of highly conserved in tropical plants. *J. Evol. Biol.* **23**, 1327-1330 (2010).

505 23.Puurtinen, M. et al. Temperature-dependent mutational robustness can explain faster molecular evolution at warm temperatires, affecting speciation rate and global patterns of species diversity. *Ecography* **39**, 1025-1033 (2016).

507 24.Kellner, A., Ritz, C. M., Schlittenhaedt, P. & Hellwig, F. H. Genetic differentiation in the genus Lithps L. (Ruschoideae, Aizoaceae) revelas a high level of convergent evolution and reflects geographic distribution. *Plant Biol.* **13**, 368-380 (2011).

510 25.Jwa, N. S. & Hwang, B. K. Convergent evolution of pathogen effectors toward reactive oxygen species signaling networks in plants. *Front. Plant Sci.* **8**, 1687 (2017).

512 26.Molina-Montenegro, M.A et al. Is the causes of plant invasions the result of rapid adaptive evolution in seed traits? evidence a latitudinal rainfall gradient? *Front. Plant Sci.* **9**, 208 (2018).

514 27.Anacker, B. L. & Strauss, S. Y. Ecological similarity is related to phylogenetic distance between species in a cross-niche field transplant experiment. *Ecology* **97**, 1807-1818 (2016).

516 28.Reich, P. B. & Oleksyn, J. Global patterns of plant leaf N andP in relation to temperature and latitude. *Proc. Natl. Acad.Sci. USA* **101**, 11001-11106 (2004).

518 29.Ordoñez, J. C. et al. A global study of relationships between leaf traits, climate and soil measures of nutrientfertility. *Glob. Ecol. Biogeogr.* **18**, 137-149 (2009).

520 30.Kerkhoff, A. J., Enquist, B. J., Elser, J. J. & Fagan, W. F. Plantallometry, stoichiometry and the temperature-dependenceof primary productivity. *Glob. Ecol. Biogeogr.* **14**, 585-598 (2005).

522 31.Yuan Z. Y. & Chen, H. Y. H. Global trends in senesced-leafnitrogen and phosphorus. *Glob. Ecol. Biogeogr.* **18**, 532-542 (2009).

524 32.Vitousek, P. M., Porder, S., Houlton, B. Z. & Chadwick, O. A. Terrestrial phosphorus limitation: mechanisms, implications and nitrogen-phosphorus interactions. *Ecol. Appl.* **20**, 5-15 (2010).

526 33.McGroddy, M. E., Daufresne, T. & Hedin, L. O. Scaling of C/N/P stoichiometry in forest worldwide: implications afterrestrial Redfield-type ratios. *Ecology* **85**, 2390-2401 (2004).

528 34.Townsend, A. R., Cleveland, C. C., Asner, G. P. & Bustamante, M. M. C. Controls over foliar N:P ratios in tropical rainforest. *Ecology* **88**, 107-118 (2007).

530 35.Lovelock, C. E., Feller, I. C., Ball, M. C., Ellis, J. & Sorell, B. Testing the growth rate vs. geochemical hypothesis for latitudinal variation in plant nutrients. *Ecol. Lett.* **10**, 1154-1163 (2007).

532 36.Marschner, H. *Mineral nutrition of higher plants*, London. Academic Press (1995).

533 37.Zhang, Y., Mathur, R., Bash, J. O., Hogrefe, C., Xing, J. & Roselle, S. J. Log-term trends in total inorganic nitrogen and sulfur deposition in US from 1990 to 2010. *Atm. Chem. Phys.* **18**, 9091-9106 (2018).

535 38.Horn, K, J. et al. Growth and survival relationships of 71 tree species with nitrogen and sulfur deposition across the conterminous U.S. *PlosOne* **14**, e0212984 (2019).
537 39.Papanikolaou, N., Britton, A. J., Helliwell, R. C. & Johnson, D. Nitrogen deposition, vegetation burning and climate warming act independently on microbial community structure and enzyme activity associated with decomposing litter in low-alpine heath. *Glob. Change Biol.* **16**, 3120–3132 (2010).

540 40.Marklein, A. R. & Houlton, B. Z. Nitrogen inputs accelerate phosphorus cycling rates across a wide variety of terrestrial ecosystems. *New Phytol.* **193**, 696–704 (2012).

542 41.Sardans, J. et al. Foliar and soil concentrations and stoichiometry of nitrogen and phosphorous across European *Pinus sylvestris* forests: relationships with climate, N deposition and tree growth. *Funct. Ecol.* **30**, 676-689 (2016).

544 42.Peñuelas J. et al. Human-induced nitrogen-phosphorus imbalances alter natural and managed ecosystems across the globe. *Nature Com.* **4**, 2934 (2013).

546 43.Peñuelas, J. et al. Increasing atmospheric CO₂ concentrations correlate with declining nutritional status of European forests. *Com. Biol.* **3**, 125 (2020).

548 44.Ahmad, N. & Mermut, A. Vertisols and Tachnologies for their Management. Development in Soil Science 24. Elsevier. Amsterdam. 302 pp (1996).

550 45.Nishiue, A., Nanzyo, M., Kanno, H. & Takahashi, T. Properties and classification of volcanic ash soils around Lake Kuwanuma on the eastern footslope of Mt. Funagata in Miyagi prefecture, northeastern Japan. *Soil Sci. Plant Nutr.* **60**, 848-862 (2014).

553 46. de la Riva, E. G. et al. Biogeochemical and Ecomorphological Niche Segregation of Mediterranean Woody Species along a Local Gradient. *Fron. Plant Sci.* **8**, 1242 (2017).

555 47.Yu, Q. et al. Stoichiometry homeostasis of vascular plants in the inner Mongolia grassland. *Oecologia* **166**, 1–10 (2011).

557 48.Sardans, J., Albert Rivas-Ubach, A. & Peñuelas, J. The elemental stoichiometry of aquatic and terrestrial ecosystems and its relationships with organismic lifestyle and ecosystem structure and function: a review. *Biogeochemistry* **111**, 1–39 (2012).

560 49.Gracia, C., Burriel, J. A., Ibàñez, J. J., Mata, T. & Vayreda, J. *Inventari ecològic i forestal de Catalunya. Mètodes*, Vol. 9. CREF, Bellaterra (2004).

562 50.Fick, A. E. & Hijmans, R. J. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. *Int. J. Climatol.* **37**, 4302-4315 (2017).

564 51.Lang, R. Verwitterung und Bodenbildung als Einfuehrung in die Bodenkunde. Schweizerbart Science Publishers, Stuttgart (1920).

566 52.Köppen, W. Klassification der Klimate nach Tempertur, Niederschlag and Jahreslauf. *Petermanns Geog. Mitt.* **64**, 243–248 (1918).

569 53.De Martonne, E. Nouvelle carte mondiale de l'indece d'aridité. Annales de Géographie 51, 242-250 (1942).

570 54.Emberger, L. La vegetation de la région méditerranéenne, essai d'une classification des groupements vegetaux. *Rev. Gén. Bot.* **648** (1930).

571 55.Vorösmarty, C. J. et al. Global threats to human water security and river biodiversity. *Nature* **467**, 555-561 (2010).

573 56.R Development Core Team. *R: a language and environment for statistical computing*. R Foundation for Statistical Computing, Vienna. Available at: http://www.R-project.org/. (2011).

575 57. Qian, H., Jin, Y. An updated megaphylogeny of plants, a tool for generating plant phylogenies, and an analysis of phylogenetic community structure. *J. Plant Ecol.* **9**, 233–239 (2016).

577 58.Paradis, E., Claude, J. & Strimmer, K. (2004) APE: analyses o fphylogenetics and evolution in R language. *Bioinformatics* **20**, 289–290 (1999).

579 59.Revell, L. J. Phytools: An R package for phylogenetic comparative biology (and other things). *Methods Ecol. Evol.* **3**, 217–223 (2012).

581 60. Blomberg, S.P., Garland, T., Ives, A.R. Testing for phylogenetic signal in comparative data: behavioral traits are more labile. *Evolution* **57**, 717–45 (2003).

583 61.Münkemüller, T. et al. How to measure and test phylogenetic signal. *Methods Ecol. Evol.* **3**, 743-756 (2012).

584 62.Pagel, M. Inferring the historical patterns of biological evolution. *Nature* **401**, 877–884 (1999).

585 63.Felsenstein, J. Phylogenies and the comparative method. *Am. Nat.* **125**, 1–15 (1985).

586 64.Revell, L. J. Two new graphical methods for mapping trait evolution on phylogenies. *Methods Ecol. Evol.* **4**, 754–759 (2013).

588 65.Ramos, S. D. S. & Liow, S. R. . Discriminant Functional Analysis. Ed: Chapelle, C.A. The Encyclopedia of Applied Linguistics, Chapter: Discriminant Function Analysis, Publisher: John Wiley & Sons (2012), Ed: Chapelle, C.A.66. Hadfield, J. D. MCMC methods for multi-response generalised linear mixed models: the MCMCglmm R package. *J. Stat. Softw.* **33**, 1–22 (2010).

592
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Phylogenetic diagrams of the foliar N concentration in the phylogenetic tree obtained using the *contMap* function of the *phytools* package in R, representing the value structure used to estimate the phylogenetic signals in the variables. The *contMap* function estimates the ancestral characters at internal nodes using maximum likelihood and assuming Brownian motion as a model for trait evolution⁶¹ and then interpolates the ancestral condition along the branches of the tree⁶². DW Dry weight.</figcaption></figure>

593

594

595

596

597

598

599

600

601

602

603

604

605

606

607

608
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Phylogenetic diagrams of the foliar P concentration in the phylogenetic tree obtained using the *contMap* function of the *phytools* package in R, representing the value structure used to estimate the phylogenetic signals in the variables. The *contMap* function estimates the ancestral characters at internal nodes using maximum likelihood and assuming Brownian motion as a model for trait evolution⁶¹ and then interpolates the ancestral condition along the branches of the tree⁶². DW = Dry weight</figcaption></figure>

609

610

611

612

613

614

615

616

617

618

619

620

621

622

623

624

625
<figure><img src="image_6.png" /><figcaption>Fig. 3. Phylogenetic diagrams of the foliar N:P ratio in the phylogenetic tree obtained using the *contMap* function of the *phytools* package in R, representing the value structure used to estimate the phylogenetic signals in the variables. The *contMap* function estimates the ancestral characters at internal nodes using maximum likelihood and assuming Brownian motion as a model for trait evolution⁶¹ and then interpolates the ancestral condition along the branches of the tree⁶².</figcaption></figure>

626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Phylogenetic diagrams of the PC1 scores in the phylogenetic tree obtained using the *contMap* function of the *phytools* package in R, representing the value structure used to estimate the phylogenetic signals in the variables. The *contMap* function estimates the ancestral characters at internal nodes using maximum likelihood and assuming Brownian motion as a model for trait evolution⁶¹ and then interpolates the ancestral condition along the branches of the tree⁶².</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Plot of the first two roots of the functional discriminant analysis using *Pinus pinaster*, *P. halepensis*, *Quercus 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. The plot depicts the 95% confidence interval of the mean corresponding to each species and the variable loadings on Roots 1 and 2.</figcaption></figure>

663

664

665

666

667

668
<figure><img src="image_9.png" /></figure>

669

670

671

672

673

**Fig. 6.** Plot of the 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. Different letters on the arrows at the right indicate significant differences (*P*<0.05). Error bars indicate the 95% confidence intervals.
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691

692
<figure><img src="image_10.png" /></figure>
<figure><img src="image_11.png" /></figure>

length=176.162

0.022 trait value 0.329
<figure><img src="image_12.png" /></figure>

length=176.162

3.723

trait value

57.372
<figure><img src="image_13.png" /></figure>
<figure><img src="image_14.png" /></figure>
<figure><img src="image_15.png" /></figure>