The role of climate, foliar stoichiometry and plant diversity on ecosystem carbon balance

Global change is affecting terrestrial carbon (C) balances. The effect of climate on ecosystem C balance has been largely explored, but the roles of other concurrently changing factors, such as diversity and nutrient availability, remain elusive. We used eddy‐covariance C‐flux measurements from 62 ecosystems from which we compiled information on climate, ecosystem type, stand age, species abundance and foliar concentrations of N and P of the main species, to assess their importance in the ecosystem C balance. Climate and productivity were the main determinants of ecosystem C balance and its stability. In P‐rich sites, increasing N was related to increased gross primary production and respiration and vice versa, but reduced net C uptake. Our analyses did not provide evidence for a strong relation between ecosystem diversity and their productivity and stability. Nonetheless, these results suggest that nutrient imbalances and, potentially, diversity loss may alter future global C balance.

The large repercussions of all these impacts on the Earth system and, therefore, on human societies Peñuelas et al., 2013) have triggered intensive research into the effects of global change on ecosystems. One of the most important features of ecosystems, given the current climate emergency, is their carbon (C) balance, this is: gross primary production (GPP), ecosystem respiration (Re) and net ecosystem production (NEP = GPP − Re).
The ecosystem-level C-use efficiency (CUEe), defined as the NEP:GPP ratio, is another important feature to understand how C is used in ecosystems (Fernández-Martínez, . While the role of climate and the alteration of the N cycle (i.e. N deposition) on ecosystem C balance have been intensively investigated (de Vries et al., 2014;Fernández-Martínez et al., 2017;Janssens et al., 2010;van Dijk et al., 2005), we barely have any evidence regarding the role of diversity and elemental composition of ecosystems on C balance. Filling in this knowledge gap is crucial given that different facets of the global change are currently reducing biodiversity and, therefore, directly (altering species elemental concentrations) and indirectly (changing species pools) shifting the elemental composition of ecosystems. On the other hand, both diversity and elemental composition of organisms have been demonstrated to play an important role on ecosystem functioning Sterner & Elser, 2002), thus entailing ecosystem changes if species are lost or stoichiometric imbalances occur.
It is now widely accepted that ecosystems with higher diversity also present higher productivity (often biomass production) and lower temporal variability Musavi et al., 2017;Oehri et al., 2017) in both, aquatic and terrestrial environments, across the globe (Liang et al., 2016). However, despite the fact that decomposition rates have also been suggested to increase with the diversity of decomposers and plant litter (Handa et al., 2014), it is not well understood how diversity relates to respiration of entire ecosystems, and, hence, to net C uptake. It is, therefore, imperative to understand the role of diversity on the ecosystem net C balance. On the other hand, there is evidence pointing towards an important role of nutrients on biomass production (Vicca et al., 2012) and ecosystem net C balance (Fernández-Martínez, . These findings suggest higher productivity and net C uptake in nutrient-rich ecosystems, where concentrations of foliar N and P of primary producers are also often higher (Firn et al., 2019;Sardans et al., 2016;Van Sundert et al., 2018), potentially leading to enhanced photosynthesis (Domingues et al., 2010;Field et al., 1983;Wright et al., 2004) and lower C losses by the emission of volatile organic compounds (Fernández-Martínez, Llusià, et al., 2018). However, the direct link between community-weighted elemental composition and C balance still remains elusive.
Hence, the goal of this study was to assess the importance, for the first time, of the effects of climate, endogenous factors (i.e. ecosystem type, productivity and stand age), diversity and elemental stoichiometry on ecosystem C balances. We here hypothesized that, once the effect of climate and endogenous factors are accounted for, larger and less interannually variable C fluxes (GPP, Re, NEP, their interannual variabilites: GPP IAV , Re IAV and NEP IAV , and CUEe) will be found in sites with higher diversity (species and phylogenetic diversity [PD]). Similarly, we hypothesized that larger and less variable C fluxes will be associated with higher foliar concentrations of N and P at the community level. To test these hypotheses, we used 62 eddy-covariance sites included in the FLUXNET 2015 database, comprising forest, shrubland, savanna and grassland ecosystems ( Figure S1; Table S1). For all 62 sites, we compiled information regarding their species abundance and their community-weighted foliar N and P concentrations, either measured or inferred for the species present.

| Carbon fluxes and climate data
We downloaded monthly weather and C flux (GPP, Re and NEP) data from eddy-covariance towers from the FLUXNET 2015 data set (Tier 2; http://fluxn et.fluxd ata.org/data/downl oad-data; Pastorello et al., 2020). From the total of 93 sites for which we could gather species abundance data, croplands excluded (see Section 2.1.2), only 62 sites met the conditions for data quality (NEE_VUT_REF_QC ≥ 0.5, following Bond-Lamberty et al., 2018), had data for at least 60 months (5 years) and did not experience strong management or disturbances ( Figure S1; Table S1). The 62 sites that were selected included boreal, temperate and Mediterranean biomes and 10 different vegetation types according the IGBP-DISCover classification (Bonan et al., 2002): evergreen needle-leaved forests (19 sites), evergreen broadleaved forests (3), deciduous broadleaved forests (13), mixed forests (6), savannas (2), woody savannas (3), open shrublands (3), closed shrublands (1), wetlands (2) and grasslands (10). We grouped these vegetation types into three categories in order to facilitate statistical analyses, being forests (evergreen and deciduous, needle-leaved and broadleaved, and mixed forests), shrublands and savannas (savannas and woody savannas, open and closed shrublands and wetlands) and grasslands. We used all monthly values of GPP and Re estimated with the daytime partitioning method. We then calculated annual values of the C fluxes as the sum of GPP, Re and NEP average monthly values. Additionally, we calculated average CUEe (Fernández-Martínez,  as the ratio of mean NEP and GPP per site. We also calculated mean annual temperature (MAT) and precipitation (MAP) per site to control for climate effects in our statistical models.
Interannual variability of C fluxes per site was estimated using the proportional variability index (PV; Fernández-Martínez, Heath, 2006), as the average of the interannual PV for each month (e.g. Interannual variability of climate variables (MAT IAV and MAP IAV ) was estimated using the same method. We collected stand age of forests from ancillary data and a global forest database (Campioli et al., 2015;Luyssaert, Inglima, et al., 2007).

| Species composition data
Information on species abundance per site was extracted from the study by Musavi et al. (2016Musavi et al. ( , 2017, ancillary data for each site, published literature containing information about the selected sites and information provided after contacting the PI's of the sites. To avoid underestimation of species richness in sites where inventories were not very accurate, we only used those species that accounted for up to 95% of the abundance of the site, similar to the study by Musavi et al., (2017). Abundances were estimated as the fraction of leaf area occupied by each species in the sites for which this information was available. The relative density of individuals or direct reports of relative abundance were used in the rest of sites. When abundances per species were not reported, we assigned equal abundances to each species within their strata (e.g. overstorey or understorey) based on the average percentage of over-and understorey of each vegetation type (e.g. evergreen needle-leaved forests have 85% overstorey and 15% understorey, so if a site would have two species in the overstorey, they both would have 42.5% relative abundance). When only vague indications about species abundance was reported, we assigned differential abundances accordingly. Species per site and their abundances can be found in Data S1 (online only-spec_abund_full.csv). Vascular plants, ferns and bryophytes were used in our analyses given that in some sites non-vascular plants could account for a considerable proportion of the species present (>10%).

| Elemental composition data
We collected all available information about foliar N and P concentration per species at each of the selected sites, including data from previous compilations (Musavi et al., 2016). We searched for that information in the main literature of the sites and by directly requesting this information to the PI's of each site. Overall, we gathered data of N and P foliar concentrations for 21 and 18%, respectively, of the relative abundance of species per site in our database (see specs_abund_full. csv file in Data S1). For the rest of the species for which we did not find any information on foliar N and P concentrations, we attributed values of foliar N and P based on average values of the species. To that aim, we used N and P data per species from the study by Tian et al. (2019). For the purpose of this study, we calculated median concentrations of N and P per species. We used the Rphylopars (Goolsby et al., 2017) R package to estimate the elemental composition of N and P from phylogenetic information (Penone et al., 2014) for those species for which we did not find literature data on their foliar N and P concentrations. Imputation of N and P for missing species was performed based on the phylogenetic tree of all species present in the selected sites and in our database of elemental composition. We used the plant megaphylogeny provided by Qian and Jin (2016). Species not found in the base phylogeny were included into our phylogenetic tree using the V.PhyloMaker (Jin & Qian, 2019) R package. Correlation between imputed and observed values (in the selected sites) was very good for both N and P (R = .91 and .83; 100 and 94 replicates, p < .001 respectively) indicating low uncertainty in our estimates of foliar N and P concentrations. All species names were checked and corrected using The Plant List database, using the R package Taxonstand (Cayuela et al., 2016). We finally calculated the community-weighted mean N and P concentrations based on the species abundance per site and their corresponding N and P concentrations.

| Diversity indices
We first calculated three different metrics of PD per site using the PhyloMeasures package in R (Tsirogiannis & Sandel, 2016. We estimated PD, mean pairwise distance (MPD) and mean nearest taxon distance (MNTD) per site; we standardized (mean = 0, SD = 1) all three variables and calculated their average as a measure of PD.
PD measures phylogenetic diversity as the total branch length in the minimum spanning subtree for a given set of species (Faith, 1992), MPD measures the MPD along the tree between all species in a given set, and MNTD measures the mean distance from each taxon to its nearest neighbour in the set of selected species (Tsirogiannis & Sandel, 2016;Webb, 2000). We also estimated species diversity as the average of the standardized values of species richness per site, Shannon, Simpson, inverse Simpson, unbiased Simpson and evenness indices (Daly et al., 2018;Hill, 1973;Hurlbert, 1971). All these indices, except species richness per site, were calculated using the vegan (Oksanen et al., 2018) package in R.

| Statistical models
We used a model averaging approach based on linear models to investigate the effects of endogenous factors, climate, diversity and foliar nutrient concentrations on annual C fluxes (GPP, Re and NEP), their interannual variability (GPP IAV , Re IAV and NEP IAV ) and CUEe.
Model averaging is a statistical technique based on multimodel inference that computes an average model from the estimates of a set of models including different predictors and weighting their relative importance according to the difference of the second-order AIC (AICc) between each model and the best model . For annual fluxes, we fitted models containing ecosystem type (forest, shrublands-savannas or grasslands), MAT, MAP and their interaction (MAT × MAP), species diversity and PD, community-weighted foliar N and P and their interaction. For Re and NEP, we also included GPP as a predictor to control for the productivity of the site as an endogenous factor. Models explaining variability in CUEe included the same predictors but without including C fluxes. For interannual variability, we repeated the same models as for C fluxes described above but including also the GPP, Re and NEP as predictors. For GPP IAV , only GPP was included, for Re IAV we included GPP and Re, and for NEP IAV we included GPP, Re and NEP. Diversity and foliar nutrient variables were log-transformed to account for non-linear relationships.
We then estimated all possible models nested within the saturated (see above) models using the dredge function in MuMin (Barton, 2018;) R package. Models including both MAT and MAT IAV were discarded from the set of models because of the high collinearity between these two predictors (variance inflation factor > 10). We selected those models with ΔAICc < 4 to estimate the conditional average model and calculate the relative importance of each predictor.
For each of these models we also calculated the variance explained by each of the predictors with the 'lmg' metric in the R package 'relaimpo' (Grömping, 2006). We then calculated the average variance explained by all predictors weighting by the relative importance of each of the models with ΔAICc < 4. We finally used models with ΔAICc < 4 to visualize the relationships between response and predictor variables using partial residual plots by means of the visreg (Breheny & Burchett, 2015) package in R.
We repeated all statistical analyses using (a) all sites (N = 62), (b) only forests (N = 41) and (c) only savannas, shrublands and grasslands (N = 21). In models using only forests, we also included stand age (log-transformed) and its interaction with C fluxes as a predictor (see models described above). All analyses were performed in R statistical software v. 3.6.2 (R Core Team, 2019).

| Main controls of C fluxes and their temporal variability
Our results indicated that climate was the main factor explaining variability in GPP across sites (Figure 1). Ecosystems with higher F I G U R E 1 Contribution of climate, endogenous factors, community-weighted foliar nutrients and diversity on terrestrial C balance. (a) Relative (R, within each model) and absolute (A) variance explained by four different type of variables (climate, endogenous [site type, productivity], nutrients and diversity) for models with ΔAICc < 4. (b) Relative importance of the most relevant variables selected by the models (at least one relationship with relative importance >0.5) and the direction of their effect. Empty and filled dots indicate p < .1 and p < .05 in the average model respectively (see Data S1 -Models). Diversity Sp , species diversity; GPP, gross primary production; IAV, interannual variability; MAP, mean annual precipitation; MAT, mean annual temperature; NEP, net primary production; Re, ecosystem respiration MAP and MAT and less variable climate (MAT IAV and MAP IAV ) had generally higher GPP. The main climatic effect on GPP propagated into Re and NEP, for which the direct climate effect was very low (and potentially negative for Re) once site productivity (GPP) was taken into account. This effect was similar to that occurred in interannual variability in C fluxes, which was strongly negatively related to the annual C fluxes (Figure 2). Across all ecosystem types, warmer sites presented higher Re IAV and NEP IAV , while across forests MAP was positively related to higher Re IAV and NEP IAV ( Figure S2). Instead, MAT IAV often presented a negative relationship with variability in C fluxes, especially in savannas and grasslands ( Figure S3).
Community-weighted foliar N and P explained a relatively large proportion of the spatial variability in GPP and NEP annual sums and their interannual variability. We found a statistically significant interaction between community-weighted foliar N and P concentrations explaining variability in GPP, Re and NEP (Figure 3). Our results showed that increasing foliar N translated into increases in GPP only in sites with higher concentrations of foliar P. On the other hand, increasing foliar P had a positive effect on GPP in sites with high foliar N concentration, while the opposite relationship was found across sites with low foliar N. Similarly, Re only increased with increasing foliar N and P when their counterpart nutrient was also high. As a result of these interactions between nutrients, GPP and Re, we found a negative interaction between foliar N and P explaining variability in NEP. Our analyses indicate that increasing foliar N and P is only slightly related to increased NEP in sites with low foliar N and P concentrations. These interactions, however, were either lost or had a very low importance in the average model when forest and non-forest ecosystems were analysed separately (Figure 1; Figures S2 and S3). We also found that variability in C fluxes was positively related to the interaction between foliar N and P, pointing towards higher variability when both nutrients were high. Across non-forest ecosystems, those findings remained only for foliar P.
Instead, the role of diversity was very small on annual C fluxes and their interannual variability when comparing all ecosystems ( Figure 1). Across all ecosystem types and across forests, sites with higher species diversity tended to present higher GPP and Re but lower NEP and higher GPP IAV . Across forests, instead, diversity did not seem to affect C fluxes nor their variability ( Figure S2). PD, however, was found to be negatively related to temporal variability in NEP across non-forest ecosystems ( Figure S3). We additionally found that stand age was an important predictor of GPP across forests. Older stands are more productive (GPP) and exhibit more temporarily stable net C uptake (lower NEP IAV ; Figures S2 and S4).

F I G U R E 2 Partial residuals plots showing the predictors of interannual variability in annual fluxes. Y-axes in (a) and (d) show GPP IAV , (b)
and (e) show Re IAV , and (c) and (f) show NEP IAV . GPP, gross primary production; IAV, interannual variability; MAP, mean annual precipitation; MAT, mean annual temperature; NEP, net primary production; PV, proportional variability; Re, ecosystem respiration F I G U R E 3 Partial residual plots showing the effect of the interaction between foliar N and P concentrations on C fluxes. Models can be found in Data S1. Nutrient levels that are compared correspond to percentiles 20 and 80 in the nutrient concentration ranges. Panels (a) and (b), (c) and (d) and (e) and (f) show, respectively, the effect of the interaction between foliar N and P on GPP, Re and NEP. DW, dry weight; GPP, gross primary production; NEP, net primary production; Re, ecosystem respiration

| Controls of carbon-use efficiency at the ecosystem level
Across sites, CUEe, defined as the NEP:GPP ratio, was mainly explained by climate and the type of ecosystem (Figure 4). On average, CUEe was higher in forests than in grasslands, but neither were statistically different from shrublands and savannas. Higher CUEe was consistently found in sites with higher MAP but low MAT and vice versa, especially across forests. Species diversity was negatively related to CUEe in non-forest ecosystems. Foliar N was mainly positively correlated with higher CUEe, but showed a strong negative correlation across non-forest ecosystems.

| Mechanisms behind the effects of foliar elemental composition on ecosystem C balance
As hypothesized, and consistent with findings in the literature (Domingues et al., 2010;Field et al., 1983;Fleischer et al., 2013;Reich et al., 2009;Wright et al., 2004), our results show that high concentrations of foliar N and P are associated with larger C fluxes (both GPP and Re), but only when both elements were found in high concentrations ( Figure 3). This result may emerge because of the inability of plants to increase their growth when one nutrient is limiting (Peng et al., 2019). Growth is, therefore, often boosted by a synergistic effect of both N and P. However, the effect is stronger on Re than GPP, and as a result net C uptake (NEP) presented the opposite behaviour. This is unexpected considering the previously reported positive effect of nutrients on C sinks (Fernández-Martínez, . Differences in these results, however, may emerge from the different approaches used: foliar nutrients in the present study and site fertility (mainly soil) in the previous one (Fernández-Martínez, . Adaptations to take up and retain nutrients shown by different plant species often mask the relationships between site fertility and foliar elemental concentrations (Van Sundert et al., 2019). In this sense, we did not find a direct statistical relationship between foliar N and P and site fertility used in previous studies    productivity], nutrients and diversity) for models with ΔAICc < 4. (b) Relative importance of the most relevant variables selected by the models (at least one relationship with relative importance >0.5) and the direction of their effect. Empty and filled dots indicate p < .1 and p < .05 in the average model respectively (see Data S1: Models). (c) Partial residual plots showing differences in CUEe among ecosystem types. GPP, gross primary production; IAV, interannual variability; MAP, mean annual precipitation; MAT, mean annual temperature; NEP, net primary production; Re, ecosystem respiration; SSG, savannas, shrublands and grasslands communities and, therefore, organic matter of better quality (Enríquez et al., 1993;Kozovits et al., 2007). The hypothesis that nutrient imbalances (C:N:P) and retarded decomposition may increase C sink capacity has actually also been suggested as one of the potential explanations for the recent increases in NEP in Asia  given the higher positive effect of the declining sulphur (S) deposition on Re than GPP found in forests (Fernández-Martínez et al., 2017). However, our results indicated that foliar nutrients played a very limited role explaining CUEe, presenting a negative relationship between foliar N and CUEe across shrublands, savannas and grasslands once MAP and species diversity were taken into account (Figure 4). These results are counter-intuitive considering the positive effect that foliar N generally has on GPP (Fleischer et al., 2013) for which they could only be explained because of a proportionally higher increase in Re than in GPP due to higher N availability in the litter (Norby et al., 2001), albeit that was not found in our analyses ( Figure S3). Across forests, however, we found weak evidence suggesting that increasing foliar N translates into higher CUEe. Further research is therefore needed to corroborate these findings and identify the mechanisms driving them.
Despite the fact that temporal variability in C fluxes are negatively related to C fluxes themselves ( Figure 2) and the latter tend to be larger when both foliar N and P are high, our initial hypothesis stating that higher foliar nutrients would be related to more stable C fluxes was not confirmed by our results. Contrary to our expectations, when controlling for the magnitude of the C flux, we found evidence for a positive effect of foliar N and P on interannual variability of C fluxes (Figure 1; Figures S2 and S3). These positive relationships could emerge as a result of a stronger link between C fluxes and environmental variability in sites with higher foliar nutrient concentrations. Fluctuating environmental conditions (i.e. temperature, water availability) could then be one of the few metabolic constraints, allowing larger peaks of photosynthesis and decomposition when nutrient requirements are satisfied, but we have no evidence in support of this mechanism.

| The influence of plant diversity on ecosystem C balance
Our results provided weak evidence that higher plant species diversity is associated with higher GPP and Re (Figure 1). Thus, our initial hypothesis stating that larger and less variable C fluxes are found at sites with larger diversity finds insufficient support from this data set. The small effect of diversity on NEP and CUEe was even found to be negative (Figure 4; Figures S2 and S3). This suggests that any positive effect of diversity on GPP is being overcompensated by an increase in Re. Nonetheless, the positive relationship between species diversity and C fluxes (GPP and Re) may emerge because of selection and complementarity effects . The selection effect refers to the higher probability of diverse communities to include one or more highly productive species. The complementarity effects contribute to the increased productivity and lower variability over time by promoting temporal asynchrony in the productivity of the species of a community (Ammer, 2019;Cardinale et al., 2012;Gessner et al., 2010;Loreau & de Mazancourt, 2013). Larger diversity of plant species may then entail a larger diversity of decomposers, which would hence contribute to a positive effect on decomposition of organic matter and, therefore, respiration. Instead, the negative relationship between PD and temporal variability in NEP across non-forest ecosystems was stronger and in agreement with our initial hypotheses. These results may also indicate that complementary effects and niche partitioning may be more relevant in non-forest ecosystems because the usual inference is that the larger the PD at a site, the more likely it is that the species are functionally different. Temporal asynchrony is then easier to achieve if species are phylogenetically, and thereby likely, functionally distant, because they are more likely to respond differently to environmental variability among each other.

| Climate and stand age effects on ecosystem C balance
Climate played a much more important role in controlling annual C fluxes and their temporal variability than either elemental composition or diversity (Figures 1 and 4; Figures S2 and S3). We found a general trend in increasing variability in C fluxes in warmer climates and a very interesting negative effect of climate variability on GPP and Re. The latter relationship may be indicative of lower performance under periods of stress due to too high or too low temperatures or precipitation. In forest ecosystems, stand age also presented a relevant role on the C balance ( Figure S4), supporting previous research (Carey et al., 2001;Luyssaert et al., 2008;Musavi et al., 2017). We found that GPP increases and C sinks become more stable as forests age, suggesting that an increase in ecosystem complexity may dampen the effects of environmental fluctuations by making ecosystems more resistant.

| Implications for the global carbon balance
Our results provide the first assessment of the combined role of climate, endogenous factors, diversity and ecosystem foliar stoichiometry on terrestrial ecosystem C balances. We found an important effect of climate and endogenous factors on ecosystem C balance, similar to that was found in previous studies (Fernández-Martínez, Luyssaert, Inglima, et al., 2007). Furthermore, our analyses also found a negative effect of climate interannual variability on GPP and Re (Figures 1 and 2), a facet of climate barely explored in C balance studies. Variability in precipitation was also likely to increase variability in C fluxes. These findings suggest that increasing climate variability, as predicted for some regions, will alter C fluxes and potentially increase their temporal variability. The important effect of climate variability on ecosystem C balance highlights the need to take into account the temporal variability of climate in addition to average values, when exploring the impacts on biological processes.
We found that foliar stoichiometry and especially plant species diversity explain a relatively low amount of the spatial variability in C fluxes and their variability over time, compared with the variance explained by climate. The effect of diversity found here seems to be lower than that shown in previous studies on small and artificial plots (Cardinale et al., 2004;Gross et al., 2014) or using forest inventories (Liang et al., 2016). The effect of the ecosystem foliar stoichiometry may also seem to be slightly lower than those in previous studies using proxy metrics of soil fertility Vicca et al., 2012). There may be several explanations for the relatively low effects of diversity and foliar nutrients that we found on ecosystem C balance.
First, C fluxes estimated by eddy-covariance towers inherently contain a rather large uncertainty themselves, estimated to be in the order of 30% of the total annual flux (Luyssaert, Inglima, et al., 2007). This uncertainty may disproportionally mask relatively small effects, operating at small scales and short gradients, such as those of diversity or elemental composition, compared to strong effects at much larger scales, such as climate. Second, our assessment of species abundance was based on literature research and contact with the principal investigators of the sites because data from standardized protocols do not exist yet. Therefore, species abundance used here is far, in terms of quality, from the information that standardized botanical surveys would provide. Openly available standardized botanical surveys at each site would definitely help to better elucidate the role of diversity on ecosystem C balances and flux networks should promote their broad implementation. Third, the elemental composition of the species within a site was reconstructed based on species identity for the majority of the sites. For several species, we had to use elemental composition data inferred from phylogeny, because foliar chemistry for the target species was not present in our database. The fact that differences in elemental composition is higher among species than within species , however, supports the validity of our approach. Additionally, all these above-mentioned uncertainties in our data are systematic across sites, which means that it would be very unlikely that any of our findings emerge as a result of them. Despite these limitations, we found statistically significant differences in ecosystem fluxes and their interannual variability in response to differences in diversity and ecosystem elemental composition. Hence, we hypothesize that the effect of biodiversity and foliar elemental composition on ecosystem C balance may be stronger than the one that we report here.
Our results, point out that global change is affecting the global C cycle, not only through climate change, but also through shifts in foliar stoichiometry of autotrophs, such as increasing N and decreasing P concentrations (Penuelas et al., 2020;Peñuelas et al., 2013).
Whether diversity loss will have a strong impact on global C balance still deserves further attention.

ACK N OWLED G EM ENTS
This work used eddy-covariance data acquired and shared by the

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available