Plant community functional composition and diversity drive fine-scale variability in carbon cycling in the tundra

1. The functional composition and diversity of plant communities are globally applicable predictors of ecosystem functioning. However, exactly how traits influence carbon cycling is yet unclear, as are the implications in a warming world. 2. To study how traits affect carbon cycling in the tundra environment, we built a hierarchical model that included abiotic conditions (summer air and winter soil temperatures, and soil resources), plant community functional composition and diversity (plant size and leaf economics), and carbon cycling (above-ground and soil organic carbon stocks, and photosynthetic and respiratory fluxes). We also simulated warming effects on peak-season ecosystem CO ​ 2 ​ budget. 3. Plant size was the strongest predictor of all carbon cycling variables except soil carbon stocks. Communities of larger plants were associated with larger CO ​ 2 ​ fluxes and above-ground carbon stocks. Communities with fast leaf economics had higher rates of photosynthesis and soil respiration, but lower above-ground biomass. 4. Diversities on axes of size and leaf economics affected ecosystem functions differently. Leaf economic diversity increased CO ​ 2 ​ fluxes and soil organic carbon stocks, while size diversity increased the above-ground carbon stock. The contributions of functional diversity metrics to ecosystem functioning were about as important as average leaf economic traits. 5. Simulations suggested that warmer summer air temperatures increase plant size, while warmer winter soil temperatures increase plant size and accelerate leaf economics. All these changes would enhance CO ​ 2 ​ uptake during peak growing season. Synthesis: We show how traits mediate the link between abiotic conditions and carbon cycling. Community composition and diversity on the two axes of the global spectrum of plant form and function have clear and separate effects on ecosystem functioning. Warmer temperatures increase plant size and accelerate leaf economics, which leads to faster net assimilation of carbon during peak growing season. More research on soil carbon stocks is needed.


Introduction
The Arctic is warming two to four times faster than the planet on average (Holland and Landrum 2015) , leading to changes in plant growth (Myers-Smith et al. 2015) , height (Bjorkman, Myers-Smith, Elmendorf, Normand, Rüger, et al. 2018) , and species distributions (Steinbauer et al. 2018) .It is unclear how such changes in vegetation will influence the carbon balance of these ecosystems.On one hand, the expansion of larger plants might increase net carbon uptake (Cahoon et al. 2012;Sørensen et al. 2019) .On the other hand, vegetation-driven changes in soil conditions might accelerate decomposition and, in turn, increase carbon losses to the atmosphere (Parker, Subke, and Wookey 2015) .Arctic soils contain more than half of the global soil organic carbon stock (Hugelius et al. 2014) , thus changes in tundra carbon cycling have the potential to influence atmospheric CO 2 concentrations at the global scale.
Despite the importance of vegetation for climate change-related feedbacks in the carbon cycle, most carbon cycling studies have described the functional properties of vegetation with coarse proxies, such as plant functional groups or vegetation classification (Virkkala et al. 2018) .Recent research has shown that classifications might ignore large variation in vegetation properties that are relevant for ecosystem functioning (Thomas et al. 2019;Cadotte, Carscadden, and Mirotchnick 2011) .Furthermore, the use of classification methods makes it impractical to compare the results of independent studies if their classifications differ.Understanding the links between plant community composition and carbon cycling in the Arctic remains challenging due to the lack of studies that use comparable and informative methods to characterize tundra vegetation.
Using globally applicable quantitative plant functional traits provides a solution to this problem.Functional traits, which represent "morpho-physio-phenological traits which impact fitness indirectly via their effects on growth, reproduction and survival" (Violle et al. 2007) , are a means to link community composition to environmental conditions and ecosystem functioning (Lavorel and Garnier 2002) .Above-ground traits of plant species and the functional composition of plant communities have been found to vary primarily along two axes (Díaz et al. 2016;Bruelheide et al. 2018) .The first axis describes the trade-offs between resource acquisitive and conservative strategies, and has often been called the leaf economics spectrum (Wright et al. 2004) which is measured with traits such as specific leaf area (SLA), leaf nitrogen content, and leaf dry matter content (LDMC) (Bruelheide et al. 2018) .The other axis is related to plant size, characterizing the trade-offs between investments to light competition, and photosynthesis and reproduction, and can be measured with, for example, plant vegetative height or stem specific density (Díaz et al. 2016) .These trait axes explain differences between species' vital rates that define population processes (Adler et al. 2014) , and ecosystem functioning from local to global scales (Diaz et al. 2004;Michaletz et al. 2014) .
Carbon cycling is a set of ecosystem functions that comprises fluxes (flows) and stocks (quantities), both above-and below-ground (Table 1).These functions are most often linked to traits in the context of the mass-ratio hypothesis (Grime 1998) , which states that ecosystem processes are determined by the trait values of dominant species in the community.Today, this hypothesis is most often operationalized using abundance-weighted average trait values (Garnier et al. 2004) .In the tundra, photosynthesis and ecosystem respiration have been shown to correlate positively with fast leaf economics, indicated by, for example, high average SLA (Sørensen et al. 2019) .Plant size is often described using biomass or its proxies, such as total plant cover or leaf area index (LAI) (Oberbauer et al. 2007;Street et al. 2007;Shaver et al. 2007;Marushchak et al. 2013) .Studies have shown that larger biomass leads to higher photosynthesis, net carbon uptake, and soil organic carbon stocks (Lafleur and Humphreys 2018;Gagnon, Domine, and Boudreau 2019) .Plant vegetative height is more rarely used to explain carbon cycling in the tundra, although it has several advantages over biomass proxies.Vegetative height is directly relevant for species niches via effects on e.g.light competition, as it describes how much of an individual's leaf area can actually intercept light (Yachi and Loreau 2007) , and is thus more suitable for modelling the effects of environmental selection on community composition.Moreover, species-specific plant height can be used to estimate functional diversity, unlike measurements of total biomass or LAI.This last point is important because it is well known that in addition to average traits, the diversity of communities also influences ecosystem processes (Tilman 1997) .Both experimental and observational studies strongly agree that ecosystem functioning generally increases along with the diversity of local species assemblages (Tilman, Isbell, and Cowles 2014;Duffy, Godwin, and Cardinale 2017) .Using functional traits allows partitioning the diversity effect to contributions from complementarity along different niche axes (Cadotte, Carscadden, and Mirotchnick 2011) .We might expect that the diversity of leaf economic strategies is linked to different ecosystem functions compared to plant size diversity.For example, height diversity correlated positively with increased biomass in Canadian experimental fields (Wilsey and Potvin 2000) , whereas leaf nitrogen content diversity increased net ecosystem exchange and ecosystem respiration in experimental grassland mesocosms (Milcu et al. 2014) .To our knowledge, there is a lack of studies investigating functional diversity effects on ecosystem functioning in the tundra.
Plant functional composition and diversity as well as carbon cycling are driven by a range of abiotic environmental conditions, such as vital soil resources (i.e.soil moisture and nutrients) and microclimatic conditions during the summer and winter .Soil resource availability has been established as an important control of not only leaf economics and plant size (Garnier et al. 2004;Spasojevic and Suding 2012;Pérez-Ramos et al. 2012) , but also photosynthesis and respiration (Arens, Sullivan, and Welker 2008;Dagg and Lafleur 2011) , and soil organic carbon stocks (Campeau, Lafleur, and Humphreys 2014;Siewert 2018) .Furthermore, higher summer temperatures increase plant size and accelerate leaf economics in the tundra (Bjorkman, Myers-Smith, Elmendorf, Normand, Rüger, et al. 2018) , and instantaneous summer temperature is one of the most important drivers of photosynthesis and respiration due to its control on enzymatic processes (Natali et al. 2011;López-Blanco et al. 2017) .
Winter temperatures play an important role in the tundra ecosystems too.Increases in plant-experienced winter temperatures due to increased snowfall (Hagedorn et al. 2014) or winter warming (Harsch et al. 2009) have been linked to treeline advances, that is, increased plant size.On the other hand, shallow snow-packs, which correlate with lower plant-experienced winter temperatures, favor slower leaf economics (Choler 2005) .Low winter temperatures might also damage plants and therefore decrease growing-season photosynthesis (Bokhorst et al. 2011;Phoenix and Bjerke 2016) .However, the relative importance of these abiotic variables on plant functional composition and diversity as well as carbon cycling is not yet fully understood.
In this study, we examine the relationships between abiotic environmental conditions, plant community functional composition and diversity, and carbon cycling in a tundra landscape using a systematic study setting of up to 220 intensively measured plots and regression models.Moreover, we simulate how summer and winter warming could cause cascading effects on the carbon uptake capacity of the ecosystem, measured as the peak-season CO 2 budget.We explain the functional properties of communities and carbon cycling with environmental variables describing soil resources, and summer and winter temperatures which have been shown to explain functional composition in this system (Happonen et al. 2019) .To understand carbon cycling comprehensively, we focus not only on CO 2 fluxes (gross primary productivity, ecosystem respiration, soil respiration), but also on carbon stocks (above-ground and soil organic carbon stocks; for abbreviations used in Materials and Methods, and Results, see Table 1).Fluxes were measured during the peak growing season, which for plants is the most active season with the largest net CO 2 uptake (Belshe, Schuur, and Bolker 2013) .Stocks, on the other hand, reflect the balance of carbon uptake and losses over a longer period of time.

Study site
The field observations for this study were collected in 2016-2018 in an oroarctic tundra environment (Virtanen et al. 2016) in Kilpisjärvi, northwestern Finland (Fig. 1).The study area is located in an elevational gradient between two mountains, Saana (1029 m.a.s.l) and Korkea-Jehkats (960 m.a.s.l), and the valley in between (ca.600 m.a.s.l.).The study area is above the mountain birch ( Betula pubescens var.pumila (L.) Govaerts) forest line, and is predominantly dwarf-shrub heath, with graminoid-and herb-rich meadows concentrated around meltwater streams.Empetrum nigrum L., Betula nana L., Vaccinium myrtillus L., Vaccinium vitis-idaea L. and Phyllodoce caeruleae (L.) Bab. are highly abundant in the area with graminoids and herbs such as Deschampsia flexuosa (L.) Trin.and Viola biflora L.
The study design consisted of 220 locations in a 1.5 x 3.0 km area (Fig. 1).The distance between two adjacent locations was a minimum of 23 m (average 101 m).Individual study locations were thus in separate vegetation patches.Due to practical reasons, the full set of environmental, trait, and carbon variables were measured at only 80 locations.Although we could have used this data set from those 80 locations alone in this study, we wanted to make use of all the field data collected to characterize the environment as accurately as possible.
Thus, the number of sampling locations for each analysis varied between 80 and 200 (Supplementary Table 1, Fig. S1-S2).At the sampling locations where all variables were collected, soil moisture, air and soil temperature, and organic and mineral layer depth were measured at the central plot of the measurement scheme.Plant functional traits, CO 2 fluxes, and above-ground carbon stocks were measured in another plot 1-3 meters away from the central plot to avoid perturbing these permanent plots.Finally, soil samples were collected from the vicinity of the trait-carbon plot.Below we describe our sampling scheme, but to keep the Methods section concise and clear, some of the details (e.g., the number of sampling locations or the name of the measurement device) are only given in the Supplementary Material.
Abiotic variables: PAR, microclimate, soil resources At the study area, we measured Photosynthetically Active Radiation (PAR) as photosynthetic photon flux density (PPFD, µmol m 2 s -1 ) continuously every 10 minutes from July 8th to August 17th 2017 (Methods S1 and Supplementary Table 2).Soil moisture was measured during the growing seasons of 2016-2018 3-6 times per year as volumetric water content (%) ( (Happonen et al. 2019;Kemppinen et al. 2018) , Supplementary Tables 1-2 ) .In each soil moisture plot, three measurements were taken, and the average was used in further analyses.
Temperature loggers were installed in 2016 to monitor temperatures 10 cm above and below ground, at 2-4 h intervals (Supplementary Tables 1-2).These single time-point measurements were consequently aggregated to monthly averages.
We assumed that the composition and diversity of plant communities and their potential to store and cycle carbon reflect average rather than instantaneous environmental conditions.
While instantaneous carbon fluxes certainly do respond to changes in e.g.temperature, fluxes were normalized to common light and temperature levels before modelling to break the link between fluxes and environmental conditions during their measurement (see Light response model and flux normalization).To quantify long-term mean summer and winter temperatures and soil moisture conditions, we averaged July air temperatures, February soil temperatures, and growing season soil moisture levels over 2016-2018 for each location using random effects models with the R package lme4 (version 1.1-18-1, (Bates et al. 2015) , Methods S2).
We chose July and February because they are respectively the warmest and coldest months in the area.We used soil instead of air temperature to characterize winter temperature conditions because during that time, most plants are protected by a thick layer of snow.
Samples of the organic soil layer were collected, after which the samples were freeze-dried, and pH was analyzed in the laboratory of the University of Helsinki, following the ISO 10390 standard (Supplementary Tables 1-2).We used soil pH as a surrogate for soil nutrient status (See eg.Björk et al. 2007;Hobara et al. 2016;Gough et al. 2000) .The average soil moisture and pH of the organic layer were highly correlated (ρ > 0.8).Hence, we reduced them to their first principal component, which we hereafter refer to as the soil resource axis.The final predictors used in the analysis were thus named summer temperature (average July air temperatures), winter temperature, (average February soil temperatures) and soil resources.

Vegetation data: trait and community data
We quantified vascular plant community species composition in the majority of the plots where CO 2 flux measurements were conducted using the point-intercept method during the peak growing season in 2017 (Supplementary Table 1).These data have previously been used to study the relationship between the environment and functional composition (Happonen et al. 2019) .A circular frame (20 cm diameter) with 20 evenly spaced pinholes was placed over the soil collar (see section CO 2 flux data) after CO 2 flux measurements.The total number of times each species touched pins (3 mm in thickness) lowered into the vegetation through the frame was counted and used as a measure of species abundance.
We measured plot-specific plant height to indicate plant size, and LDMC as well as SLA to represent the leaf economics spectrum for all species (Happonen et al. 2019) .Due to time constraints, we did not measure other traits.However, these three traits have been shown to be good proxies of the axes of plant trait variation (Bruelheide et al. 2018;Díaz et al. 2016) .
LDMC and SLA were highly correlated (⍴ = 0.8) .Thus, results for LDMC and SLA were highly similar, and to reduce redundancy, we chose LDMC to represent the leaf economics spectrum in this study.Trait measurements were done separately for all locations, and always on plant individuals within the soil collars in order to account for intraspecific trait variation.
Height was quantified as the height of the highest leaf on two random ramets within the collar.LDMC was measured from two leaf samples taken from two different ramets (Methods S3).Location-specific trait values for each species were quantified as the average of the two individual trait values.
As a measure of functional composition, we used community weighted means (CWM) of these two traits, which have been shown to correlate with environmental conditions and ecosystem functioning (Garnier et al. 2004) .We chose coefficient of variation (CV, standard deviation/mean) as the measure of functional diversity.Many traits have a log-normal distribution, meaning that natural variation within a community increases with average trait values (Bjorkman, Myers-Smith, Elmendorf, Normand, Thomas, et al. 2018) .The CV represents trait variation in relation to that which would be expected for a given CWM value, and can be thought of as normalized functional variation in relation to functional composition.In the analysis, we refer to these variables as height, LDMC, CV height and CV LDMC .
CO 2 flux data We measured CO 2 exchange using a static, non-steady state non-flow-through system (Livingston and Hutchinson 1995) composed of a transparent acrylic chamber (20 cm diameter, 25 cm height; Supplementary Table 1) during the growing season in 2017.The chamber included a small ventilator, a carbon dioxide probe, an air humidity and temperature probe, and a measurement indicator (Vaisala, Vantaa, Finland, Supplementary Table 2).In the chamber, CO 2 concentration, air temperature and relative humidity were recorded at 5-s intervals for 90 s.PAR was logged manually at 10-s intervals during the same period using a quantum sensor with a hand-held meter (Supplementary Table 2).
Steel soil collars (21 cm in diameter and 6-7 cm in height) were inserted in the soil at least 24 h before the measurements to avoid potential CO 2 flush from soil due to the disturbance caused by the installation of the collars.The soils in the study area are relatively rocky and have long horizontal roots growing through them, thus our collars were embedded only ca. 2 cm into the soil.To guarantee an air-tight seal, we sealed the edges of the collar using inert quartz sand.The chamber was placed on top of the collar and ventilated after each measurement.We progressively decreased the light intensity of NEE measurements from ambient conditions to ca. 80%, 50% and 30% PPFD by shading the chamber with layers of white mosquito net (n = 7-10).ER was measured in dark conditions (0 PPFD), which were obtained by covering the chamber with a space blanket (n = 3).To measure SR, all above-ground vascular plant biomass was clipped.The clipping was done ≥24 h before the measurements to avoid disturbance.We conducted three SR measurements in dark conditions in each plot.The soil CO 2 emissions consist partially of moss and lichen respiration, as we were not able to remove entirely all miniature cryptogams on the soil due to their tight integration with the soil particles.Each plot was visited twice at midday during the peak season in 2017.During the first time, plots were visited to measure NEE and ER between the 26 th of June and 27 th of July, and during the second time, plots were visited to measure SR between the 7 th and 15 th of August.Fluxes were calculated with linear regression (Methods S4).

Carbon stock data
We measured the depth of the soil organic and mineral layers from three points on each central plot using a metal probe (Supplementary Table 1).In the analyses, we used the mean value of the three point measurements to represent the organic and mineral layer depths in each sampling location.We collected samples of roughly 1 dl from the soil organic and mineral layers with metal soil core cylinders (4 to 6 cm diameter, 5 to 7 cm height).The organic samples were collected from the top soil, and mineral samples directly below the organic layer.The soil samples were taken on the 1 st to 31 st of August 2016 and 2017.We used a freeze-drying method to remove moisture from the samples.Bulk density (kg m -3 ) was estimated by dividing the dry weight by the sample volume.Total carbon content (C%) was analyzed using a Vario Elementar -analyzer (Supplementary Tables 1-2) or was derived with a loss of ignition method (Method S5).Before C% analysis, mineral samples were sieved through a 2 mm plastic sieve.Organic samples were homogenized by hammering the material into smaller pieces.
SOC for organic and mineral layer was estimated using Equation 1: SOC = C% / 100% ⨉ bulk density (kg m -3 ) ⨉ layer depth (m) (Eqn 1) Finally, organic and mineral layer stocks were summed together to calculate the total SOC stock.Plots with no soil were excluded from the analysis.For more details, see Methods S5.
Above-ground vascular plant biomass was collected between the 1 st and 10 th of August, 2017 from the collars.Biomass samples were oven-dried at 70°C for 48 h and weighed after drying.Above-ground carbon stocks (AGC) were estimated by multiplying the total biomass by 0.475 (Schlesinger, 1991) .In the analysis, we refer to these carbon stock variables as SOC We also used this model together with PAR and temperature logger data to simulate plot-specific NEE with 10 minute intervals for the 30 day time period between the 8 th of July and 7 th of August 2017.Temperature measurements were first interpolated from 2-4 h to 10 minute resolution.We summed the predictions to create an estimate of peak-season CO 2 budget.
SR was normalized to 20℃ with a Bayesian linear mixed model.The model included a plot-specific random intercept.A linear instead of an exponential temperature response was used, because it reduced unexplained variance.All priors were left as the brms defaults.We refer to this variable hereafter as SR.

Hierarchical model of tundra carbon cycling
We first explored whether carbon-cycling models should include both the direct effect of environmental conditions and the indirect effects of environmental conditions via trait composition and diversity.This is because we speculated that average environmental variables correlate strongly with trait composition and diversity, and would thus not provide any additional information to carbon-cycle models.We compared the predictive accuracy of models explaining carbon cycling with direct and indirect environmental effects against those with indirect effects only using pareto-smoothed importance sampling leave-one-out cross-validation (Vehtari, Gelman, and Gabry 2017) , using the R-package loo version 2.1.0(Vehtari et al. 2019) .Predictive accuracy was measured as expected log pointwise predictive density (ELPD).The models in this comparison were fitted to a common dataset, where none of the variables had missing observations (n=80).
After verifying that direct environmental effects were not needed to explain ecosystem functions, we built a hierarchical model of tundra carbon cycling by combining five submodels in a similar framework that is used in structural equation models.The submodels describe 1) environmental effects on trait composition and diversity (n = 87), 2) trait effects on CO 2 fluxes (GPP, ER, SR, n = 129), 3) trait effects on above-ground carbon stocks (n = 134), 4) trait effects on soil organic carbon stocks (n = 129), and 5) the sensitivity of peak-season CO 2 budget to photosynthesis and respiration (n = 102) (Fig. 2).The sample sizes of these models differed, because variables were measured with different intensity (Supplementary Table 1).All response variables were modelled with multiple linear regression.For example, in the first submodel, soil resources, and summer and winter temperatures were used as predictors for trait composition and diversity.Plant height, CO 2 fluxes, and carbon stock variables were log e -transformed before analyses to satisfy the assumption of homoscedasticity.When a submodel had more than one response variable, we also modelled residual correlations among the response variables.For these models, priors were left as the brms default weakly informative priors.All models were fitted with four MCMC chains of 2000 iterations.The first 1000 samples of each chain were discarded as warmup, leaving 4000 MCMC samples of each parameter of each model.
We used the hierarchical model to simulate warming effects on ecosystem capacity to sequester carbon via changes in plant community functional composition and diversity (Methods S7).We predicted the peak-season CO 2 budget for median environmental conditions, a winter warming scenario (+1℃ February soil), a summer warming scenario (+1℃ July air temperature), and a combination of the two warming scenarios.For the scenarios with summer warming, we ran the simulations with and without direct temperature effects on summer respiration.Direct temperature effects were included in the models by multiplying predicted ER by the exponential of the temperature sensitivity parameter (β Temperature , Eqn S2) from the previously fitted light-response model.We did not take into account residual uncertainty in the predictions.We propagated uncertainty by evenly sampling 500 predictions from the posterior distribution of each submodel, and used these values as the explanatory variables in the next model.Finally, we also thinned the predictions from the last model to 500 samples.For each scenario, we calculated the probability of the modelled ecosystem to be a net CO 2 sink.
Data and scripts used in this study are deposited into Zenodo (Happonen et al. 2020) .

Light response model
The light response model used to calculate the normalized CO 2 flux variables seemed to converge; no MCMC chains showed erratic behaviour under visual inspection, and rank-normalized potential scale reduction factors for all parameters stayed below 1.01.Visual inspection of the posterior predictive distribution showed a very good fit to the data, also evidenced by a Bayesian R 2 (Gelman et al. 2019) of 0.95.

Abiotic versus plant functional trait controls on carbon cycling
Models with only trait-related covariates predicted carbon cycling as accurately -or even more accurately -than models with both environmental and trait covariates (Fig. 3).This means that after conditioning on functional composition and diversity, long-term average environmental conditions did not add new information about carbon cycling.By the principle of parsimony, we therefore focus our analysis on the simpler models.

Environmental controls on traits
Environmental conditions were better predictors of functional composition than of functional diversity.Vegetation height responded positively to both summer and winter temperatures (Figs.4a, 5).There was less certainty about the effect of soil resources on vegetation height, but the link seemed to be negative with a posterior probability of 96% (Fig. S3).LDMC was negatively associated with soil resources and winter temperatures (Fig. 4a).
The only environment-functional diversity connection whose 95% credible interval did not overlap zero was the positive link between winter temperatures and CV LDMC .The residual correlation between LDMC and CV Height was negative (ρ = -0.26,CI = -0.46--0.05).No other residual correlations were found among trait variables.Bayesian R²s for the models explaining vegetation height, LDMC, CV Height and CV LDMC were 0.27, 0.41, 0.07, and 0.14, respectively.
Trait controls on CO 2 fluxes and carbon stocks GPP, ER and SR had very similar responses to plant community functional composition and diversity.All fluxes responded positively to vegetation height and CV LDMC , and negatively to LDMC (Figs. 4b,5), although the link between LDMC and ER was negative with only 83% probability (Fig. S4).Residual correlations between log-transformed flux pairs GPP-ER, GPP-SR, and ER-SR were 0.6, 0.3, and 0.5, respectively.Bayesian R²s for the models explaining GPP, ER and SR were 0.43, 0.40 and 0.36, respectively.
SOC was positively related to CV LDMC (Fig. 4d, Fig. S6).LDMC and SOC had positive residual correlations with a 97.5% probability.Bayesian R²s for the models explaining AGC and SOC were 0.39 and 0.15, respectively.
CO 2 budget sensitivity to warming The simulated peak-season CO 2 budget was 1.4 times more sensitive to changes in ER than in GPP (Fig. 4e).Fluctuating light and temperature levels and the nonlinear light response curve did not cause major nonlinearities in the determination of the budget on the basis of temperature-and PAR-standardized GPP and ER, as indicated by the 0.98 R² of the linear model.The expected peak-season CO 2 budget was positively affected by the functional changes brought about by simulated summer and winter warming, even when the direct effects of respiration were taken into account (Fig. 6).The probability of an ecosystem patch being a summer net CO 2 sink grew in all but one scenario.In the simulation with summer warming and a direct temperature effect, the probability of being a sink fell by one percentage point.

Discussion
Arctic vegetation is being reshuffled at unprecedented speed (Post et al. 2009) , both in predictable and surprising ways (G.J. Jia, Epstein, and Walker 2003;Phoenix and Bjerke 2016) .Thus, it is imperative that we form a solid understanding of the causal forces shaping Arctic ecosystems and their functioning.Here, by using a hierarchical modeling framework, we provide new insights into the relative roles of environmental conditions and plant functional traits in determining CO 2 fluxes and carbon stocks in the tundra.A functional trait-based approach proved very informative, since we were able to partition the effects of vegetation composition and diversity on carbon cycling to size and resource use strategy-related components.Functional composition and diversity on the recently identified axes of plant above-ground functional trait variation (Díaz et al. 2016;Bruelheide et al. 2018) thus have clear and separate effects on tundra carbon cycling.Since these community trait attributes also respond differently to environmental conditions, they are excellently suited to inform us about the impacts of global change and vegetation shifts on carbon cycling and resulting climate feedbacks.
Carbon cycling was most parsimoniously explained by trait composition and diversity alone, suggesting that average environmental conditions did not provide any new information for predicting carbon cycling (i.e., traits already included this information).However, instantaneous environmental conditions also had a strong impact on photosynthesis and respiration, as evidenced by the temperature sensitivity of the light-response curve.This

Plant community functional composition regulates fine-scale carbon cycling
Plant size had a strong and positive relationship with all fluxes and pools except for soil organic carbon stocks.The relationship between plant size and photosynthesis was stronger than that of plant size and ecosystem respiration.Evidence suggests that plant size is increasing around the Arctic in response to climate change.Sometimes this is in part due to "shrubification", the increase in the biomass, cover, and abundance of woody species (Myers-Smith et al. 2011) , but other large species are increasing as well (Bjorkman, Myers-Smith, Elmendorf, Normand, Rüger, et al. 2018) .Our study agrees with the results of Lafleur and Humphreys (2018) , who discovered that taller shrubs had larger growing season net uptake than lower shrubs in the low Arctic.They argued that this was likely due to colder soils limiting heterotrophic respiration and taller plants having a larger leaf area to capture sunlight.However, our result disagrees with the findings by Gagnon et al. , (2019) who showed that vegetation and soil organic carbon stocks were highest in tall shrub vegetation, as we found no link between plant size and soil organic carbon stocks.Nevertheless, from our results it follows that the documented changes in plant height around the tundra (Bjorkman, Myers-Smith, Elmendorf, Normand, Rüger, et al. 2018) have already caused large shifts in the functioning of Arctic ecosystems, especially on carbon fluxes and above-ground carbon storage, but the effects on total ecosystem carbon storage require further investigation.
The effects of the leaf economics spectrum on carbon cycling were more variable.Our study supports the well-established strong positive relationship between faster economic traits and photosynthesis (Williams et al. 2006;Street et al. 2007;Shaver et al. 2007) , but suggests that the connection to ecosystem respiration is not as strong.The relationship between ecosystem respiration and economic traits is more complicated due to the various sources of ecosystem respiration (above-ground plant tissues, roots, symbionts and heterotrophs) which all have different drivers (Segal and Sullivan 2014;Barba et al. 2018) .The positive effect of faster leaf economics on soil respiration was stronger and more certain than on ecosystem respiration.We speculate that there are three potential reasons for this.First, faster growing plants often produce more litter that is rich in nutrients and more easily broken down by soil microbes, resulting in higher soil respiration rates (Cornwell et al. 2008;Freschet, Aerts, and Cornelissen 2012) .Second, faster plants produce more root exudates, which are, together with root respiration, suggested to be a main source of soil respiration in the tundra (Illeris, Michelsen, and Jonasson 2003) .Third, plants with high leaf nutrient concentrations also have high nutrient concentrations in the roots, as suggested by the whole-plant economic spectrum hypothesis (Freschet et al. 2010) , which might also accelerate soil respiration by facilitating root respiration (S.Jia et al. 2013) .Thus, there are complex background processes that could be associated with the relationships that we have observed.
Although above-ground traits were linked to soil respiration, they were poorly linked to soil organic carbon stocks, and a large part of variation in soil organic carbon stocks remained unexplained.It is possible that these weaker links are in part due to soil organic carbon being measured outside the vegetation collar, albeit in visually similar vegetation and soil conditions (see Fig. S2).Soil organic carbon stocks could also be weakly related to traits because they are a result of carbon inputs from above-and below-ground plant tissues, outputs from ecosystem and soil respiration, and lateral transport, which all respond differently to trait variation.Trait effects on soil organic carbon stocks are also more indirect (e.g. via litter decomposition and accumulation) compared to, for example, the relationship between leaf economics and photosynthesis (De Deyn, Cornelissen, and Bardgett 2008) .In addition to above-ground vascular plant traits, soil organic carbon stocks could be explained by the traits of mosses, which can play an important role in wetter landscape positions (Douma et al. 2007) , or plant litter and root traits which are more directly linked with soil organic carbon stocks (Orwin et al. 2010) .Soil organic carbon stocks are also likely influenced by additional abiotic drivers, such as aeolian and fluvial transportation of litter (Klaminder, Yoo, and Giesler 2009) , or soil development and carbon accumulation in the past (Palmtag et al. 2015) .Thus, comprehensive understanding on the causes of variation in soil organic carbon stocks can probably only be achieved by measuring static environmental and biological conditions along with the effects of matter flows across the landscape.
Earlier studies in the tundra have often characterized the effects of vegetation functional composition on carbon cycling with LAI (Street et al. 2007;Shaver et al. 2007;Dagg and Lafleur 2011;Marushchak et al. 2013) .Measuring LAI does not require information on species assemblages, making it useful for landscape-scale or temporal predictions of carbon cycling (Shaver et al. 2007;Cahoon, Sullivan, and Post 2016) .However, LAI is a composite of several trait axes as, for example, increases in plant size and SLA can both increase LAI (Fletcher et al. 2012) .Consequently, understanding the mechanisms driving the relationships between the environment and LAI, or LAI and carbon cycling, can be challenging.Our study adds to the results of these earlier studies by partitioning plant community functional composition to community position on the two main axes of trait variation, thus improving our mechanistic understanding about vegetation-carbon linkages in the tundra.

Functional diversity increases CO 2 fluxes and carbon stocks
All fluxes and the soil organic carbon stock responded positively to the functional diversity of leaf economics, but not to the diversity of plant size.Thus, within-community variability on the slow-fast strategy axis increased ecosystem functioning, supporting the niche complementarity hypothesis.There are several mechanisms for how soil organic carbon stocks could increase with plant functional diversity.One such mechanism is that increased productivity could lead to higher litter inputs from leaves and roots (DeMarco et al. 2014) .
Since leaf economic diversity did not increase above-ground carbon stocks, our interpretation is that the positive diversity effect on soil organic carbon stocks is probably due to increased root litter production, indicative of higher total investments in acquiring below-ground resources.These observations would be consistent with more efficient partitioning of below-ground resources, which would also explain the positive effects of leaf economic diversity on CO 2 fluxes.
Functional diversity of plant size was positively associated with above-ground carbon stocks, indicating that multi-layered vegetation stores more carbon compared to less structurally diverse plant communities.Similar effects of size-related trait diversity on above-ground biomass have been reported for temperate semi-natural grasslands (Schumacher and Roscher 2009) .On the other hand, Conti and Díaz (2013) reported a negative correlation between structural diversity and above-ground carbon stocks in subtropical forests, but note that this is probably due to confounding factors.Our results provide strong evidence that structurally layered tundra communities do indeed contain more carbon than vertically homogeneous vegetation patches, because we modelled the effects of structural and resource-use trait diversities separately, and because we controlled for the direct effects of functional composition in constructing our diversity measures.
Effects of warmer temperatures on peak-season CO 2 budget Simulated warming without including a direct temperature effect on respiration increased the capacity of ecosystems to sequester carbon during peak growing season.This held true for both summer air and winter soil temperatures.These effects manifested through different mechanisms: summer warming increased photosynthesis by increasing plant size, while winter warming both increased plant size and accelerated leaf economics.This demonstrates how environmental conditions outside the growing season are important for understanding changes in high-latitude ecosystem functioning (Happonen et al. 2019) .
Our results suggest that +1℃ higher summer air temperatures increase summer ecosystem respiration by ca.3%.Including this direct temperature effect on peak-season CO 2 budget did not reverse the expected summer warming effects on the budget, but it did slightly increase the probability of vegetation turning into a CO 2 source.Results from other parts of the tundra have shown that short-term warming can in fact decrease net carbon uptake during the growing season by increasing respiration more than photosynthesis (Biasi et al. 2008;Parmentier et al. 2011) .Based on our results, this is also true in Scandinavian mountain tundra, but in the absence of other drivers, the long-term effects of warming on peak-season CO 2 budget are more likely to increase net sink strength due to warming effects on plant functional traits.These positive indirect effects are not certain to be realized, however.For example, reindeer grazing has been shown to buffer tundra communities against the expansion of taller species (Vuorinen et al. 2017;Kaarlejärvi, Eskelinen, and Olofsson 2013) .
Even without biotic buffering, dispersal limitation and other inertia in the population dynamics of long-lived tundra species can cause time lags in vegetation transitions (Alexander et al. 2018) , whereas accelerated respiration in response to temperature increase is more or less immediate.Taken together with the fact that respiratory carbon losses are bound to increase also beyond the peak growing season (Natali et al. 2019) , our results cannot be interpreted to mean that future warming will increase the annual CO 2 sink of the tundra.

Conclusions
We found that in a tundra landscape, fine-scale carbon cycling was explained by plant community functional composition and diversity, which mediated the effects of average environmental conditions on carbon cycling.Average plant size was the strongest predictor of most carbon cycling variables, but leaf economic diversity mattered the most for soil organic carbon stocks.Plant size increased CO 2 fluxes and above-ground carbon stocks, whereas fast leaf economics were associated with higher CO 2 fluxes and lower above-ground carbon stocks.All carbon cycling variables were positively affected by the diversity of some functional trait.Both warmer summer air and winter soil temperatures would lead to functional changes in plant communities that would enhance carbon uptake during the peak season.By utilizing globally applicable plant functional traits, our findings facilitate forming an integrated view on the causes and ecosystem consequences of climate change-related vegetation transitions in the tundra.

Photosynthesis
The process by which plants build carbohydrates from light and CO 2 .GPP Gross primary productivity normalized at a common irradiance (600 PPFD, partly sunny conditions)

Ecosystem respiration
The process by which plants, soil animals and microbes convert sugars into energy, thus releasing CO 2 to the atmosphere.

ER
Ecosystem respiration normalized at a common air temperature (20 ° C, average air temperature)

Soil respiration
The sum of respiration of all heterotrophic organisms (decomposition by animals, fungi, and microbes) and root respiration.

SR
Soil respiration normalized at a common air temperature (20 ° C, average air temperature) Peak-season CO 2 budget CO 2 sequestration by gross primary productivity minus CO 2 release from ecosystem respiration accumulated over a certain period of time.
Peak-season CO 2 budget GPP and ER were modeled over a 30-day period in July and August.Then, cumulative fluxes were calculated and NEE was estimated as GPP-ER.

Above-ground carbon stock
The carbon content in the living plant biomass.AGC Above-ground carbon stocks were estimated by multiplying the above-ground vascular plant biomass by 0.475.

Soil organic carbon stock
Soil organic carbon stocks in the entire soil horizon excluding roots.

SOC
Soil organic carbon stocks were estimated using the information of soil organic and mineral layer depth, bulk density, and carbon or organic matter content.
and AGC.Light response model and flux normalizationAll models were run in R version 3.5.3with the package brms version 2.8.0 (Bürkner, 2018) , which is an interface to the Bayesian modelling platformStan (Carpenter et al. , 2017)  .We fitted plot-specific light-response curves of NEE using the Michaelis-Menten equation in a non-linear hierarchical Bayesian model.The intercept term (dark respiration) was fitted with an exponential temperature response.For a full description of the model, see Methods S6.We used this model to predict the rate of NEE at dark (0 PPFD) and average light (600 PPFD) conditions, and a temperature of 20℃ at each plot.Normalizing removes the effects of variation in instantaneous light and temperature conditions within the landscape from the resulting variables.We refer to the flux normalized to dark conditions as ecosystem respiration (ER).We then subtracted ER from the NEE normalized to average light conditions to arrive at an estimate of normalized GPP.GPP, ER, and SR are always shown with positive signs.Positive numbers for NEE and CO 2 budget indicate net CO 2 gain to the ecosystem (i.e.CO 2 sink) and negative numbers indicate net CO 2 loss to the atmosphere (i.e.CO 2 source).
apparent mismatch is expected, since our normalized carbon fluxes were cleaned of the effects of varying temperature and light levels.Separating realized carbon fluxes to potential (normalized fluxes governed by average environmental conditions and plant community characteristics) and realized components (governed by flux potential and instantaneous environmental conditions) was necessary to allow us to study the relative effects of instantaneous and average environmental conditions as well as plant functional community composition and diversity on carbon fluxes.
Figure 1.The study site.The vegetation is heterogenous.On the rocky slopes, the vegetation is sporadic (a).Intermediately productive sites are dominated by both deciduous (b: Betula nana, Vaccinium myrtillus) and evergreen (d: Empetrum nigrum) ericoid-graminoid heaths.Lush meadows can be found in the valley bottom (c).Orange dots represent the 220 study locations.The orthophoto (0.5 m resolution) is provided by the National Land Survey of Finland under a Creative Commons Attribution 4.0 licence ( https://www.maanmittauslaitos.fi/en/maps-and-spatial-data/expert-users/product-descriptions/orthophotos , accessed on 2019-10-28.).Photographs by Julia Kemppinen.

Figure 2 .
Figure 2. Our theoretical framework, where traits act as mediating links between environmental conditions and carbon cycling.Large arrows represent the submodels, which describe environmental effects on trait composition and diversity, trait effects on CO 2 fluxesand carbon stocks, and the sensitivity of peak-season CO 2 budget to photosynthesis and respiration.When a submodel had more than one response variable, we also modelled residual correlations among the response variables, shown in thin arrows.Submodels used all explanatory variables (e.g.all carbon cycling submodels included all traits as covariates), except for the peak-season CO 2 budget model (marked with an asterisk), which omitted soil respiration, since it is a subcomponent of ecosystem respiration.

Figure 3 .
Figure 3. Leave-one-out predictive accuracy of CO 2 flux (CF), above-ground carbon stock (AGC) and soil organic carbon stock (SOC) models with traits and environmental variables or with traits alone as predictors.Predictive accuracy was measured as expected log pointwise predictive density (ELPD).The lines are 95% credible intervals around the point estimate.

Figure 4 .
Figure 4. Modelled responses of a) plant community functional composition and diversity to the environment, b) CO 2 fluxes, c) above-ground carbon stocks and d) soil organic carbon stocks to plant community functional composition and diversity, and e) peak-season CO 2 budget to GPP and ER.Each panel uses 500 draws from the posterior distribution of the relevant parameter to simulate the marginal effect of one predictor when all other predictors are at their means.Height = community-weighted mean plant height, LDMC = community-weighted mean leaf dry matter content, CV height = coefficient of variation of plant height, CV LDMC = coefficient of variation of leaf dry matter content, Summer temperature = average July air temperatures, Winter temperature = average February soil temperatures, Soil resources = first principal component of soil moisture and pH), GPP = gross primary production normalized at a common irradiance (600 PPFD), ER = ecosystem respiration normalized at a common air temperature (20℃), SR = soil respiration normalized at a common air temperature (20℃), peak-season CO 2 budget = cumulative CO 2 balance over July.

Figure 5 .
Figure 5.A dichotomized representation of the causal network implied by our hierarchical carbon cycling model.Lines show statistically significant relationships between the variables (α = 0.05, red = positive effect, black = negative effect).Variables on the left explain those on the right.Double-headed arrows represent residual correlations.Positive numbers for CO 2budget indicate net CO 2 gain to the ecosystem (i.e.CO 2 sink).Low LDMC (leaf dry matter content) implies fast leaf economics, hence for ease of interpretation, all connections involving LDMC have been inverted for this figure (e.g. from negative to positive).

Figure 6 .
Figure 6.The effects of simulated warming on peak-season CO 2 budget with and without a direct temperature effect.Positive numbers for CO 2 budget indicate net CO 2 gain to the ecosystem (i.e., CO 2 sink).