Intraspecific trait variability is a key feature underlying high Arctic plant community resistance to climate warming

In the high Arctic, plant community species composition generally responds slowly to climate warming, whereas less is known about the community functional trait responses and consequences for ecosystem functioning. The slow species turnover and large distribution ranges of many Arctic plant species suggest a significant role of intraspecific trait variability in functional responses to climate change. Here we compare taxonomic and functional community compositional responses to a


INTRODUCTION
In the rapidly warming Arctic (IPCC, 2021), warminginduced vegetation changes may mediate tundra ecosystem functioning (Mekonnen et al., 2021), potentially with positive climate feedback consequences.In general, climate warming has stimulated Arctic plant growth in recent decades, reflected by plant community changes, increased plant abundance, and taller plants across the Arctic in local plots (e.g., Bjorkman et al., 2018;Elmendorf, Henry, Hollister, Björk, Boulanger-Lapointe, et al., 2012a), as well as at larger scales (Bhatt et al., 2017;Goetz et al., 2005).However, Arctic plant communities do not respond evenly to climate warming across the tundra biome or across habitats within sites (Bhatt et al., 2017;Bjorkman et al., 2020;Elmendorf, Henry, Hollister, Björk, Boulanger-Lapointe, et al., 2012a;Myers-Smith et al., 2020).What drives heterogeneity in warming responses in the tundra and the consequences of this variability for ecosystem functioning and resilience to climate change are not fully understood.
Overall, warming-induced changes in plant community structure and species composition have been more prevalent at relatively warm low Arctic sites than at colder high Arctic sites (e.g., Elmendorf, Henry, Hollister, Björk, Bjorkman, et al., 2012b;Elmendorf, Henry, Hollister, Björk, Boulanger-Lapointe, et al., 2012a).This may seem surprising given the higher rates of warming in the high Arctic (IPCC, 2021), but it can be attributed to several interacting biotic and abiotic factors.The small species pools in the high Arctic relative to the low Arctic consist of species under more severe life-history constraints (Arft et al., 1999) that are adapted to more significant interannual climatic variability and shorter growing seasons with longer photoperiods (Hudson & Henry, 2010;J onsd ottir, 2005).At smaller spatial scales, greater plant community responsiveness has been observed in wet relative to dry tundra habitats within both low Arctic (Bjorkman et al., 2020;Elmendorf, Henry, Hollister, Björk, Boulanger-Lapointe, et al., 2012a) and high Arctic sites (Edwards & Henry, 2016;Hudson et al., 2011).It is unknown, however, whether and how the observed resistance of species composition to climate warming in the high Arctic is associated with responses in plant community functional composition.
Trade-offs limit the extent of variation in functional traits that affect plant fitness, resulting in relatively few essential trait combinations (Díaz et al., 2016).Globally, most of the aboveground trait combinations (i.e., the functional space) can be captured by two trait dimensions, size (of whole plants as well as their parts) and leaf economics spectrum, which balance leaf construction costs against growth potential (Díaz et al., 2016).In turn, plant functional traits may affect ecosystem carbon sequestration, litter decomposition, and nutrient cycling (e.g., Díaz et al., 2004;Sørensen et al., 2019).Although tundra plant species only occupy a small part of the global trait space along the size dimension, they demonstrate almost the full trait variation along the leaf economics spectrum.However, this variability also becomes increasingly constrained toward the colder parts of the tundra (Thomas et al., 2020).
Recent climate warming of the tundra biome has been associated with a general shift in Arctic plant communities toward taller plants.In contrast, the strength and direction of warming responses of other size traits and leaf economics are more variable, governed by local conditions such as soil moisture (Bjorkman et al., 2018).Less is known about whether and how other aspects of the Arctic environment influence plant functional trait responses to warming and to what extent intraspecific trait variation (ITV) contributes to those responses.For instance, snow regimes in Arctic landscapes are strongly shaped by topography, affecting both growing season length and moisture (Rixen et al., 2022;Walker, 2000) and may, in turn, affect tundra plant community diversity and functional responses to warming (e.g., Niittynen et al., 2020).
To date, most studies of plant functional traits in tundra systems have relied on species-level trait information and so ignore ITV (e.g., Bjorkman et al., 2018;Kemppinen et al., 2021).ITV can reflect both phenotypic plasticity within an individual and genetic differentiation within a population (Albert et al., 2010;Bolnick et al., 2011), offering two potential mechanisms for functional responses to environmental change in systems with limited taxonomic change (i.e., plastic responses and differential mortality, respectively).Globally, ITV accounts for around 25% of plant community trait variation (Siefert et al., 2015), with an even more significant relative contribution locally (Messier et al., 2017;Thomas et al., 2020).In general, ITV should be relatively more important in communities with low species richness (Siefert et al., 2015;Thomas et al., 2020) and in communities dominated by species with broad environmental niches (Bolnick et al., 2011;Sides et al., 2014) or high phenotypic plasticity (Albert et al., 2010).According to the climate variability hypothesis (Spicer & Gaston, 1999), greater temperature variability selects organisms with broad thermal tolerances.There is evidence that geographic range size is related to niche breadth (Slatyer et al., 2013).In the high Arctic, species richness is relatively low, and most plant species have wide distribution ranges and are adapted to high seasonal climate variability (e.g., J onsd ottir, 2005), characteristics that predict relatively large ITV.Furthermore, most Arctic plant species are long-lived (J onsd ottir, 2011), which may slow down species turnover in response to climate warming, thereby enhancing the resistance of plant community species composition to warming and further increasing the relative importance of ITV in functional responses.
One of the most critical climate-related services that permafrost-affected ecosystems provide relates to carbon storage (IPCC, 2021).Abiotic drivers are usually considered the primary determinants of the ecosystem functions associated with this climate feedback property (Schuur et al., 2015).However, plant functional trait composition also affects carbon balance.At the individual level, leaf chemical and structural traits strongly influence carbon-cycling processes as plants toward the fast-return end of the leaf economics spectrum produce thin and nutrient-rich leaves (e.g., leaves with high specific leaf area [SLA] and N content; Wright et al., 2004) that promote higher photosynthesis rates (Reich et al., 1997) and fast decomposition by soil microbes (Funk et al., 2017).Leaf economics also scales to the community level, where canopy N and leaf area index correspond well with primary production rates (Reich et al., 2014).In Arctic tundra, habitats dominated by different plant functional groups exhibit differences in carbon cycling processes.Specifically, tall deciduous shrubs accelerate CO 2 exchange and decomposition rates either directly or indirectly through internal ecosystem feedback (e.g., Christiansen et al., 2018;Kropp et al., 2020;Lafleur & Humphreys, 2018;Myers-Smith et al., 2011), suggesting that functional traits should relate well to ecosystem function.Indeed, moving beyond functional group classifications, recent studies showed that community means of leaf economics and plant height affect tundra ecosystems in their ability to fix CO 2 (Happonen et al., 2022;Sørensen et al., 2018) owing to (1) the trade-offs associated with being a "fast" or "slow" species (Reich, 2014;Wright et al., 2004) and (2) plant height scaling positively with size and, therefore, photosynthetic tissue, i.e., leaf area (Lafleur & Humphreys, 2018).Consequently, climate warming-induced plant community changes, leading to tundra plant communities with taller plants (Bjorkman et al., 2018), will likely cause significant shifts in ecosystem properties related to carbon cycling.These changes should be quantifiable and predictable using functional trait approaches-alone or in concert with environmental drivers (Díaz et al., 2007).Nevertheless, the utility of functional approaches in predicting ecosystem responses has recently been questioned (Funk et al., 2017), and few studies have investigated the direct effects of functional trait change and climate warming on tundra CO 2 exchange rates.
This study investigated the functional responses of high Arctic vegetation to long-term experimental warming.Specifically, we asked whether and how the apparent resistance of high Arctic plant communities to warming was mediated by ITV and to what extent such trait variability within species has consequences for ecosystem functioning.To address these questions, we analyzed plant community composition data from a long-term (17-year) warming experiment replicated across three major habitats shaped by topography and contrasting snow regimes in high Arctic Svalbard.In the final year, we augmented these community data with detailed plot-scale measurements of 12 size-and leaf economics-related plant functional traits for all vascular plant species present.We also measured peak-season carbon fluxes in the experimental plots in the final study year to explore the potential roles of the environment, plant community taxonomy, and plant functional trait composition in mediating ecosystem responses to warming.
We hypothesize that (1) in line with other high Arctic sites, the taxonomic composition of the plant community does not change, or only weakly changes, in response to long-term experimental warming.Based on the foregoing literature review, we further hypothesize that (2) the plant functional trait composition is affected by experimental warming and (3) if both 1 and 2 are true, then we expect most community-level functional trait responses to experimental warming to be driven by ITV.Finally, we expect (4) rates of photosynthesis and respiration to be governed by the realized trait variation across habitats and climate warming treatments.Incorporating ITV should thus improve model performance and explanatory power for models of ecosystem carbon fluxes.

Study area
Annual temperature (1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) at Svalbard Airport (10 km northwest of the experimental site) averaged À4.6 C, and annual precipitation averaged 191 mm (Førland et al., 2011).The linear trend in mean annual temperature during 1889-2018 was an increase of 0.32 C per decade and by 1.66 C per decade from 1991 to 2018 (Nordli et al., 2020).Similarly, summer temperatures (June-August) increased by 0.13 C per decade during 1889-2018 and by 0.66 per decade during 1991-2018.The experimental site (78 11 0 N, 15 45 0 E) is situated at the south-southeast-facing slope of the Endalen valley at approximately 80 m above sea level, where distinct habitats are shaped by topography and snow accumulation in winter.The experiment includes three common Arctic habitats: a snowbed where deep snow accumulates (>100 cm) causing late snow melt (on average in mid-June), a Cassiope heath at intermediate exposure, snow depth, and snowmelt time, and a more exposed Dryas heath with only shallow snow (up to 10 cm) and early snowmelt (on average in mid-May).Consequently, the growing season duration varies between habitats and is ~2.5, 3, and 3.5 months in the snowbed, Cassiope, and Dryas heath, respectively.The soils are typical Cryosols with a thin organic layer on top of inorganic sediments (Jones et al., 2010).Soil moisture is related to topography and increases with snow accumulation in winter.
In the snowbed, total vascular plant cover is low, characterized by the forb Bistorta vivipara, the deciduous prostrate dwarf shrub Salix polaris, and the grasses Poa arctica and Festuca richardsonii.The moss cover is around 70%, relatively deep (>5 cm), and dominated by Tomentypnum nitens and Sanionia uncinata.In the Cassiope heath, the erect evergreen dwarf shrub Cassiope tetragona dominates.Other abundant vascular plant species are S. polaris and B. vivipara, and the most abundant moss species are T. nitens and S. uncinata.In the Dryas heath, the evergreen dwarf shrub Dryas octopetala dominates.Other abundant vascular species are Carex rupestris, B. vivipara, and S. polaris.The moss layer is well developed but shallower than in the other habitats, with S. uncinata as the most abundant species.

Experimental design
In 2001, ten 75 Â 75-cm (0.56 m 2 ) permanent plots were established in each of the three habitats.Five plots in each habitat were randomly assigned to a warming treatment in the following year.The other five served as controls.Climate warming was simulated by the use of hexagonal Plexiglas open-top chambers (OTCs) (Molau & Mølgaard, 1996), with a basal diameter of 150 cm (ca.1.8 m 2 ) and 47 cm in height, leaving ample space inside the chamber for plant sampling for trait measurements without disturbance to the permanent plot or risking edge effects.Following rain-on-snow events in the winter of 2008-2009, a large proportion of the Cassiope shrubs were killed by basal ice formation in two control plots and two OTC plots in the Cassiope heath.
Such extreme events are interesting aspects of climate change.Still, since the paper primarily addresses the effects of experimental warming on the plant communities, these plots were excluded from all subsequent analyses.

Environmental variables
Climate data were obtained for the study area during the last 3 years of the study by an automatic weather station (HOBO H21-002, Bourne, MA, USA) placed in the Dryas heath habitat, measuring air temperature (HOBO S-THB-M008) and photosynthetic radiation (PAR) (HOBO PAR S-LIA-M003) at 2 m height above ground and soil water content (HOBO S-SMC-M005) at 5 cm depth.Surface temperature and soil temperature were monitored over the year in each plot during two different time periods using TinyTag Plus data loggers (Gemini Data Loggers, Chichester, UK) in 2004-2005 (soil temperature at 10 cm) and iButton data loggers (DS1922L-F5 thermochrons, Maxim Integrated, San Jose, CA, USA) in 2015-2018 (soil temperature at 5 cm).Volumetric soil moisture was obtained from all 30 plots in 2018 with a handheld ML3 ThetaProbe (Delta-T Devices, Cambridge, UK).

Plant community assessment in taxonomic space
In all plots, a detailed vegetation analysis using the point intercept method was performed in 2003, 2009, and 2015, following standard protocols of the International Tundra Experiment (Molau & Mølgaard, 1996).We used 100 points per plot and recorded all hits (intercepts) through the canopy in each point, down to the bryophyte and lichen layer.If no bryophytes or lichens were present, the last hit was litter, biocrust, rock, or soil surface.Canopy height, including reproductive structures, was recorded in each point in 2015.Vascular plants were all recorded to the species level and most of the bryophytes and lichens (some only identified to genus level).

Plant community assessment in functional space
We focused on 12 plant functional traits to investigate plant community responses to warming along the two main dimensions of the global plant functional trait space (e.g., Díaz et al., 2016).Three traits reflected plant size (plant height, leaf dry mass, leaf area), and nine reflected the leaf economics spectrum (SLA, leaf thickness, leaf dry matter content [LDMC], nitrogen content [%N], carbon content [%C], phosphorus content [%P], carbon: nitrogen ratio [C:N], nitrogen isotope [δ 15 N], and carbon isotope ratio [δ 13 C]).The isotope ratios were included because they may indicate the water status of the plants (C) and nutrient availability along environmental gradients (N) (Pérez-Harguindeguy et al., 2013), and both ratios were therefore expected to respond to the warming treatment.For each plot, up to three individuals of all plant species covering more than 1% of the plot were selected in mid-July 2018.To avoid destructive sampling within the permanent long-term monitoring plots, the plant material for trait measurements was collected outside these plots but within the much larger OTCs, and in the surroundings close to each control plot.If more than three individuals of a species were available, the individuals to be sampled were selected using a randomization procedure.The height of each sampled plant individual was measured from the ground to the top of its highest photosynthetic leaf (standing, unstretched height), excluding inflorescences.The plant's aboveground parts were then sampled, placed in a ziplock plastic bag, labeled by plot ID, date, height, taxon, and sample ID, and transported to the lab.There was no precipitation during sampling that could complicate subsequent fresh weight measurements.The samples were stored outdoors and in the shade (temperature around 6 C) to ensure they stayed water saturated until the measurements started, within 24 h of collection.The procedure for trait measurements followed standard trait measurement protocols (Pérez-Harguindeguy et al., 2013).
Leaves (including petioles and stipules for Dryas octopetala, if present, but excluding sheaths for graminoids) were separated from the plants using tweezers and scalpels and sorted into paper envelopes with a standardized labeling system.Since leaf sizes were generally small and shape varied between plant species within the study site, standardized rules were applied.For each plant, up to three healthy and mature leaves were sampled and stored in a separate envelope with a unique sample ID.If the leaves of a species were tiny (e.g., Dryas octopetala, Saxifraga oppositifolia, or Salix polaris), several leaves equivalent to an area of ~3 cm 2 were collected for each sample.Fresh, moist leaves were weighed (in grams, rounded to four decimals; Mettler AE200, Mettler TOLEDO, and AG204 DeltaRange [0.1 mg precision]) for wet mass.To estimate leaf area, leaves were scanned (Canon Lide 220, resolution 300 dpi), and the leaf area was calculated using ImageJ (Schneider et al., 2012) and the LeafArea package (Katabuchi, 2015) in R version 4.0.2.Three measurements of leaf thickness per sample were taken, when possible, using a digital caliper Neuss,Germany).Measurements on the leaf veins were avoided, although this was not possible for tiny leaves.The leaf samples were then dried at 60 C to a constant mass (24-28 h) in a drying cupboard (Thermo Scientific Heraeus, USA) and weighed for dry mass.We calculated SLA as the leaf area divided by the dry mass (cm 2 /g) and LDMC as the dry mass divided by the fresh mass (g/g).
Following the size and mass measurements, the leaves were ground to a powder and analyzed for nutrients (P, N, C) and isotope ratios (δ 15 N, δ 13 C) at The University of Arizona.Total phosphorus concentration was determined using persulfate oxidation followed by the acid molybdate method (APHA, 1992).Phosphorus concentration was then measured colorimetrically with a spectrophotometer (Thermo Scientific Genesys20, Waltham, MA, USA).Carbon (C), nitrogen (N), and their stable isotope ratios were measured by the Department of Geosciences Environmental Isotope Laboratory at the University of Arizona on a continuous-flow gas-ratio mass spectrometer (Finnigan Delta PlusXL, Waltham, MA, USA) along with an elemental analyzer (Costech, Valencia, CA, USA).Samples of 1.0 AE 0.2 mg were combusted, and standardization was based on acetanilide for N and C concentration, NBS-22 and USGS-24 for δ 13 C, and IAEA-N-1 and IAEA-N-2 for δ 15 N.The carbon-to-nitrogen ratio (C:N) was also calculated and analyzed.

Ecosystem CO 2 flux measurements
In mid-July 2018, we measured peak growing season net ecosystem CO 2 exchange (NEE) and ecosystem respiration (R eco ) using an infrared gas analyzer (Li-840, LI-COR Biosciences, Lincoln, NE, USA) connected to a custom-made Plexiglass chamber (headspace volume = 25 L).We used a polyethylene skirt and a heavy chain to create an airtight seal between the chamber and the atmosphere.Two fans fitted within the chamber ensured sufficient air circulation during measurements.Daytime fluxes were measured twice per experimental plot, on separate days.For each plot and day, we measured NEE during ambient light conditions, followed by thorough aeration of the chamber headspace for 10-15 s before measuring R eco by covering the chamber with an opaque polyethylene cloth.Each flux measurement lasted for at least 120 s, logging CO 2 concentration values every second.During each flux measurement, we collected data on microclimatic parameters, i.e., soil temperature at 5 cm depth, canopy temperature, volumetric soil moisture integrated across 0-10 cm depth, and PAR using handheld thermometers (Fisher Scientific, Oslo, Norway), an infrared thermometer (Biltema, Åsane, Norway), an ML3 ThetaProbe (Delta-T Devices, Cambridge, UK), and a Li-190 quantum sensor (LI-COR Biosciences, Lincoln, NE, USA), respectively.All measurements were visually inspected before analyses, and only measurements that showed a consistent linear relationship between CO 2 and time for at least 60 s were used.PAR ranged from 493 to 1480 and 211 to 650 μmol photons m À2 s À1 during the two measurements, respectively.Fluxes were calculated as the rate of change in [CO 2 ] over time using the following formula (Jasoni et al., 2005): where δ CO 2 δt is the slope of the CO 2 concentration against time (μmol mol À1 s À1 ), P is the atmospheric pressure (kPa), R is the gas constant (0.008314 kPa m 3 K À1 mol À1 ), T is the air temperature inside the chamber ( C), V is the chamber volume (m 3 ), and A is the ground surface area (m 2 ).
Gross ecosystem productivity (GEP) was calculated as NEE + R eco and standardized to PAR = 700 μmol photons m À2 s À1 using a rectangular hyperbolic relationship (Thornley & Johnson, 1990): where αGEP is the initial slope of the rectangular hyperbola, or apparent quantum yield of GEP, and GEP max is the asymptotic maximum GEP at high light intensities.We used the mean flux rate across measurement days for each plot in the subsequent statistical analyses.
In addition to our static chamber measurements, we calculated July (peak growing season) soil respiration (R soil ) rates based on data from passive, forced-diffusion dynamic chambers (Eosense, Dartmouth, Canada) (Risk et al., 2011), which recorded continuous soil respiration rates at 4-to 6-h intervals in the Dryas heath during the 2015-2017 growing seasons (n = 3 per treatment).For each plot, we used the mean July R soil flux rate from the three measurement years in the statistical analyses.

Data analyses
All analyses were done in R version 4.2.1 (R Core Team, 2021).

Changes in taxonomic composition
To assess the effects of habitat, year, and the warming treatment on plant taxonomic community composition, we used principal component analysis (PCA), as implemented in the rda function in the vegan package (Oksanen et al., 2020).These analyses were based on square-root-transformed abundance data, and we tested whether community composition varied among treatments, habitats, and years (treatment Â site Â year) using the adonis function.For a more targeted test of the effect of experimental warming over time, we conducted PCAs for each habitat separately.This allowed us to zoom in on the treatment-control contrast over time (treatment Â time), which may otherwise be lost in the much larger compositional differences between habitats compared to the treatments.
We calculated the change in Bray-Curtis distances, i.e., species compositional differences between plots, and between the first to last survey from each plot, using the vegdis function in the vegan package (Oksanen et al., 2020).We also calculated change in the summed abundances of all vascular plants; summed abundances for the functional groups, forbs, graminoids, shrubs, bryophytes, and lichens; and species richness (number of species), diversity and species evenness (Shannon diversity/log[richness]) from the first to last surveys.We tested whether change in these community metrics varied by treatment and habitat using ANOVA, where the response was modeled as a function of habitat, treatment, and their interaction.We also tested whether a change in habitat-treatment combinations differed from zero using t-tests.

Changes in functional trait composition and the importance of ITV
To assess inter-and intraspecific trait variability, we calculated the community-weighted mean trait values for each plot using community species composition data from 2015 and trait values from 2018.Since community composition changes over time in this system are slow, the time difference in collecting the last species community data (2015) and trait (2018) data was assumed minor for the trait analyses.We imputed trait data if any species present in a plot did not have at least two measurements for each trait from that plot to be able to generate trait distributions (see following discussion).For this, we first used measurements from that species in the same habitat and treatment, and if that was not available, we used measurements from any plot in the same habitat to fill in the trait coverage for the whole community in each plot.Thus, the majority of trait measurements used in estimating community trait distributions for a plot came from traits measured directly in the plot of interest.For the traits plant height, dry mass, leaf area, thickness, SLA, and LDMC, 90% of the plots had at least 90% coverage at the plot level.For the nutrient-related traits, 75% of the plots had at least 76.2% coverage at the plot level.The community trait distribution was generated using a bootstrapping approach with the R package traitstrap (Telford et al., 2021) by randomly sampling measured trait values with replacement for each trait and each species from each plot proportional to the abundance of that species in that plot (Enquist et al., 2015;Wieczynski et al., 2019).This procedure was repeated 100 times, and for each repetition we calculated the mean of the trait distribution for each trait.These 100 trait distribution means were then averaged to determine each trait's overall distribution mean in each plot.This approach has the distinct advantage of allowing us to include intraspecific variability and the hierarchical structure of the experimental design in our estimates of community trait distributions.For each trait we assessed the differences in trait means between habitat types and treatments in all cases using ANOVAs where the community mean trait values response of interest were modeled as a function of habitat type, warming treatment, and their interaction.
In parallel to the taxonomic compositional analyses, we used PCA, as implemented in the rda function in the vegan package (Oksanen et al., 2020) to assess the effects of habitat and the warming treatment on plant trait community composition.We based these analyses on centered and scaled data and tested whether trait community composition varied between habitats and treatments using the adonis function.For a more targeted test of the effect of experimental warming, PCA was also done for each habitat separately to zoom in on the treatment-control contrast, which could otherwise be lost in the much more significant differences between habitats compared to the treatments.The effect of treatment was tested using the same test as described earlier.
To quantify the relative importance of inter-versus intraspecific variation in the functional community response to warming, we compared the plot-level trait distribution means using the bootstrapping method described earlier, computed with and without the ITV.First, measures including ITV, i.e., based on all trait measurements for each species in each plot, were calculated as described in the preceding section (referred to as the specific mean).Second, measures excluding ITV were estimated by calculating the average trait value for each species across all control plots, and using these values to calculate community-weighted mean trait values for each plot, thereby ignoring plot-, habitat-, and treatment-specific trait variability within species (referred to as the fixed mean).Differences between the specific (including both species turnover [i.e., interspecific variation] and ITV), and fixed (only species turnover) trait means reflect the contribution of ITV to the total observed trait variation across treatments and habitats.To quantify the contribution of ITV and species turnover to the variation in community mean trait values, habitats, and treatments, separate ANOVAs with fixed, specific means and their difference as response and the interaction of habitat and treatment as predictor were performed.We then decomposed the sum of squares from the ANOVAs as described by Lepš et al. (2011).To test whether the difference between specific and fixed community-weighted traits in each treatment from each habitat was different from zero, a t-test was used.We also determined the relative importance of intraspecific variation versus interspecific variation in traits by partitioning trait variance into between-and within-species components using a mixed model where the fixed effect was only an intercept and a random effect for species intercepts based on methods from Messier et al. (2010) using the lme (Bates et al., 2015) and varcomp (Qu, 2017) functions.

Ecosystem CO 2 fluxes
To examine whether taxonomic and functional trait changes lead to shifts in ecosystem carbon fluxes, we tested the effect of experimental warming on GEP 700 , R eco , and NEE using general linear models.We specified treatment and habitat as fixed main effects for all models, including their interaction term and plot as random effect.In addition, we used a general mixed-effects model (lme4 package; Bates et al., 2015) to analyze treatment effects on our multiyear R soil data, with treatment as fixed main effect and plot as random effect to account for our repeated measures design.
We investigated the effects of three groups of variables: (1) environmental (plot-scale microclimate measurements), (2) taxonomic (plant functional group abundance [graminoid, forb, bryophyte, evergreen shrub, deciduous shrub, lichen], evenness, richness, diversity, canopy height), and (3) plant functional traits (community-weighted traits including and not including intraspecific variability: plant height, leaf area, leaf thickness, LDMC, SLA, P, C, N, C:N) on ecosystem carbon flux (i.e., GEP 700 , R eco , and NEE) across habitats and treatments.We assessed the effect of individual predictor variables on ecosystem carbon flux using linear models on scaled predictor variables (Appendix S1: Figure S8).Next, we constructed multiple linear models for all three CO 2 flux components (GEP 700 , R eco , and NEE) across habitats and treatments of the different predictor groups (environment, taxonomic, and plant functional traits), resulting in a total of nine models.Variables were selected based on the stepwise backward selection using the stepAIC function (MASS package; Ripley et al., 2020).Variables were only excluded if they reduced the Akaike information criterion (AIC) score by more than 2 (Burnham & Anderson, 2002).The resulting models represented the overall effects of the three predictor groups on different ecosystem carbon fluxes.To assess the effect of ITV on ecosystem carbon fluxes, we compared the three trait models selected by backward selection including ITV with models containing the same trait predictors but with ITV excluded from the trait predictors.Last, we combined the three models in one full composite model, consisting of combined environmental, taxonomic, and trait variables, for each CO 2 flux component.This allowed us to assess the overall variation explained by each group of variables and shared variation explained across models.

Environmental characteristics and community structure
There was substantial seasonal and annual variation in the measured environmental parameters (Appendix S1: Figure S1A).PAR reached a maximum (>700 μmol m À2 s À1 ) in late June in all measured years.The monthly mean air temperature (2 m above ground) was highest in July and ranged between 6.9 AE 0.26 C (in 2017) and 8.7 AE 0.30 C (in 2016).Temperature fluctuations were much more extensive during winter than summer, with a few warm spells (temperature above 0 C) each winter (Appendix S1: Figure S1A).At the plot level, surface mean annual temperatures in control plots ranged between À1.57AE 0.96 C and À0.03 AE 0.87 C and July temperatures between 9.1 AE 0.54 C and 12.0 AE 0.45 C across years and habitats, with larger variation during winter (Appendix S1: Figure S1B).Mean July soil temperatures were ~3 C lower than at the surface, ranging between 5.8 AE 0.41 C and 8.7 AE 0.04 C across years and habitats (Appendix S1: Figure S1B).Owing to the sporadic data and considerable interannual variation, it was not possible to detect trends in the timing of snowmelt or temperature change in the control plots over the 17-year study period.
There was a trend toward increased annual surface temperature in response to the warming treatment by 1.25 C across all years and habitats (F 1,486 = 3.5, p = .062),but not for soil temperature (0.63 C higher, not significant).The warming treatment enhanced July surface temperatures by, on average, 1.6 C across all years and habitats (F 1,45 = 21.0,p < .001)and tended to increase soil temperature (Appendix S1: Figure S1B).This trend was also reflected in the July 2018 spot measurements of soil temperature and a reduction in soil moisture by the warming treatment (significant only in the Dryas heath, reduction from 15% to 10% volumetric water content; Appendix S1: Figure S2).The warming treatment did not affect surface and soil temperatures in winter.
Canopy height was relatively low but varied among the habitats.The low vascular plant cover in the snowbed resulted in low average plot canopy height above the moss layer (1.1 AE 0.0 cm in 2015, Appendix S1: Figure S3).The canopy height was highest in the Cassiope heath, where the erect growing Cassiope tetragona dwarf shrub dominated (average plot canopy height 1.3 AE 0.1 cm; Appendix S1: Figure S3), while the dominance of the prostrate growing Dryas octopetala dwarf shrub resulted in lower canopy height in the Dryas heath (average canopy height 1.1 AE 0.0 cm; Appendix S1: Figure S3).The warming treatment significantly increased canopy height in the Cassiope and Dryas heaths but not in the snowbed.

Taxonomic compositional change over time and in response to warming
Plant community species composition differed over habitats, treatments, and time, with significant differences in treatment effects and temporal dynamics between the habitats (Figure 1a).Most of the variation in community composition was found between habitats, which could account for about two-thirds of the explained variation in the data (Appendix S1: Table S1).In all habitats, there was a change over time in plant community species composition from the first plant community analysis in 2003 to the last in 2015, both in the ambient controls and treatment plots (Figure 1b-d).Also, there was a significant treatment effect in the Cassiope and Dryas heaths, mainly reflecting somewhat more pronounced community shifts over time in the warmed plots than in the controls (Figure 1c,d).A few other community metrics were affected by warming, most strongly in the Dryas heath, where evenness decreased and total vascular plant and evergreen shrub abundance increased (Appendix S1: Figure S4).

Functional community compositional variation and response to warming
The three communities were also well differentiated in functional trait composition.The snowbed plots occupied the resource-acquisitive part of the trait space, whereas the other two communities were characterized by more resource-conservative trait values and were, in turn, separated by size-related traits, the Cassiope heath being characterized by taller plants than the Dryas heath (Figure 2a; Appendix S1: Table S2).Overall, experimental warming significantly affected community trait composition (Appendix S1: Table S2).Plants tended to be taller, having larger leaves and higher LDMC and carbon content in the warmed plots, and this trend was most pronounced (but not significant) in the Dryas heath (Figure 2b-d).
In univariate analyses, all measured functional traits except LDMC and %P differed between habitats (Appendix S1: Figure S5; Appendix S1: Table S3) and overall, habitat type explained between 13% and 80% of the total variation in individual traits (Appendix S1: Figure S6).In most cases, this was because the snowbed differed from the other two habitats, being characterized by higher values for the traits SLA, dry mass, leaf area, % N, and δ 15 N, along with lower values for leaf thickness, % C, and C:N ratio (Figure 2; Appendix S1: Figure S5).Overall, the warming treatment explained much less of the trait variation than habitat (Appendix S1: Figure S6).Warmed plots had consistently higher leaf dry mass and leaf area and lower δ 15 N than control plots, with the strongest responses in the Dryas heath.The warming treatment did not significantly affect other traits (Appendix S1: Figure S5).

Intraspecific trait variability
ITV accounted for between 30% (plant height) and 71% (%P) of the total trait variation in the whole community across all habitat types (Appendix S1: Table S4).The trait variation attributable to the warming treatment or the interaction between habitat and treatment tended to be due to intraspecific variation (Appendix S1: Figure S6).
In both control and warmed plots, including the intraspecific variation shifted the mean community-weighted values of most traits and habitats (Figure 3).Within the control plots, ITV significantly increased %N, decreased % C, and marginally decreased C:N in snowbed controls, and ITV significantly decreased dry mass, leaf area, and SLA and tended to reduce %N and δ 15 N and increased C:N (p < 0.1) in Dryas heath controls (Figure 3).Within the warming plots, ITV contributed to increased dry mass and leaf area in the snowbed and Dryas heath, and increased SLA, decreased P%, and decreased δ 15 N values in warming plots in the Dryas heath (Figure 3).Overall, we found that the importance of ITV (i.e., the effect of including the intraspecific variation in the trait mean calculations) differed between warmed and control plots for four traits in the Dryas heath (Appendix S1: Table S5).The trait mean increased when ITV was included for dry mass, leaf area, and SLA, whereas for δ 15 N it decreased under warming relative to the controls.None of the other traits or habitat types showed significant differences in the importance of ITV for the trait mean between warming and control treatments.

Habitat and warming effects on ecosystem CO 2 fluxes
Peak growing season GEP 700 and R eco tended to differ between the habitats (Appendix S1: Figure S7).Fluxes measured from the Dryas heath were generally greater than fluxes in the snowbed, whereas fluxes within the Cassiope heath were intermediate between the snowbed and the Dryas heath (Appendix S1: Table S6, Figure S7).In contrast, NEE rates were similar between the Dryas heath and the snowbed.Both habitats were net sources of CO 2 to the atmosphere during our daytime measurements, whereas the Cassiope heath was neither a sink nor a source (i.e., NEE was not different from 0).Across all habitats, experimental warming increased peak growing season R eco and, similarly, the warming treatment enhanced July R soil within the Dryas heath (Appendix S1: Table S6, Figure S7).There was also a significant effect of the warming treatment on GEP 700 , but this effect was driven primarily by a strong positive response within the Cassiope heath (significant treatment Â habitat interaction; Appendix S1: Table S6) and negligible change across treatments within the other two habitats.

Effects of microclimate, taxonomic composition, and functional trait composition on ecosystem CO 2 fluxes
Ecosystem CO 2 fluxes were related, to different degrees, to single variables representing the environmental, taxonomic, and functional trait variable groups.Specifically, whereas GEP 700 was affected by a multitude of different taxonomic and functional trait variables, fewer single variables significantly affected R eco and NEE (Appendix S1: Figure S8).Surprisingly, the fixed and specific (i.e., including ITV) versions of the functional trait variables often differed significantly both in the magnitude and direction of their single trait effects on ecosystem CO 2 fluxes (Appendix S1: Figure S8).Nevertheless, more consistent trait patterns emerged when models were constructed that included the impact of all environmental (microclimate), taxonomic, and plant functional trait variables (Appendix S1: Table S7).For GEP 700 , environmental and taxonomic variables accounted for more than half of the explained variation, with canopy temperature, canopy height, and leaf N variables retained in the final model (Appendix S1: Table S7; Figure 4a,b).For R eco , the variance that can be explained by environmental, taxonomic, or functional trait variables (i.e., shared variance due to the collinearity of these groups of predictors) dominate (Figure 4a,b).Therefore, canopy temperature, plant species diversity, and traits associated with both leaf economics and plant size were retained in the final model (Appendix S1: Table S7).For R eco , environmental variables alone explained more variance than taxonomic or functional trait variables alone, although not by much (Appendix S1: Table S7).In contrast, for NEE, the resource acquisition-related traits, leaf area, and plant height were clearly the most important group of explanatory variables (Appendix S1: Table S7; Figure 4a,b).All three flux components retained some analog of temperature, plant size, and leaf economics variables in their final models.
Including ITV enhanced the explanatory power of traits for all ecosystem CO 2 fluxes, especially for GEP 700 and R eco (Figure 4c; Appendix S1: Table S7).For NEE, including ITV also increased the total variance explained by microclimate, taxonomic community, or functional traits, whereas for GEP 700 and R eco it did not change substantially (Figure 4b vs. a).Negative covariance was relatively minor for the three combined models.

DISCUSSION
Community species composition changed over time in all three habitats.Although species compositional responses to long-term experimental warming were found in two habitats, the effects were minor relative to the overall habitat and temporal variability.Plant community functional trait composition was modestly affected by the warming treatment, detected in three out of 12 traits and mainly in one of the habitats.The trait responses to warming were primarily driven by intraspecific trait variability.Both peak season rates of photosynthesis (GEP 700 ) and respiration (R eco ) increased in response to the warming treatment across habitats.Together, environmental, taxonomic, and functional trait variables explained a large proportion of the total variation in peak season CO 2 fluxes, and this proportion increased when ITV was accounted for.

Plant community change over time in taxonomic space and effects of warming
The change in taxonomic plant community composition in the ambient control plots over the study period reflected both directional response to recent climate warming in Svalbard (Nordli et al., 2020) and dynamic plant community responses to the much more significant climatic interannual variability.Similar to other long-term studies from high Arctic regions (Hollister et al., 2015;Hudson & Henry, 2009;Hudson & Henry, 2010), the directional changes were slow in comparison with dynamic temporal changes (e.g., van der Wal & Stien, 2014) and compared with community changes found in warmer parts of the tundra (Elmendorf, Henry, Hollister, Björk, Boulanger-Lapointe, et al., 2012a).Two resurvey studies in central Svalbard after 70 and 85 years also reported no directional changes or indications of "greening" in response to ongoing climate change (Kapfer & Grytnes, 2017;Prach et al., 2010).
The warming treatment only slightly modified species compositional changes, and these effects were confined to the Cassiope and Dryas heaths, our driest habitats, where the abundance of evergreen shrubs increased.This contrasts with earlier studies that found the strongest responses in moist habitats (Edwards & Henry, 2016;Elmendorf, Henry, Hollister, Björk, Bjorkman, et al., 2012b).Furthermore, the commonly detected decrease in soil moisture in response to the OTC treatment (e.g., Bokhorst et al., 2013) was only found in the Dryas heath, suggesting an environmental factor that overrides soil moisture limitation.In Arctic landscapes, topography determines snow accumulation and hydrology (Walker, 2000) and, thus, the length of the growing season, which is particularly critical in the high Arctic.The early snowmelt in the Dryas heath resulted in a 1.4 times longer growing season, on average, compared to the late-melting, moist snowbed, allowing a longer cumulative response time to the warming treatment, which may explain why we saw the most significant responses within this habitat.
Several factors might contribute to the resistance of the plant community composition to both ambient and long-term experimental warming in the high Arctic.First, high Arctic regions are characterized by small species pools, often substantial dispersal barriers, which may delay the migration of more thermophilic and responsive species, as illustrated by the substantial lagged responses to climate in their subsequent establishment into the plant communities observed throughout the Arctic (Elmendorf et al., 2015).Second, most Arctic plant species are long-lived, and many rely on various forms of clonal growth for population maintenance (e.g., J onsd ottir, 2011), enhancing their ability to persist in the face of environmental change.Third, the plants are adapted to sizeable interannual climate variability, and it may take substantial, long-term warming to push the system beyond that variability range (e.g., Hudson & Henry, 2010;J onsd ottir, 2005).For example, in our study, the 3 C variability range for surface mean July temperatures across study years exceeded the 1.6 C experimental temperature increase.However, the severe damage caused by basal ice formation in some of the Cassiope heath plots during our study indicates that increased frequency of extreme events might accelerate community changes in the high Arctic.Finally, as suggested by our results, plant functional trait variation, particularly ITV, may substantially enhance community functional resilience by enabling plastic responses by individuals (as demonstrated by dynamic community vegetation responses to interannual climate variability) or adaptation to long-term warming, thereby retaining fitness within populations.

Community variation in functional space and effects of warming
The first two axes of the plant community-weighted functional trait space of our results primarily reflected the two global dimensions of leaf economics spectrum and plant size, respectively (Bruelheide et al., 2018;Díaz et al., 2016), albeit with reduced trait ranges, particularly for size traits, comparable to the "extreme" tundra plants in a study by Thomas et al. (2020).Snowbeds were characterized by resource-acquisitive traits, whereas the Dryas and Cassiope heaths were characterized by resource-conservative traits, and were again separated by traits related to size and carbon dynamics.
Community functional trait composition responded to experimental warming and involved increased values of two size-related traits, leaf area and leaf dry mass, and decreased δ 15 N, with the strongest responses in the Dryas heath.This partly contradicts previous studies, where leaf area decreased with warmer temperatures at dry sites but increased with warming at moist sites (Bjorkman et al., 2018).Furthermore, plant height did not respond to warming, again contradicting other tundra studies on individual and community levels (Baruah et al., 2017;Bjorkman et al., 2018;Hollister et al., 2015;Hudson et al., 2011).The inherently low stature of most high Arctic plants may be an adaptation to the harsh climate.Such constraining adaptation in combination with slow within-habitat species and genotype turnover may explain this lack of community plant height response and why size responses to experimental warming were expressed by leaf traits (increased leaf area and dry mass) rather than by increased plant height in this high Arctic community.
The significant decrease in community foliar δ 15 N in response to warming contradicts shoot-level results for five plant species at a high Arctic site in Canada, where no warming responses were found (Hudson et al., 2011).An increase in δ 15 N has been observed in relation to N-availability gradients (Craine et al., 2015) and has also been positively related to mean annual temperature and precipitation (Amundson et al., 2003).It is possible that the decrease in this trait along the snow-moisture gradient of our study, from relatively high in the snowbed to lower levels in the Dryas heath, signals such an N-availability gradient.Experimental warming generally increases N-mineralization rates in cold regions (Salazar et al., 2020), and we would therefore expect an increase in leaf δ 15 N in response to warming.However, the opposite happened in the Dryas heath.The decrease in soil moisture may have played a role.Still, such general trends may also be modulated by other functional traits of the responding species at a local scale, in this case by increased input of recalcitrant litter by the abundant evergreen dwarf shrub Dryas octopetala.Indeed, Vowles and Björk (2019) suggested that, in contrast to the accelerating impact of expanding deciduous shrubs on ecosystem processes in response to climate warming, such as rates of N mineralization (Mekonnen et al., 2021;Myers-Smith et al., 2015), expanding evergreen shrubs will decelerate ecosystem processes through their low nutrient and recalcitrant litter, thereby counteracting the direct effects of warming on those processes (e.g., Cornelissen et al., 2007).
It has been suggested that the frequently observed drop in leaf nutrient levels in response to short-term experimental warming (2-5 years; e.g., Doiron et al., 2014;Tolvanen & Henry, 2001) will level out in the long term as increased soil N mineralization satisfies plant nutrient demand (Michelsen et al., 2012).The lack of leaf nutrient responses to warming in our study is consistent with other long-term warming experiments in the high Arctic (Hudson et al., 2011) and more likely reflects either no direct warming impacts on soil N mineralization or an indirect decelerating impact by increased input of recalcitrant litter, or a combination of these.This is supported by a recent study showing no effects of long-term experimental warming on decomposition rates within our Svalbard communities (Björnsd ottir et al., 2022).

The role of ITV
The contribution of ITV to the total community-level trait variation in our study, ranging between 30% and 71% for different traits, was relatively large compared to both global and tundra averages (Siefert et al., 2015;Thomas et al., 2020).This supports the general expectation that intraspecific variation should be relatively significant in species-poor regions (Siefert et al., 2015) and harsh environments (Niu et al., 2020).
It is usually assumed that a large ITV reflects a greater ability of organisms to exist along broad environmental gradients (large niche breadth) and to adjust to environmental changes (Sides et al., 2014).In all the habitats of our study, the most dominating plant species have relatively large geographic ranges and are apparently adapted to sizeable interannual climate variability (J onsd ottir, 2005).Accordingly, high intraspecific variability in size-related traits may contribute to dynamic annual variation in aboveground plant biomass, as has been observed in Svalbard plant communities in response to interannual variability in summer temperatures (Petit Bon, 2020;van der Wal & Stien, 2014), enhancing plant community resilience to climate variability.
Although the overall community functional trait response to the long-term warming treatment of our study was modest and confined to three of the 12 measured functional traits, it was primarily explained by ITV.This supports our hypothesis that the plant community responses to long-term warming in functional space are mainly driven by intraspecific variation in traits.These results indicate that studies that rely solely on changes in taxonomic composition or community-weighted means that do not incorporate intraspecific variation are likely underestimating the effects of warming.Leaf size traits measured at the individual level have been reported as responsive to warming experiments for a range of Arctic species (Baruah et al., 2017;Hudson et al., 2011), indicating plasticity.A transplant experiment in China revealed a high relative plasticity in leaf δ 15 N among alpine plants (Henn et al., 2018).That was most likely also the case for this trait in the predominantly long-lived clonal plants in response to warming at our Svalbard site.

Effects of warming and functional traits on peak season ecosystem CO 2 fluxes
Apart from the Cassiope heath, the studied habitats appeared to be net sources of atmospheric CO 2 at peak growing season under ambient temperatures, as is often observed in high Arctic regions (Christiansen et al., 2012;Welker et al., 2004).Experimental warming increased GEP 700 and R eco , resulting in negligible overall changes to net ecosystem CO 2 balance, NEE.As there was relatively little warming-induced change in the plant community in either taxonomic or functional trait space, we attribute the largely offsetting warming responses in GEP 700 and R eco to the direct kinetic impact of temperature on plant and soil microbial physiological activities.Although our flux measurements were obtained over just two measurement days, and the limited sample size therefore warrants caution, the observed flux responses are similar to reports from more flux-focused studies in high Arctic Canada and Greenland (Lupascu et al., 2014;Marchand et al., 2004;Welker et al., 2004).Therefore, we believe our flux data adequately reflect the relative importance of environmental, plant taxonomic, and functional trait variables on CO 2 exchange rates, i.e., ecosystem functionality, across habitats during this critical time for annual plant-carbon uptake.
Overall, the largest proportion of the variation in GEP 700 and R eco that was uniquely attributed to only one group of variables, i.e., variation not shared between groups, was accounted for by the environment (canopy temperature).This likely reflects the short response time by the main mechanisms responsible for these CO 2 fluxes.Nevertheless, variables belonging to the taxonomic and functional trait groups also explained considerable variation, either alone or as the shared variance between two or more variable groups.For example, canopy height or the size-related traits of plant height and leaf area, as well as the leaf economics trait leaf N content (Wright et al., 2004), consistently remained in the final, reduced models for GEP 700 and R eco , reflecting the importance of these traits for ecosystem CO 2 exchange (Reich et al., 1997).
Overall, all three variable groups accounted for broadly similar amounts of variation in GEP 700 and R eco rates (R 2 range across individual groups alone = 0.2-0.37),although with functional traits situated at the lower end (functional trait group R 2 = 0.2-0.25,when including ITV).In contrast, plant functional traits uniquely explained a large proportion of the overall variation in NEE (>40%), primarily driven by productivity-related traits, plant height, leaf area, and leaf N. The explanatory power of our combined group models was moderate to high (R 2 range across the three flux components, GEP 700 , R eco , and NEE: 0.4-0.55),which is similar to or greater than models from two recent subarctic tundra studies also attempting to link traits to ecosystem CO 2 exchange rates (Happonen et al., 2022;Sørensen et al., 2018).Although studies that relate functional community composition directly to measures of ecosystem function have generally been lacking (Funk et al., 2017), an increasing number of studies report strong effects of community-scaled traits on ecosystem properties (Grigulis et al., 2013;Reich et al., 2014).However, even if not all ecosystem properties are explained well by functional community composition (van der Plas et al., 2020), plant-related ecosystem functions, such as carbon-cycling processes, consistently show at least moderate and often good relationships (Funk et al., 2017;Happonen et al., 2022;Reich et al., 2014;van der Plas et al., 2020;this study).In addition, for all three CO 2 flux components, we consistently found that a more significant proportion of the variance explained by the functional traits variable group alone was accounted for when ITV was included.Thus, variability in size and leaf economics traits-i.e., the functional traits that constrain the two primary dimensions along which most individual and community trait assembly occurs (Bruelheide et al., 2018;Díaz et al., 2016)-is an essential component in determining community functionality at our high Arctic site.Similarly, Happonen et al. (2022) found positive relationships between subarctic tundra CO 2 fluxes and within-community variability in SLA and LDMC (but not for plant height), suggesting that increasing functional diversity is linked to ecosystem functionality (Cadotte et al., 2011).Adding to this, we found that the total explained variance in our final, reduced NEE model increased when including ITV, indicating an effect of trait plasticity on ecosystem carbon balance.Taken together, these relatively straightforward effects of including trait variation within species fit well with recent studies showing that intraspecific variation in traits repeatedly impacts ecosystem functionality at least as much as species turnover effects (e.g., Des Roches et al., 2018).
Although our analyses of the ecosystem CO 2 fluxes shed new light on how plant functional traits and especially ITV can affect ecosystem functioning, we also note that ecosystem CO 2 fluxes are highly variable in time and space.Consequently, more extensive flux measurements are required to capture the ecosystem characteristics.Nevertheless, we show similar linkages, and explanatory power, between tundra plant functional composition and ecosystem CO 2 fluxes as other studies with more flux measurements (Happonen et al., 2022;Sørensen et al., 2019), and so our study adds to the evidence that trait-functionality relationships, at least regarding CO 2 exchange rates, exist and are consistent across tundra habitats and sites.

CONCLUSIONS
Our study provides new insights into the impacts of climate warming on plant communities across high Arctic tundra landscapes and the role of plant functional trait variation on community resilience and ecosystem functioning.The studied plant communities were strongly differentiated among habitats in taxonomic and trait composition and showed a modest but significant functional response to experimental warming.Consequently, functional trait variation explained considerable variation in ecosystem CO 2 exchange rates, broadly similar in magnitude to environmental and taxonomic variables.However, warming responses in CO 2 fluxes were mainly driven by the direct kinetic impacts of temperature on plant physiology and biochemical processes.The results also provide evidence that ITV adds to community functional resilience, which plays a vital role in vegetation resistance to climate warming in high Arctic Svalbard.
Here we have focused on climate warming impacts in relation to habitat differences in snow regimes.It is likely, however, that, combined with other aspects of predicted climate change, such as changes in precipitation and snow accumulation and increased frequency of extreme events, for instance, rain on snow, climate change impacts will eventually exceed community resilience in some habitats and cause substantial community shifts in the future with consequences for ecosystem functioning.Finally, as climate warms, the immigration of more thermophilic and responsive plant species may eventually facilitate faster community change.

F
I G U R E 1 (a) First two axes of principal component analysis (PCA) of plant community species composition in control plots (open symbols) and experimentally warmed plots (closed symbols) in the three studied habitats (snowbed, Casiope heath, Dryas heath) over the years 2003 (large symbols), 2009, and 2015, with the amount of variation explained by each PCA axis indicated.The species with the best fit (length of vector) are visualized.(b-d) First two PCA axes for community species composition change of plots within each habitat.
First two axes of principal component analysis (PCA) of plot-scale plant community functional trait composition (i.e., community-weighted trait means) of plots in different habitats (snowbed, Cassiope heath, Dryas heath) and treatments (control, warming) with the amount of variation explained by each PCA axis indicated.Colored points indicate habitat type; point shape indicates treatment type.Vectors indicate how traits are related to the first two axes and are labeled by the trait they represent.(b-d) First two PCA axes for functional composition of plots within each habitat.

F
I G U R E 3 Contribution of intraspecific trait variation (ITV) to trait variation calculated as difference between mean community functional trait distributions with (specific average) and without (fixed average) ITV in different habitats and treatments (control, warming).The zero line indicates that ITV did not change the community trait distribution mean value, whereas positive values indicate that ITV increased the community trait distribution mean, and negative values indicate that ITV decreased the community trait distribution mean.Asterisks above each habitat and treatment combination indicate whether the mean difference is different from zero (t-test).+ p < 0.1; *p < 0.05.

F
I G U R E 4 Variance of ecosystem function, i.e., peak growing season GEP 700 , R eco , and net ecosystem exchange flux rates, across habitats and treatments explained by three groups of variables: (1) environmental (plot-scale microclimate measurements), (2) taxonomic (plant functional group abundance [graminoid, forb, bryophyte, evergreen shrub, deciduous shrub, lichen], evenness, richness, diversity, canopy height), and (3) plant functional traits (community weighted traits with and without including intraspecific trait variability: plant height, leaf area, leaf thickness, leaf dry matter content, specific leaf area, P, C, N, C:N).(a) Variance explained by unique and combined group effects excluding intraspecific trait variation (ITV).(b) Variance explained by unique and combined group effects including ITV.(c) Variance explained by plant functional traits with and without ITV (gray + black and gray bars, respectively).Note that negative values can occur when combining two groups in a model that explains less of the variation than the individual groups.