Climate, soil resources and microbial activity shape the distributions of mountain plants based on their functional traits

We demonstrate tight associations between microbial decomposition activity, plant functional traits associated to different resource acquisition strategies and plant distributions. This highlights the importance of plant–soil linkages for mountain plant distributions. These results are crucial for biodiversity modelling in a world where both climatic and soil systems are undergoing profound and rapid transformations.


Introduction
Human-induced global change not only alters climate, but also soil properties and functioning (Rillig et al. 2019). Climate and soil jointly influence plant species distributions through their effect on plant growth, survival and reproduction (Buri et al. 2017). Moreover, there is a growing body of evidence that the interplay of plant functional traits with these environmental gradients shapes plant distributions (Pollock et al. 2012, Carboni et al. 2018. The framework of the 'fast-slow' plant economics spectrum predicts that the environment selects plants with specific functional traits, resulting in the spatial segregation of plants with different ecological strategies (Choler 2005, Reich 2014). Favorable environmental conditions select for functional trait values related to an exploitative strategy (i.e. rapid nutrient acquisition and fast growth, Aerts 1999, Reich 2014. Plants with exploitative traits lead to a soil organic matter that is easier to decompose, which increases nutrient availability, and in return favors exploitative plants. On the contrary, limiting environmental conditions promote a conservative strategy (i.e. slow plant growth but high nutrient retention capacity, Aerts 1999). Conservative plant traits induce low decomposability of soil organic matter, which slows down nutrient release in the soil, and thus selects plants with conservative trait values (Wardle and Bardgett 2002, Grigulis et al. 2013, Legay et al. 2016). Although it is well known that plants and soils are linked by feedback loops at local scales, the influence of plant-soil linkages on plant distributions is under-studied at large biogeographic scales.
One of the best-studied environmental drivers of plant success is climate (e.g. temperature, Grace 1987, Choler 2018. In mountain environments, plant success depends on the ability to grow and reproduce quickly within the short growing season. Moreover, climatic extremes such as freezing events or drought can damage plant tissues and thus limit survival (Grace 1987). In addition to climate, soil physicochemical properties, such as pH, soil organic matter and soil C/N drive plant growth because they are related to nutrient availability (Lee 1998, Buri et al. 2017. However, measuring nutrient supply for plants is not straightforward at large spatial scales. Plant nutrient supply depends on the decomposition rate of soil organic matter by microbial communities. Many different microorganisms, mainly bacteria and fungi, produce extracellular enzymes (exoenzymes) to recycle the complex macromolecules constituting the debris of dead living organisms such as plants themselves. Different exoenzymes target different molecules and at the end of the decomposition chain, some of them enable the release of soluble monomers and finally nutrient mineralization (Sinsabaugh et al. 2008, Burns et al. 2013. Nutrient supply for plants should thus be related to the activities of exoenzymes catalyzing the terminal reactions of decomposition. The sum of the activities of the main exoenzymes in the top-soil (hereafter total EEA) has been suggested as a relevant indicator of the microbial decomposition activity (Burns et al. 2013) that makes nutrients available for plant uptake.
In soils, both the overall nutrient availability and the element stoichiometry (here we focus on the balance between carbon, nitrogen and phosphorus) are driven by the quantity and nutritional quality of plant litter inputs, and affect plant growth in return (Li et al. 2019). Microbes modulate the production of element-specific exoenzymes according to soil element stoichiometry (Sinsabaugh et al. 2008) Thus, the relative investment of microbial community in the acquisition of each type of nutrients can be used as a proxy of nutrient limitation (e.g. high activity of N-targeting exoenzymes, EEN, in comparison to total EEA, EEN/total EEA, reflects a low availability in environmental nitrogen relative to carbon and phosphorus). Finally, microbes compete with plants for nutrients dissolved in the soil solution. The relative demand of the microbial community for nitrogen compared to phosphorus can be assessed through ratios of exoenzymes (i.e. EEN/EEP, Piton et al. 2019). For example, a high EEN/ EEP ratio means that plant competition with microbes will be stronger for nitrogen than for phosphorus, and might reduce nitrogen availability for plants. In sum, considering soil exoenzymes and their ratios promises an important way forward to include soil nutrient availability driven by reciprocal plant-soil linkages in plant distribution models. Moreover, accounting for plant functional traits related to nutrient acquisition strategies will allow accounting for the fact that conservative species respond in a different way to nutrient limitation than exploitative species.
In this paper, we aim to disentangle how climate, soil physico-chemical properties and microbial decomposition activity (the two latter are hereafter summarized as 'soil properties') influence the distributions of 44 plant species in the French Alps, while accounting for the modulating effects of plant traits. Our approach is based on a hierarchical multispecies distribution model and addresses the two following questions: 1) what are the relative influences of climate, soil physico-chemical properties and microbial decomposition activity (measured via ratios of exoenzymatic activities) in explaining plant distributions in mountain ecosystems? 2) Do species with different functional traits respond differently to gradients of climate and soil properties?

Study sites and species data
We studied the distributions of plant species along 18 elevation gradients in the French Alps. The gradients belong to the long-term observatory ORCHAMP (<www.orchamp. osug.fr>, Supplementary material Appendix 1) and were selected according to the following criteria: 1) continuous gradient from about 900 m to 3000 m (Supplementary material Appendix 1 Fig. A1-A), 2) exposure and slope along the gradient is homogeneous, 3) vegetation is typical for the elevation stages with forests dominating the lower parts and alpine meadows the higher parts of the gradient and 4) all gradients together are representative for the environmental and topographical variability of the French Alps (Fig. 1, Supplementary material Appendix 1 Fig. A1-B). Between 2016 and 2018, a minimum of five sampling plots (30 × 30 m, see Fig. 1 for details) were installed along each elevation gradient, with an average of 200 m elevation difference, resulting in 99 plots in total (Supplementary material Appendix 1 Fig.  A1). Each plot included a subplot for collecting plant data and three subplots for soil sampling (Fig. 1). Plant and soil data were always sampled in the same year. All plant species in the subplot were identified by professional botanists. For the analyses, we only kept the 44 species for which we had sufficient presences (i.e. > 20 plots, Supplementary material Appendix 2 Table A1) to build reliable species distribution models (Guisan et al. 2017). For each plant species, we selected plant height and three leaf traits, the ratio between leaf carbon and nitrogen contents (leaf C/N), leaf dry matter content (LDMC) and specific leaf area (SLA) because they represent important aspects of plant nutrition and light competition (Cornelissen et al. 2003). Traits were retrieved from our own trait database, and averaged at the species level. On average 10 individuals were sampled over different environmental conditions typical for the French Alps following standard protocols (Supplementary material Appendix 2 Table  A1; Cornelissen et al. 2003).

Climatic data
To best characterize species' climatic niches, we calculated a set of bioclimatic variables as averages over 1988-2018. More specifically, we used solar radiation and growing degree days (GDD, Choler 2018) to represent the average length and intensity of the growing season, intensity of freezing events (freezing degree days, FDD, Choler 2018) and water stress (climatic water deficit, CWD) in each plot (see Table 1 for details). These variables were calculated from the SAFRAN-SURFEX/ISBA-Crocus-MEPRA reanalysis (Durand et al. 2009, Vannier 2012. This corresponds to a version almost identical to the dataset made available by Vernay et al. (2019). The last corrections made before the publication of this dataset are sufficiently minor to not modify the conclusions of this study. A key advantage of this model is that it combines observed weather station data and the output of numerical weather prediction models, and provides a joint product addressing meteorological and snow conditions in mountainous regions, based on large-scale topographical features (e.g. elevation, slope, aspect). This is of particular importance in mountain environments, where the starting point and the length of the growing season is determined by the presence of snow (Choler 2005, Carlson et al. 2015.

Soil sampling, soil physico-chemical properties and exoenzymatic activities
We sampled soil in each of the three 2 × 2 m soil subplots located 5 m below the inner transect of each sampling plot ( Fig. 1). In each subplot, we sampled ten soil cores (10 cm depth × 5 cm diameter) and immediately pooled and homogenized them. We sieved the resulting composite samples with 5.6 mm-sieves. We froze 2.75 g per sample for subsequent exoenzyme analyses (soil sub-samples B), 5 g to quantify water content of each sample, and kept the rest for physicochemical analyses (soil sub-samples A).
From soil sub-samples A, we measured pH, soil C/N, total nitrogen content, soil organic matter. These soil physico-chemical properties also represent long-term effects of plants on soils. For example, soil C/N, total nitrogen content and soil organic matter depend on the quantity and quality of the plant litter inputs, and top-soil pH is also determined by interactions with the plant rhizosphere. We first further homogenized, dried and 2 mm-sieved the soil. A part of the soil was grounded to a particle size below 250 μm with an ultra-centrifugal grinder ZM 200 (Retsch ZM200) to determine total soil carbon and nitrogen contents using an elementary analyzer (Flash EA1112, Thermo Scientific). We gently crushed the rest of the samples, measured pH following the ISO 10390:2005 norm using a pH-meter (pH7110, inoLab) in a 1/5 solution (1 volumic part soil sample and 4 volumic parts distilled water), and quantified soil moisture content (after drying a subsample at 70°C for

h) and organic matter (SOM) content by Loss on Ignition (after 4 h incubation in a muffle furnace at 550°C).
From soil sub-samples B, we estimated exoenzymatic activity of seven exoenzymes involved in the degradation of C-rich substrates (EEC: α-Glucosidase (AG), β-1,4-Glucosidase (BG), β-D-Cellobiosidase (CB) and β-Xylosidase (XYL)), N-rich substrates (EEN: β-1,4-N-acetylglucosaminidase (NAG) and leucine aminopeptidase) and P-rich substrates (phosphatase (PHOS)) using standardized fluorimetric techniques (Bell et al. 2013). We homogenized soil during 1 min in a Waring blender in 200 ml of a sodium acetate buffer solution. Then soil slurry was added in technical duplicates to a 96-deep-well microplate with 200 µl of substrates specific to each enzyme at the saturation concentration. Exoenzymatic activities were assayed at pH 5 (mean pH across the sites, equivalent to meta-analysis from Sinsabaugh et al. 2008). For each soil sample duplicated standard curves were prepared by mixing 800 µl of soil slurry with 200 µl of 4-methylumbelliferone (MUB) or 7-amino-4-methylcoumarin (MUC) in 96-deep-well microplates with growing concentrations from 0 to 100 µM concentration. Plates were incubated at 25°C in the dark (3 h) on a rotary shaker (150 rpm) before centrifugation at 2900 g (3 min). Finally, fluorescence was measured on 250 µl of supernatant on a microplate reader (Varioscan Flash, Thermo Scientific) with excitation wavelength set to 365 nm and emission set to 450 nm. Exoenzymatic activities were expressed as nmol g soil −1 h −1 , after correcting for negative controls. We calculated the total potential extracellular enzymatic activity (total EEA) as the sum of the potential enzymatic activities of the seven exoenzymes as well as the ratio of the potential activities of exoenzymes targeting nitrogen vs phosphorus (EEN/EEP), the relative potential activity of exoenzymes targeting nitrogen (EEN/total EEA) and phosphorus (EEP/total EEA, see Table 2 for details), respectively.

Statistical analyses
We aimed to build a model in which species traits can influence the species response to environmental variables. We did so in two steps. First, we determined the shape of the species response to environmental variables using a flexible datadriven model (i.e. generalized mixed-effect additive model) and used a stepwise selection procedure to select the best set of variables. Second, with the selected environmental variables, and for each single trait, we built a final hierarchical species distribution model that includes species response to environment modulated by species traits.

Step 1 -selection of environmental variables and associated response shapes
Our list of environmental variables was initially as follows: climate (i.e. GDD, FDD, CWD, solar radiation), soil physico-chemical properties (i.e. soil C/N, pH, total nitrogen content, SOM) and microbial decomposition activity (i.e. total EEA, EEP/total EEA, EEN/total EEA, EEN/EEP). All variables were centered and scaled, while plant height was also log-transformed. 1) We removed variables that were too correlated using the 'findCorrelation' of the package 'caret' in R. If two variables are highly correlated, this function removes the variable with the highest mean absolute correlation. Average correlations are re-evaluated each time a variable is removed (Kuhn 2008). We also checked whether these choices made sense from an ecological point of view. SOM, EEN/EEP, CWD and solar radiation were dropped, in line with our expectations that they were moderately associated with plant species distributions (see Supplementary material Appendix 2 Fig.  A3 for details).
2) To determine the shape of species response to environmental variables, we then used a generalized additive mixedeffect model (GAMM) with species identity as random effects. Only the plant response to GDD was clearly unimodal and therefore we included GDD as a squared term in the subsequent analyses (Supplementary material Appendix 2 Fig. A4).
3) To select the best combination of environmental variables, we then performed a forward selection using Bayesian information criteria (BIC). Since the final hierarchical model is  (Allen et al. 1998, Vannier 2012) Sum of daily climatic water deficit: the sum of the negative daily climatic water balance, over the growing season, averaged over 1988-2018.
based on generalized linear mixed-effect models (GLMM), we here used a GLMM instead of a GAMM. The stepwise algorithm selection was constrained to stop for ΔBIC > 10 or after six abiotic variables were included. We used the final model resulting of this selection procedure (model 1) as the basis for the subsequent analyses (models 2-5 including traits).
Step 2  where the X c,i were the fixed environmental variable c ∈ {1…6} at site i selected by model 1, γ s was the random intercept for species identity, ε i was a normally distributed error term. The B c,s parameters were modeled as: where β c,s is the random slope for species identity, β c,0 can be interpreted as the coefficient for a hypothetical species with an average trait and β c,t described how trait T modulates the partial response of the species to the environmental variable X c . In this paper, trait T was averaged at the species level. The interaction between the mean trait T and the environmental variable X c provides information on how far the species is from its environmental optimum, and in that sense 'modulates' species probability of presence along environmental gradients. A positive trait-environment interaction (i.e. positive β c,t ) means that a species s with high trait value T s has a high probability of presence (high P s,i ) for a site i with high values of the considered environmental variable X c,i .
The four hierarchical GLMMs were fitted using a binomial response and a logit link. Species identity was included as a random effect, for which both slope and intercept were allowed to vary.
Parameters of all models were estimated by Laplace approximation of maximum likelihood using the 'Bobyga' optimiser of the 'lmer' function (package 'lme4' in R, R Core Team). We ensured model convergence by prohibiting the estimation of the correlations between random effects.

Model evaluation
To evaluate the hierarchical GLMM, we calculated the true skill statistic (TSS), which has the advantage to account both for the model sensitivity (i.e. proportion of observed presences predicted as presences) and specificity (i.e. the proportion of observed absences predicted as absences, Allouche et al. 2006). TSS can vary from −1 to 1, where +1 indicates perfect fit and values of zero or less indicate a performance no better or worse than random (Allouche et al. 2006). We also calculated TSS values per species. Since the TSS requires a threshold to transform species' probability of presence into binary presence-absence data, we selected the threshold that maximizes the TSS values.

Relative roles of climate, soil physico-chemical properties and microbial decomposition activity on plant distributions
The most parsimonious GLMM model identified with the stepwise selection procedure had a good predictive accuracy and included a single climatic variable together with three soil variables (TSS = 0.436, model 1; Fig. 2).
Overall, soil properties were key to predict plant distributions. Soil C/N, a proxy for organic matter decomposability, had the strongest impact on plant presence, followed by growing degree days (GDD 2 ), pH and total EEA. Across  Fig. 2).

Species with different functional traits respond differently to gradients of soil properties
In general, both model 1 (GLMM without traits) and the hierarchical trait-based GLMM models 2-5 (hierarchical trait-based GLMM) had good predictive accuracy for most species (average TSS = 0.46; Supplementary material Appendix 2 Fig. A5). Only three species had a TSS score lower than 0.3. Models accounting for functional traits (models 2, 3, 4 = 0.440, model 5 = 0.444) were slightly better than the model accounting for the environment only (model 1, TSS = 0.436).
Plant species showed different responses to the environment and notably a strong trade-off along the gradients of two soil characteristics (model 1, Fig. 3): Some species were more generally found under a combination of low soil C/N and high total EEA, while others were more common under a combination of high soil C/N and low total EEA. This tradeoff was not due to a sampling bias since the two gradients were independent in our study area (i.e. low correlation of soil C/N and EEA over all plots; Supplementary material Appendix 2 Fig. A3).
Species-specific responses could be partially related to their functional traits. Species with conservative traits (high leaf C/N, high LDMC) were favored where organic matter was hard to decompose (high soil C/N, Fig. 3, 4, models 2, 3). Furthermore, smaller species preferred low C/N and high total EEA sites (Fig. 3, 4, model 4). All plants benefited from microbial decomposition activity (high total EEA) but especially small species were favored on soils with high nutrient availability. Interestingly, no functional trait modulated species responses to pH or GDD (models 2, 3, 4, 5), but variability between species responses to the GDD gradient was also much lower compared to the variability along gradients of soil properties (Fig. 3, Supplementary material Appendix 2 Fig. A6-A).

Discussion
Understanding the mechanisms underlying plant distributions is a great challenge to anticipate the consequences of global changes on plant communities (Thuiller et al. 2008). Here, we develop a multi-species distribution model first to test how soil properties, especially microbial decomposition activities, shape plant species distributions in mountain environments, and second to account for how functional traits modulate species response to environmental gradients. We find that top-soil properties linked to nutrient recycling (microbial decomposition activity and organic matter decomposability) play a key role regarding plant distributions. Furthermore, we show that functional traits of plants related to resource acquisition strategies modulate plant response to soil properties. Exploitative species with Figure 2. Estimated standardized effect sizes of the environmental variables selected by the stepwise procedure using BIC (model 1). Red, green and blue points represent environmental variables related to climate, soil physico-chemical properties and exoenzymatic activities, respectively. acquisitive functional traits (e.g. low leaf C/N, low LDMC) are favored on soils with quick nutrient recycling, while conservative species (e.g. woody, with high leaf C/N and LDMC) dominate on nutrient-poor soils. In line with the 'fast-slow' plant economics spectrum (Reich 2014), we show that the spatial segregation of ecological strategies following top-soil properties helps to explain mountain plant distributions along elevational gradients.

Relative roles of climate, soil physico-chemical properties and microbial decomposition activity in explaining plant distributions in mountain ecosystems
Climate remained an important driver of plant distributions in our mountain system. Species probability of presence was maximized for intermediate values of GDD (model 1, Fig. 2). This was expected, given the importance of the length and intensity of the growing season for mountain plants (Körner 2003, Choler 2018. Interestingly, we show that top-soil properties play a key role for plant distributions, with an explanatory power comparable to that of climate. Increasing pH decreased the probability of presence of plant species (model 1, Fig. 2). In the literature, pH is consistently reported to be an important driver of plant distributions. Its influence is possibly due to the switch from calcareous to siliceous bedrock, for which plants are specialized or not (Lee 1998, Corlett andTomlinson 2020). Given that we measured pH on the first horizon, its influence might also be related to the interplay between plant litter, rhizosphere and soil physico-chemical properties, which affects microbial communities and nutrient supply to plants (Smiley 1974, Wang andTang 2018).
Our study demonstrates that soil properties associated with nutrient availability are tightly associated with plant species distributions. Overall, species presence probabilities increased in nutrient-rich conditions, such as those indicated by high microbial decomposition activity and organic matter of high quality (i.e. high total EEA and low soil C/N, Fig. 3,  4). Some studies found a significant effect of soil C/N on trees and woody species (e.g. positive on ericaceous as Vaccinium myrtillus, Coudun andGégout 2007, Walthert andMeier 2017), while others found weak effects on grassland plants (Dubuis et al. 2013, Buri et al. 2017. The varied effects of soil C/N were explained by the fact that measures of soil nutrient pools do not give information on the decomposition of the molecules involved (Dubuis et al. 2013, Buri et al. 2017. In mountain areas, soil C/N might be low although organic pools of nitrogen are tied up with soil organic matter as complex insoluble polymers that are not accessible to plants (Wardle and Bardgett 2002, Körner 2003, Buri et al. 2017. However, our results show that combining soil C/N as a measure of organic matter decomposability and total EEA as measure of the microbial decomposition activity improves the representation of nutrients available to plants both in forests and in grasslands and thus our capacity to predict plant distributions (Fig. 3).

Species with different functional traits respond differently to gradients of soil properties
Although all species respond positively to favorable environmental conditions, especially exploitative plants were advantaged in nutrient-rich places with mild climate. Under these favorable conditions, conservative species bear the costs of Figure 3. Species-specific partial response to total EEA as a function of species-specific partial response to soil C/N. Pie-charts represent the functional trait values of the species. The red and blue traits are expected to modulate the species response to soil properties and climate, respectively. The right-up box is a species (Vaccinium vitis-idaea) example with the legend for the functional traits. The left-bottom box shows the direction of increase for plant height, LDMC and leaf C/N. Full species names are given in Supplementary material Appendix 2 Table A1. 1557 being excluded by plants with a more exploitative strategy (mostly herbaceous species, Fig. 3, Supplementary material Appendix 2 Fig. A6-7). Surprisingly, while we expected traits related to light competition (i.e. tall plants with high SLA) to confer a particular advantage to plants under mild environmental conditions, we found no significant modulating effect of SLA for any of the environmental axes, and height conferred an advantage rather on nutrient-poorer soils. This is in opposition to the literature on broader grassland types (Pollock et al. 2012, Carboni et al. 2018), but can be due to the fact that alpine plants are overall selected for rather small sizes and a small range of SLA to resist to cold temperatures (Körner 2003). The advantage for tall plants on nutrient-poorer soils is due to the fact that within our dataset those plants were mostly woody species with a slower growth rate (Juniperus communis, Vaccinium myrtillus, Picea abies, Vaccinium vitis-idaea, Fig. 3, Supplementary material Appendix 2 Fig. A7).
In our alpine study, the variability among the responses of species along the gradient of soil properties was bigger than along the climatic gradient (Supplementary material Appendix 2 Fig. A6), which could suggest that spatial segregation of plants with different strategies was mostly driven by soil properties (Reich 2014). Although the probability of presence of most plants increased in places with high soil organic matter decomposability (low soil C/N), the probability of presence of some species increased on nutrientpoor soils (e.g. woody species, Fig. 3, Supplementary material Appendix 2 Fig. A6-7), reflecting a differentiated strategy with selection for conservative traits, i.e. higher leaf C/N and LDMC (Fig. 4, Supplementary material Appendix 2 Fig.  A7). Conservative leaf traits are necessary on nutrient-poorer soils with high soil C/N and give an advantage reinforced by plant-soil feedbacks (Aerts 1995, Reich 2014. Leaves with higher C/N contribute to the low decomposability level of the soil organic matter, which slows down decomposition, and reinforces high soil C/N (Grigulis et al. 2013, Legay et al. 2016). This confirms the hypothesis of Aerts (1999), which states that in nutrient-poor habitats, selection is based more on the ability of plants to limit nutrient losses than on their ability to access nutrients quickly.
Although our results suggest that exploitative plants are often excluded as a result of strong nutrient limitations of mountain soils, nutrients such as nitrogen are not only supplied to plants from the decomposition of soil organic matter. With the continuous increase of atmospheric nitrogen deposition since the industrial age (Galloway et al. 2008), nitrogen limitation might become less problematic in the future. Associated with milder temperatures, increased atmospheric nitrogen depositions in mountain areas could result in a shift towards communities dominated by exploitative plants and the competitive exclusion of conservative ones (Hautier et al. 2009, Boutin et al. 2017. As plants with exploitative traits produce a litter easier to decompose, these community shifts might drive accelerated nutrient cycling rates in alpine systems, with unknown consequences on the carbon and nitrogen cycles at larger scales. Finally, we used the mean of traits measured at different places, mostly in the French Alps, which represents the mean habitat suitability of the modeled species (Albert et al. 2012). However, within species, plants also respond to gradients of climatic and soil properties through modifications of the values of response traits by phenotypic plasticity or local adaptation (Lavorel and Garnier 2002). It has been shown that these variations around the mean traits influence the relationships of diversity measures with environmental gradients (Albert et al. 2012), and that intra-specific trait variability influences species distributions at small scales in alpine systems (Chalmandrier et al. 2017). We might thus expect that intra-specific trait distributions would influence our traitbased distribution models. For example, it might reinforce the importance of traits-environment interactions or improve the explanatory power of trait-based species distribution models by providing information at finer scale for each species. In that sense, accounting for intraspecific trait variability is an interesting perspective of trait-based species distribution models.
To conclude, by bridging soil biogeochemistry and plant biogeography we provide novel insights on how soil properties linked to nutrient availability shape mountain plant distributions in addition to climate. We show that incorporating microbial decomposition activity improves plant distribution models, probably because it allows estimating plant nutrient supply. Finally, our results support our hypothesis that plant functional traits associated with conservative-exploitative strategies modulate distributions along gradients of top-soil properties associated with nutrient availability. Future prospects include a better understanding of the influence of the interaction between climate change and atmospheric nitrogen deposition on plants and the inclusion of intraspecies trait variability in large-scale distribution models.

Data availability statement
Most of the raw data used in the manuscript are directly available on the ORCHAMP website (<https://orchamp.osug.fr/ home>). The data supporting the results are available from the Dryad Digital Repository: <https://doi.org/10.5061/ dryad.2z34tmpjg> (Martinez-Almoyna et al. 2020).