Plant trait response of tundra shrubs to permafrost thaw and nutrient addition

Plant traits reflect growth strategies and trade-offs in response to environmental conditions. Because of climate warming, plant traits might change, altering ecosystem functions and vegetation–climate interactions. Despite important feedbacks of plant trait changes in tundra ecosystems with regional climate, with a key role for shrubs, information on responses of shrub functional traits is limited.

Here, we investigate the effects of experimentally increased permafrost thaw depth and (possibly thaw-associated) soil nutrient availability on plant functional traits and strategies of Arctic shrubs in northeastern Siberia. We hypothesize that shrubs will generally shift their strategy from efficient conservation to faster acquisition of resources through adaptation of leaf and stem traits in a coordinated whole-plant fashion. To test this hypothesis, we ran a 4 year permafrost thaw and nutrient fertilization experiment with a fully factorial block design and six treatment combinations – permafrost thaw (control, unheated cable, heated cable) × fertilization (no nutrient addition, nutrient addition). We measured 10 leaf and stem traits related to growth, defence and the resource economics spectrum in four shrub species (Betula nana, Salix pulchra, Ledum palustre and Vaccinium vitis-idaea), which were sampled in the experimental plots. The plant trait data were statistically analysed using linear mixed-effect models and principal component analysis (PCA).

The response to increased permafrost thaw was not significant for most shrub traits. However, all shrubs responded to the fertilization treatment, despite decreased thaw depth and soil temperature in fertilized plots. Shrubs tended to grow taller but did not increase their stem density or bark thickness. We found a similar coordinated trait response for all four species at leaf and plant level; i.e. they shifted from a conservative towards a more acquisitive resource economy strategy upon fertilization. In accordance, results point towards a lower investment into defence mechanisms, and hence increased shrub vulnerability to herbivory and climate extremes.

Compared to biomass and height only, detailed data involving individual plant organ traits such as leaf area and nutrient contents or stem water content can contribute to a better mechanistic understanding of feedbacks between shrub growth strategies, permafrost thaw and carbon and energy fluxes. In combination with observational data, these experimental tundra trait data allow for a more realistic representation of tundra shrubs in dynamic vegetation models and robust prediction of ecosystem functions and related climate–vegetation–permafrost feedbacks.

1 Introduction

Plants have different strategies to use resources to grow, reproduce, compete with neighbour plants and defend themselves against pathogens and herbivores (Bazzaz et al., 1987; Ordoñez et al., 2010). However, as resources – nutrients, water and light – can have limited availability, plants are subject to trade-offs so that they have to allocate the resources to one function vs. another (Grime, 1977; Westoby et al., 2002; Reich, 2014). Environmental changes, such as the ones promoted by climate warming (e.g. increasing amount of resources available in the soil), may modify these trade-offs and plant strategies (Grime, 2006; Ordoñez et al., 2010).

Plant strategies and trade-offs can be identified by measuring plant traits and their correlations (Grime et al., 1997; Westoby et al., 2002). The so-called leaf economics spectrum is an example of how leaf traits show similar covariation across species according to the resource acquisition strategy, which goes from slow (i.e. resource conservation) to rapid resource acquisition (Díaz et al., 2016). Plant traits also determine plant responses to environmental factors and underpin their effects on ecosystem processes and services (Lavorel and Garnier, 2002; Kattge et al., 2011; Soudzilovskaia et al., 2013). Analysing plant trait responses to climate warming can provide insight into future ecosystem structure and functioning (Díaz et al., 2007).

In low-resource environments such as the Arctic tundra, plants adopt a conservative strategy with low rates of resource acquisition, growth and tissue turnover (Chapin, 1980). Low stature, small leaves of low specific leaf area (thick leaves with dense tissue) and long leaf life span reflect that conservative strategy (Reich et al., 1997; Cornelissen, 1999). This strategy allows plants to allocate resources to other processes, such as defence against pathogens and herbivores, and confer stress resistance (Chapin et al., 1993). Plants, however, might adopt a more acquisitive strategy under the environmental conditions projected for the Arctic in the course of this century (Post et al., 2019).

Surface air temperature in the Arctic has risen more rapidly than in other regions over the past decades and is projected to keep increasing: ∼ 3 ^∘C by the year 2100 under emission scenario RCP4.5 (IPCC, 2013). Ground heat flux and soil temperature are also expected to increase, accelerating permafrost thaw and soil organic matter mineralization (Rustad et al., 2001; Richter-Menge and Overland, 2010; Elmendorf et al., 2012). The release of nutrients trapped in the permafrost (Kokelj and Burn, 2003; Weintraub and Schimel, 2003; Schuur et al., 2009), together with an enhanced soil mineralization rate (Schmidt et al., 1999; Walther et al., 2002), will increase nutrient availability for tundra vegetation (Keuper et al., 2012). Several warming experiments (Elmendorf et al., 2012), satellite imagery (i.e. AVHRR, MODIS and Landsat multi-decadal records of the normalized difference vegetation index (NDVI); Myers-Smith et al., 2011) and repeat multi-decadal aerial photography (Tape et al., 2012) have shown effects of recent climate warming on tundra vegetation growth, productivity and distribution, especially on shrubs (Myers-Smith et al., 2015; Myers-Smith and Hik, 2018). Our current knowledge of tundra shrub responses to climate change concerns mainly their performance traits (detailed in Violle et al., 2007), especially plant height and biomass. We still know precious little about the functional traits underpinning these responses or the effects expanding shrubs may have on ecosystem functions (but see Hudson et al., 2011; Kremers et al., 2015; Barrett and Hollister, 2016). A recent pan-Arctic plant trait sampling and analysis effort has revealed a generally strong spatial temperature–trait relationship, which was, however, mediated by soil moisture (Bjorkman et al., 2018a, b). This study also highlighted the limitations of the observational space-for-time substitution method and identified the need for experimental studies to elucidate intraspecific trait responses to environmental drivers. Shrub responses to climate may have consequences for the carbon cycle (e.g. increase carbon uptake) and the surface energy budget (e.g. decrease albedo), which in turn may affect the regional climate (Eugster et al., 2000; Chapin, 2003; Beringer et al., 2005; Bonfils et al., 2012; Pearson et al., 2013; Juszak et al., 2017). A better understanding of shrub trait responses to climate and of shrub–climate interactions is fundamental to improving dynamic global vegetation models and predictions of vegetation shifts (Cramer et al., 2001; Doherty et al., 2010; Wullschleger et al., 2014).

The objective of this study is to experimentally investigate the consequences of increased permafrost thaw and nutrient addition on aboveground traits and trait coordination of tundra shrubs. We hypothesize that, under simulated future environmental conditions (i.e. permafrost thaw and soil nutrient increase), (i) shrubs will shift their strategy from efficient conservation to faster acquisition of resources through adaptation of leaf and stem traits and (ii) leaf traits, stem traits and plant height will show a coordinated response to these environmental changes as they all belong to the same overall resource economy dimension within the functional trait space. To test our hypotheses, we ran a permafrost thaw and fertilization experiment for 4 years in Siberia and measured 10 plant traits related to the leaf economics spectrum, growth and defence in tundra shrubs. Whereas most previous studies focused on one or two shrub species only, we explicitly compare the responses of four predominant species in order to find commonalities vs. idiosyncrasies of intra- and interspecific trait response, as these are critical for upscaling from site level to tundra ecosystems at larger scales.

2 Materials and methods

2.1 Study area

The study area is located in the nature reserve of Kytalyk, in the continuous permafrost region of Yakutia, northeastern Siberia (70^∘49^′ N, 147^∘28^′ E, 10 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$). Ice-rich permafrost and shallow active layers characterize the area (van Huissteden et al., 2005; Iwahana et al., 2014). The mean annual precipitation is 210 mm and the mean annual air temperature −13.1 ^∘C, with minimum and maximum monthly means of −33.5 ^∘C in January and 11.2 ^∘C in July (1980–2013, WMO station 21946, Chokurdakh, monthly summaries of GHCN-D, NOAA National Climatic Data Center).

The experimental plots were placed on a moist acidic tussock tundra area, the soil of which is classified as Gelisol (Wang et al., 2017). In the Circumpolar Arctic Vegetation Map (Raynolds et al., 2019) the vegetation type in this area is classified as tussock-sedge, dwarf-shrub, moss tundra. The main vegetation has a maximum canopy height of 25 cm and comprises sedge allies (mainly Eriophorum vaginatum), abundant deciduous and evergreen dwarf shrubs, bryophytes and lichens. The growing season lasts from the end of June to the end of August in the study area (Parmentier et al., 2011). The slightly acidic soil (pH 6) has a silty-clay texture and high organic matter content (Blok et al., 2010; Bartholomeus et al., 2012). The soil organic matter decomposition is low as indicated by the high average carbon-to-nitrogen ratio (22) and low cellulose-to-lignin ratio (2.4) (Iturrate-Garcia et al., 2016). During the mid-growing season, the mean active layer thickness is 35 cm, increasing to about 50 cm at the end of the season.

2.2 Experimental design

To test whether climate change might have effects on shrub traits, we ran a permafrost thaw and nutrient fertilization experiment from 2011 to 2014 (Wang et al., 2017). The experiment had a fully factorial block design with five blocks, each with six plots of 1.5 m × 1.5 m placed at randomly chosen locations in July 2010. Each block covered an area of approximately 10 m × 10 m. A picture of a typical block setup is contained in the supplementary material of Wang et al. (2017). Six treatment combinations – permafrost thaw (3 levels) × fertilization (2 levels) – were randomly assigned to the plots within blocks.

The permafrost thaw treatment consisted of no cable, unheated cable and heated cable. For this treatment, we buried heating cables at approximately 15 cm depth in the unheated and heated plots in July 2010 (Wang et al., 2017). One year later, we connected the cables of the heated plots to solar panels in order to increase the thaw depth. The unheated cable plots served as a reference for the permafrost thaw treatment, while plots without cable were included to quantify possible disturbance effects of the cable alone. For the fertilization treatment (nutrient addition vs. no addition), we applied slow-release NPK fertilizer tablets with micronutrients (Osmocote Exact Tablet, Scotts International, Heerlen, the Netherlands). The tablets were applied at approximately 5 cm depth at the start of the experiment (July 2011) and again in 2013 (5.6 g N, 1.4 g P and 3.7 $\mathrm{g}\mathrm{K}\cdot {\mathrm{m}}^{-\mathrm{2}}\cdot {\mathrm{yr}}^{-\mathrm{1}}$), which increased the exchangeable nutrient content mainly in the upper soil layer (Wang et al., 2017).

2.3 Soil temperature and thaw depth

Soil temperature of each plot was measured continuously in 2013/14 at four depths (0, 5, 15 and 25 cm) using temperature loggers (iButton DS1922L/DS1921G, Maxim Integrated, USA). Thaw depth was measured twice in July 2014 by vertically introducing a centimetre-scale metal rod until hitting the frozen soil (Wang et al., 2017).

2.4 Study species and sampling

We investigated the response of four shrub species dominant at the study site and present in all experimental plots: two deciduous species, Betula nana ssp. exilis (Sukaczev) Hultén and Salix pulchra Cham., and two evergreen species: Ledum palustre ssp. decumbens (Aiton) Hultén and Vaccinium vitis-idaea L. (Fig. 1). The abundance of the four species was broadly similar in all plots, except for S. pulchra, which was less abundant.

Figure 1 Study species: Betula nana (a), Salix pulchra (b), Ledum palustre (c) and Vaccinium vitis-idaea (d).

In the mid-growing season of the last year of the experiment (31 July–12 August 2014), we randomly selected six healthy-looking individuals (with less than 20 % leaf damage) of each species per plot, except for S. pulchra, for which only one to four individuals were present per plot. We cut the selected individuals 4 cm below the root collar after measuring their height. The sampling and transport of the plant samples followed the protocol for standardized trait measurements described in Pérez-Harguindeguy et al. (2013). Most of the plant traits were measured in the laboratory within a few hours.

2.5 Plant traits

We selected 10 aboveground plant traits which are related to and provide insight into shrub growth, defence and nutrient acquisition strategies, as well as into the interactions between tundra shrubs and carbon and energy fluxes. We measured the selected leaf and stem traits in each individual of the four shrub species one single time (based on sampled shrub individuals; see Sect. 2.4).

2.5.1 Height

Plant height was measured in the field as the vertical distance from the ground to the tallest vegetative tissue of the selected individuals (maximum vegetative height).

2.5.2 Leaf area (LA) and specific leaf area (SLA)

We cut two leaves per individual, including the petiole, from the top and bottom canopy layers. We scanned the leaves with a flatbed scanner (LiDE 70 Canon Inc., Japan, 300 dpi image resolution) calibrated with a 1 cm² reference. Then, we estimated LA by counting pixels using the software MatLab R2014a (The MathWorks, Inc., MA, USA). We oven-dried the scanned leaves (60 ^∘C, 72 h) and weighed them to determine SLA by dividing the LA of each leaf by its dry weight.

2.5.3 Leaf dry matter content (LDMC)

We followed a variation of the partial rehydration method to determine LDMC using the same leaves as for LA (Vendramini et al., 2002; Vaieretti et al., 2007). To assure maximum hydration, we cut whole individuals in the morning, wrapped the samples in moist paper and put them in sealed plastic bags (Pérez-Harguindeguy et al., 2013). We kept the samples in the dark at low temperatures until they were weighed within the following 6 h to obtain fresh mass. The individual leaves were re-weighed after oven-drying them (60 ^∘C, 72 h). LDMC was the dry mass of a leaf divided by its fresh mass.

2.5.4 Leaf nitrogen concentration (LNC)

Oven-dried leaves were milled, and leaf carbon and nitrogen concentrations were determined by dry combustion (TruSpec Micro-CHN analyser, Leco Corporation, MI, USA) in samples of 2 mg. Then, the carbon-to-nitrogen ratio (C : N) was calculated.

2.5.5 Leaf phosphorus concentration (LPC)

We used a colorimetric assay employing ammonium heptamolybdate to determine LPC. Milled samples of 0.05 g were combusted in a muffle furnace (B180 Nabertherm, Germany) programmed with 1 h heating ramp-up to 600 ^∘C and 2.5 h at 600 ^∘C. We added 2 mL of 0.1 M H₂SO₄ to the ashes, followed by 5 mL of distilled water, and filtered the suspension (Macherey Nagel MN615). The phosphorus in the extracts was determined using a continuous-flow analyser (Skalar Analytical B.V., the Netherlands) calibrated with KH₂PO₄ standards.

2.5.6 Stem-specific density (SSD)

We cut approximately 3 cm long sections of the main stem at one-third of the stem length and removed the bark. We measured the diameter and length of the stem sections without bark, oven-dried (60 ^∘C, 72 h) and weighed them. SSD was determined by dividing the dry mass of a section by its volume.

2.5.7 Stem water content (SWC)

We weighed the sections used for SSD before and after oven-drying them. SWC was estimated as the difference between fresh and dry weight divided by the dry weight.

2.5.8 Xylem diameter and bark thickness

Samples including the 2 cm above and below the root collar of the main stem were taken and preserved in ethanol (40 % vol. aqueous solution) until laboratory processing. We cut thin sections of 20–30 µm along the root collar of each individual and placed them on microscope slides. We photographed and measured xylem diameter and bark thickness following the protocol described in Iturrate-Garcia et al. (2017).

2.6 Statistical analysis

To test if soil temperature and thaw depth were affected by permafrost thaw and fertilization treatments, we used linear mixed-effect models fitted in ASReml (ASReml 3.0, VSN International Ltd., UK). The fixed terms of the models were block (factor with five levels), permafrost thaw treatment (two levels: heating, no heating), fertilization treatment (two levels) and treatment interaction. For the analysis, we averaged the thaw depth values per plot. The soil temperature values were averaged by growing and no-growing season per plot and depth class. Because we aggregated the data per experimental plot, we removed the random term (plot) from the analysis, in order to avoid model overfitting.

Table 1 Average soil temperature (standard deviation) and average thaw depth (standard deviation) of the five experimental blocks for plots grouped into no-fertilization (NF) and fertilization treatments (F). Data are grouped by growing season (T_Jun–Aug) and no-growing season (T_Sep–May) and depth (0, 5, 15 and 25 cm). Significant effects are in bold.

^a P < 0.05. ^b P < 0.01. ^c P < 0.001.

We also used linear mixed-effect models to test the treatment effect on plant traits. Height, LA, bark thickness and xylem diameter were log-transformed prior the statistical analysis to meet assumptions of linearity. First, we analysed plant traits of the four species together and then plant trait of each functional type (PFT; deciduous and evergreen). In the species analysis, we modelled each plant trait as a function of block (a fixed factor with five levels), permafrost thaw treatment (fixed factor with three levels), fertilization treatment (fixed factor with two levels), species (fixed factor with four levels) and the interaction between treatments and species. In addition to these fixed terms, we also considered the interaction between species and block, which was a term recognized in the course of the statistical analysis to take into account species–specific trait differences among blocks. The random terms of the model were plot (factor with 30 levels) and the interaction of plot and species. In the PFT analysis, we modelled plant traits as a function of block, permafrost thaw treatment, fertilization treatment, PFT (a fixed factor with two levels), species and the interaction between plot and species. In both cases, we assessed if the effects of the permafrost thaw treatment on plant traits were due to the disturbance of the buried cables or the treatment per se. For that purpose, we split the three-level permafrost thaw factor into two contrasts of one degree of freedom each, i.e. cable presence (heated and unheated cables vs. no cable) and heating (heating cables vs. unheated cable and no cable). We used the first contrast to assess disturbance effects (heating followed by cable presence) and the second contrast to assess treatment effects (cable presence followed by heating). After running these models for species and PFTs, we found that plant traits were significantly different among species, even between species within the same PFT. Consequently, we analysed the four species separately to maintain ecological information. In this case, we fitted block, permafrost thaw treatment, fertilization treatment and the interaction between the permafrost thaw and fertilization treatments as fixed terms and plot as a random term.

In order to explore shrub plant strategy and its change with treatments, standardized (Z-scored) plant trait data were subjected to a principal component analysis (PCA; vegan package version 2.4-0; Oksanen et al., 2016). We only considered the fertilization treatment (nutrient addition and no addition) in the PCA, as most traits were not responsive to the permafrost thaw treatment (see Results). We performed a separate analysis for leaf traits (SLA, LDMC, LNC, LPC and C : N) and one for stem traits (SSD, bark thickness, xylem diameter and SWC) and height. Scores and variable loadings resulting from the PCA were scaled for visual depiction of data.

To test for relationships between leaf economics, and stem traits and height, we used linear mixed-effect models. We extracted the loadings of the first principal component axes (PC1) of the leaf trait and the stem trait–height PCA. The response variables in our models were height and stem trait–height PC1 loadings. Block and leaf trait PC1 loadings were set as fixed terms, and plot as a random effect. The significance of the linear relationships between variables was analysed using Pearson's correlation coefficients in addition to the linear mixed-effect models.

All data were analysed using R 3.4.1 (http://r-project.org, last access: September 2017).

3 Results

3.1 Treatment effects on soil temperature and thaw depth

Permafrost thaw and fertilization treatments affected soil temperature and thaw depth (Table 1). The soil temperature during the growing season was significantly higher in heated plots than in unheated plots. The soil temperature at 5 cm depth was 0.6 ^∘C higher (F_1,29 = 9.00, p < 0.05), at 15 cm 1.1 ^∘C higher (F_1,29 = 17.97, p < 0.001) and at 25 cm 0.8 ^∘C higher (F_1,29 = 14.02, p < 0.01). The difference of surface soil temperature (0 cm) between heated and unheated plots was not significant. During the no-growing season, temperature differences were significant at all the depths. The soil temperature at 0 cm was 1.4 ^∘C higher in heated plots than in unheated plots (F_1,29 = 16.2, p < 0.01), at 5 cm 1.3 ^∘C higher (F_1,29 = 26.1, p < 0.001), at 15 cm 1.3 ^∘C higher (F_1,29 = 17.5, p < 0.01) and at 25 cm 1.2 ^∘C higher (F_1,29 = 16.9, p < 0.01). The difference in soil temperature between fertilized and unfertilized plots was also significant at all the depths, but only during the growing season. The soil temperature was lower in the fertilized plots: at 0 cm 0.9 ^∘C lower (F_1,29 = 11.6, p < 0.01), at 5 cm 1.0 ^∘C lower (F_1,29 = 19.1, p < 0.01), at 15 cm 0.4 ^∘C lower (F_1,29 = 12.7, p < 0.01) and at 25 cm lower 0.5 ^∘C (F_1,29 = 6.24, p < 0.05). The thaw depth was 10.7 cm deeper in heated plots than in unheated plots (F_1,29 = 24.6, p < 0.001), but 3.9 cm shallower in fertilized plots than in unfertilized plots (F_1,29 = 5.40, p < 0.05). Fertilization treatment effects on soil temperature and thaw depth did not depend on the permafrost thaw treatment.

3.2 Treatment effects on leaf and stem traits and plant height

The permafrost thaw treatment had no significant effect on most shrub traits. Only LA responded significantly to the permafrost thaw treatment (F_1,28 = 18, p < 0.001) when analysing all four species together. At the species level, the permafrost thaw treatment affected only LA of S. pulchra and L. palustre. Individuals of both species had greater LA on heated plots than on control and unheated plots (Table 2). The permafrost thaw treatment only increased SWC for S. pulchra (F_1,28 = 12.8, p < 0.01). Neither the effect of the combination of treatments (permafrost thaw × fertilization) nor the disturbance caused by the buried cables was significant for most measured leaf and stem traits. Exceptions were a significant treatment combination effect on bark thickness of B. nana (F_1,25 = 4.54, p < 0.05) and L. palustre (F_1,25 = 8.15, p < 0.01), and LA of S. pulchra being negatively affected by the buried cables (Table 2).

Table 2 Effects of fertilization (Fert), cable disturbance (Ca) and soil heating (H) on leaf traits of each shrub species (B. nana (Betn), S. pulchra (Salp), L. palustre (Ledp) and V. vitis-idaea (Vacv)). Treatment columns show the average and standard error of the response variables for no fertilization (NFert), fertilization (Fert), no heating cables (Ct), unheated cables (Ca) and heated cables (H). LMM columns show the Wald test outputs for our linear mixed-effect models. Significant effects are in bold. Treatment combination effect (heating × nutrient addition) was not included as it was significant only for LA of Betn (F_1,25 = 15.1, P < 0.01).

^a P < 0.001. ^b P < 0.01. ^c P < 0.05.

The fertilization treatment had a significant effect on all leaf traits, height and SWC, but not on bark thickness, xylem diameter or SSD, when the four species were analysed together (results not shown). At the PFT level, traits were significantly different between deciduous and evergreen species, except for LDMC and SSD. We also found that the fertilization effects on LA, LNC, LPC, C : N and SWC differed between PFTs (Table S2 in the Supplement). The relative increase of LA and decrease of C : N with fertilization was greater for evergreen than for deciduous species. For LNC, LPC and SWC, the increase was greater for deciduous than for evergreen species. At the species level, the fertilization effect on LA, LNC and C : N was significant for all four shrub species (Table 2). Fertilization effects were also significant for SLA, LPC and LDMC of all species except for S. pulchra (Table 2). Leaves in the fertilized plots were bigger and thinner (higher SLA), and had higher nutrient concentration (LNC, LPC) and lower LDMC and C : N than leaves in unfertilized plots (Table 2). For stem traits, the fertilization treatment significantly increased the SSD of B. nana (F_1,29 = 10.1, p < 0.01) and SWC of both deciduous species (B. nana: F_1,29 = 17.8, p < 0.001; S. pulchra: F_1,29 = 13.9, p < 0.01). Xylem diameter and bark thickness responses to nutrient addition were not significant.

3.3 Coordinated trait response to fertilization

In the leaf trait PCA with all four species combined, shrub individuals were separated into species with low overlap along the first principal component axis (PC1) (Fig. 2). PC1 explained 64 % of the variation among individuals and was mainly related to leaf nutrient content (LNC, LPC) and C : N. We found B. nana at the lower end of PC1, associated with high SLA and leaf nutrient concentrations, and V. vitis-idaea at the upper end of the axis. B. nana was the species with the widest range along PC1. The second PC axis (PC2) explained 19 % of the variation and was mainly related to LDMC. Under nutrient addition, we observed a similar trait change for all four species. Leaves on fertilized plots had lower LDMC and C : N and higher LNC, LPC and SLA than leaves on unfertilized plots (Fig. 2).

Figure 2 Principal component biplot of leaf traits for all four shrub species combined, showing change along the leaf conservative–acquisitive continuum (thick black arrows) when nutrients were added. Leaf traits included are leaf dry matter content (LDMC), carbon-to-nitrogen ratio (C : N), leaf nitrogen content (LNC), leaf phosphorus content (LPC) and specific leaf area (SLA). Points are the trait scores of individuals without fertilization (closed circles) and with fertilization (open circles). Sample scores are scaled by a factor of 15, and variable loadings by a factor of 7. Squares indicate the centre of the ordiellipses (standard error with 95 % confidence interval) of trait scores without nutrient addition (solid lines) and with nutrient addition (dashed lines). The first principal component explains 64 % of the total variance, while the second component explains 19 %.

Similar leaf trait space occupation was found when we ran the PCA for each species separately (Fig. 3). PC1 explained a slightly greater amount of total variance among individuals for the evergreen species (65 % for L. palustre and 60 % for V. vitis-idaea) than for the deciduous species (54 % for B. nana and 41 % for S. pulchra). Individuals were separated into two clusters along PC1 corresponding to individuals from fertilized and unfertilized plots. PC2 explained 17 and 18 % of the variation among individuals for L. palustre and V. vitis-idaea, respectively, and 20 % for deciduous species. The main results were maintained when we excluded LNC from the above analysis, showing that PC1 was not driven by the potential correlation of C : N and LNC (Table S3 in the Supplement).

Figure 3 Principal component biplots of leaf traits for each shrub species. Change in leaf traits when nutrients are added is shown by thick black arrows. Leaf traits included in the PCA are leaf dry matter content (LDMC), carbon-to-nitrogen ratio (C : N), leaf nitrogen content (LNC), leaf phosphorus content (LPC) and specific leaf area (SLA). Points are the trait scores of individuals without fertilization (closed circles) and with fertilization (open circles). Sample scores are scaled by a factor of 15, and variable loadings by a factor of 7. Squares indicate the centre of the ordiellipses (standard error with 95 % confidence interval) of the trait scores without nutrient addition (solid lines) and with nutrient addition (dashed lines). The total variance explained by the two first principal components (PC1, PC2) is indicated as a percentage between brackets on the axes.

Similarly, we ran a PCA for stem traits and plant height for each of the four species (Fig. 4). For these traits, individuals overlapped more on the PCA ordination plane. However, there was a trend towards taller individuals having lower SSD and higher SWC in the fertilized plots for three species, but not for V. vitis-idaea. Indeed the stem trait–height space was generally similar for all the species except V. vitis-idaea. PC1 explained slightly more variation among individuals than PC2, especially for deciduous species.

Figure 4 Principal component biplots of plant height and stem traits for each shrub species. Stem traits included are xylem diameter, bark thickness, stem water content (SWC) and stem-specific density (SSD). Change in traits when nutrients are added is shown by thick black arrows. Points are the trait scores of individuals without fertilization (closed circles) and with fertilization (open circles). Sample scores are scaled by a factor of 15, and variable loadings by a factor of 7. Squares indicate the centre of the ordiellipses (standard error with 95 % confidence interval) of the trait scores without nutrient addition (solid lines) and with nutrient addition (dashed lines). The total variance explained by the two first principal components (PC1, PC2) is indicated as a percentage between brackets on the axes.

3.4 Plant strategies – correlation of leaf traits with stem traits and plant height

We found significant correlation between PC1 of the leaf trait PCA (leaf PC1) and plant height for all species, except for S. pulchra (Fig. 5). We also found a significant correlation between leaf PC1 and stem trait–height PC1 for B. nana and V. vitis-idaea (Fig. S1 in the Supplement). Individuals found in the upper range of the stem trait–height PC1 (high values for height, xylem diameter and bark thickness) were also found on the upper extreme of leaf PC1 (high values of LNC and LPC).

Figure 5 Relationship between plant height and the first axis of the principal component analysis (PC1) of leaf traits for Betula nana (a), Salix pulchra (b), Ledum palustre (c) and Vaccinium vitis-idaea (d). Points are trait values for individuals on unfertilized (black) and fertilized plots (grey). Solid lines are values predicted by the linear mixed-effect model, and dashed lines are the upper and lower limits of the predicted value confidence interval. Pearson's correlation coefficient (r) and p value (P) are indicated on each panel. Main leaf traits comprising PC1 are indicated by grey arrows on the x axis and grouped into leaf resource acquisition (higher SLA, LNC and LPC) and conservation traits (higher LDMC and leaf C : N).

4 Discussion

We experimentally tested the effects of increased thaw depth and nutrient availability on plant traits of four tundra shrub species. While no strong responses to permafrost thaw were observed, our findings did show a coordinated response of leaf traits to fertilization, i.e. from a strategy of conservation of resources towards more rapid resource acquisition at leaf level, as we had hypothesized. Stem traits also tended towards a coordinated response to fertilization, though to a lesser extent. Moreover, one of the two deciduous (i.e. Betula nana) and one of the two evergreen species (i.e. Vaccinium vitis-idaea) showed a coordinated response of leaf and stem traits to fertilization along the same resource economics axis.

4.1 Treatment effects on plant traits

We expected that permafrost thaw and fertilization treatments would affect plant traits. However, our results showed that most of the plant traits responded only to the shallow nutrient addition. Plant growth in high-latitude ecosystems is highly nutrient-limited (Billings and Mooney, 1968; Shaver and Chapin, 1980; Epstein et al., 2000). Nutrient addition releases shrubs from this limitation and promotes their growth and biomass production (Chapin and Shaver, 1996; DeMarco et al., 2014; Iturrate-Garcia et al., 2017). Nutrient addition released shrubs from growth limitation as evidenced by the plant trait changes we found, such as greater height, SLA and leaf nutrient concentration (Hudson et al., 2011; Reich, 2014). It is notable that in this short time frame, reduced soil temperatures and permafrost thaw depth were measured under the strong fertilization treatment. Despite less favourable soil physical conditions, shrubs followed a more acquisitive growth strategy under fertilization as compared to the permafrost thaw treatment, which had higher soil temperatures and deeper thaw depth.

The fact that plant traits were less responsive to permafrost thaw than to fertilization might be explained by the relatively large amount of nutrients added to the fertilized plots. Parallel soil analyses were performed in the experimental plots by Wang et al. (2017; supplementary material). These analyses, based on buried resin bags (Giblin et al., 1991), showed no increase of exchangeable nutrients in the unfertilized unheated and heated plots for nitrogen (N) and phosphorus (P) at any of the analysed depths. In the fertilized plots, they found that the nutrient amount in the top soil layer (< 5 cm) increased by 4 (N) and 5 (P) times. A significant but smaller increase was found only for N at deeper layers (25 cm). These findings, therefore, show that the nutrient amount added to the plots with the fertilization treatment was greater than the nutrient amount released by the thawing treatment.

In plots with heated cables, enhanced nutrient availability is expected through acceleration of soil organic matter mineralization (Knorr et al., 2005; Schaeffer et al., 2013), besides permafrost thaw and related release of nutrients. Hartley et al. (1999) found effects of soil warming on subarctic shrub growth by using heating cables buried at 5 cm depth, which increased the soil temperature by 5 ^∘C. In our study, however, the heating cables were buried at 15 cm below the surface in order to increase the permafrost thaw without increasing the soil temperature of the shallower layers (< 15 cm). Consequently, most warming was in the mineral soil layers below 15 cm, whereas the increase of soil temperature in the shrub root layer was lower than the threshold (1 ^∘C or greater) needed for increasing nutrient mineralization (Schmidt et al., 1999). In addition, soil moisture conditions along the vertical profile affect the energy partitioning within the soil. Humid conditions close to the permafrost table, where the cables were buried, can promote energy partitioning towards permafrost thawing instead of towards soil heating. Under these conditions, increasing the input energy of the soil might result in a low temperature rise (i.e. high soil thermal conductivity), limiting in turn greater mineralization rates in plots with heated cables. The depth of soil layers at which nutrients were available for plants and that of shrub rooting might also explain the different trait responses to the treatments. Most of the root biomass of the shrub species studied occurs at shallow soil depth (ca. 5–10 cm), which is shallower than the permafrost thaw depth during the growing season (Churchland et al., 2010; Keuper et al., 2012; Wang et al., 2017).

4.2 Coordinated leaf trait response to nutrient addition

Resource availability is thought to be one of the main drivers of plant strategy selection (Grime, 2006; Ordoñez et al., 2010). In Arctic tundra, where resource availability is low, shrub species adopt a conservative strategy with slow growth and tissue turnover, which enhances plant survival under harsh conditions (Chapin et al., 1993). However, the “slow traits” associated with the conservative strategy are disadvantageous in the case of higher resource availability as shrub species could be outcompeted (e.g. through shading) by other species with faster growth and biomass production (Reich, 2014). In Arctic tundra, graminoid species – particularly grass species – are expected to shade and outcompete shrubs, as suggested by warming and fertilization experiments carried out at tundra sites (Dormann and Woodin, 2002: Gough and Hobbie, 2003; Wang et al., 2017). Our results show that species with similar resource economic strategies cluster into groups – deciduous and evergreen plant functional types – defined by their covarying leaf traits (Reich et al., 1997, 1999). On unfertilized plots, the deciduous shrub species B. nana and S. pulchra were characterized by leaf traits associated with faster resource acquisition: high SLA and leaf nutrient concentration and low LDMC and C : N. In contrast, the evergreen shrub species L. palustre and V. vitis-idaea were characterized by leaf traits associated with resource conservation, as expected due to a slower tissue turnover as compared to deciduous shrubs (Chapin and Shaver, 1996).

We found different plant trait responses to fertilization with PFT for most leaf traits. Despite these differences, the increase of nutrients promoted a common coordinated response of leaf traits of all species, which reflects a change in resource economics from conservation to faster acquisition, even in the case of the evergreen species. Thus, there appears to be a comparable shift towards resource acquisitiveness in the leaf economics spectrum both between PFTs, i.e. from evergreen to deciduous (Wright et al., 2004; Freschet et al., 2010; Díaz et al., 2016), and within species (this study; Aerts et al., 2012). Since deciduous shrubs have been found to expand much more than evergreen shrubs in biomass and abundance in response to fertilization, both in Eurasian and North American tundra (van Wijk et al., 2003), our findings point to a possibly important positive feedback between species turnover and intraspecific change with respect to resource economics traits.

4.3 Stem traits response to nutrient addition

Stem traits were less responsive to treatments than leaf traits, which might be explained by the relatively short time frame of the experiment. Turnover of wood tissue is slower than that of leaf tissue (Negrón-Juárez et al., 2015). Thus, stem traits might require more time to show responses. Furthermore, the age heterogeneity of the selected shrubs might mask stem trait responses. Older individuals have higher SSD and greater bark thickness than younger ones (Woodcock and Shier, 2002; Patiño et al., 2009; Poorter et al., 2014). Therefore, stem trait responses might become statistically significant when longer-term experiments are run and shrubs within the same age class (i.e. similar stem diameter) are selected.

Under nutrient addition, we found that coordinated stem trait response tended towards greater height and SWC and lower SSD. These findings are in line with previous studies showing a negative relationship between wood density and water content (Dias and Marenco, 2014). Stems with lower SSD have less space filled with cell walls than those with higher SSD, and therefore more water can be stored within the stem wood (McCulloh et al., 2011; Dias and Marenco, 2014). Woody species with denser wood grow slower, have less wood water content and produce smaller and thicker leaves, which might be associated with a whole-plant strategy (Bucci et al., 2004; Wright et al., 2004; Ishida et al., 2008; Chave et al., 2009). However, our results showed that coordination between stem height PC1 and leaf PC1 was only significant for half of the species. For Betula nana and Vaccinium vitis-idaea, the significant relationship between both axes suggests that these species coordinated stem and leaf traits (e.g. conservative trade-off at stem and leaf levels), resulting in a whole-plant strategy. The lack of coordination between stem and leaf traits for Salix pulchra and Ledum palustre suggests that, for certain species, functional trade-offs at stem and leaf levels may operate partly independently (Fortunel, Fine and Baraloto, 2012).

4.4 To grow or to defend

Our findings suggest that shrubs will grow taller, acquire more resources and allocate them to produce larger leaves at lower cost (thinner leaves with lower LDMC and C : N). These changes in plant traits, together with an expected faster growth, will come at a cost for shrubs: a decrease of their stress resistance (growth–defence trade-off) (Chapin et al., 1993; Chave et al., 2009; Iturrate-Garcia et al., 2017). The faster resource acquisition will make shrubs more vulnerable to herbivory due to higher leaf nitrogen content (Mattson, 1980; Díaz et al., 2016) and to adverse environmental conditions (i.e. low nutrient availability) as a consequence of low nutrient tissue reserves (Reich, 2014). We also found that shrubs with more rapid resource acquisition grew taller but without increasing their bark thickness and SSD, which might enhance shrub vulnerability to pests, mechanical and hydraulic failure, and extreme climatic events (Baraloto et al., 2010; Reich, 2014; Díaz et al., 2016).

4.5 Shrub–climate feedbacks

Vegetation is strongly coupled with environmental conditions (Wookey et al., 2009; Medinski et al., 2010). Shrubs will be affected by climate warming, with resultant changes in plant strategy and traits, affecting species diversity and ecosystem functions, such as carbon cycling and the surface radiation budget (Chapin et al., 1996; Beringer et al., 2005; Myers-Smith et al., 2011). The carbon uptake associated with increasing shrub growth and biomass production together with longer turnover time due to carbon storage in branches as compared to leaf material will affect the carbon cycle (Hobbie et al., 2002; Mack et al., 2004). Moreover, shrub trait changes may affect ecosystem processes as well. The production of low-cost tissues might accelerate litter decomposition because these tissues are easier to decompose than expensive ones (McLaren et al., 2017).

Our results suggest that tundra shrubs will be affected by increased nutrient availability in shallow soil layers. Deeper-rooting species, such as graminoids, may benefit more from nutrient release by permafrost thaw in deep soil layers (Keuper et al., 2017; Wang et al., 2017). In competition with graminoids, shrubs will pre-empt nutrient and light resources by growing faster and taller, producing denser canopies and leaves with greater photosynthetic area (Chapin and Shaver, 1996; Hudson et al., 2011; Elmendorf et al., 2012; Díaz et al., 2016). Bryophyte and lichen diversity is expected to decline due to the increase of shading and litter deposition associated with those changes (Cornelissen et al., 2001; van Wijk et al., 2003; Elmendorf et al., 2012; Lang et al., 2012). As a consequence of the cryptogam decrease, the thermal insulation of the permafrost might be reduced (Blok et al., 2011a), promoting permafrost thaw and the release of carbon (e.g. in the form of methane) to the atmosphere (Schuur et al., 2008; Schaefer et al., 2011). However, shrub cover increase has been reported to reduce summer permafrost thaw locally (Blok et al., 2011b; Nauta et al., 2015; Wang et al., 2017). While these studies discussed shading effects as main cause, our detailed trait analysis suggests additional mechanisms associated with water demand. Shrubs under nutrient addition showed greater SLA, lower LDMC and higher water content of leaves and stems, pointing towards enhanced water demand through higher photosynthetic potential and evapotranspiration. The higher water demand might deplete soil water resources, as suggested by the lower soil moisture and summer soil temperature in the fertilized plots (for detailed results see supplementary material in Wang et al., 2017), where deciduous shrubs increased most. This depletion might result in reduced permafrost thaw through decreasing soil moisture, thermal conductivity, heat flux and temperature, which suggest that shrub shading might not be the only driver of the reduced permafrost thaw. Water demand by plants, especially shrubs, might be at least as important, as also documented in Juszak et al. (2016). Interestingly, soil moisture has been found as a potential growth co-limiting factor of tundra shrubs (Blok et al., 2010, Myers-Smith et al., 2015). However, shrubs might be released from water limitation by the predicted concomitant increase in precipitation. Related effects on shrub growth, community composition and feedbacks with the permafrost system and the atmosphere remain to be tested.

5 Conclusions

The climatic conditions projected for the Arctic, the shrub growth sensitivity to climate and the importance of shrub–climate feedbacks for ecosystem functioning suggest that a special effort should be made to better understand future tundra changes and adaptation to the new climatic conditions. Here, we presented the response of a wide set of traits of selected dominant species in tussock tundra to permafrost thaw and increased nutrient availability. This response can be considered a step towards more realistic dynamic global vegetation models, although generalization should be considered cautiously due to the short time frame of the response, the spatial heterogeneity of Arctic regions and the complexity of shrub–climate feedbacks. According to our results, coordinated trait responses representing the whole plant (including wood and bark traits, as in our study, and ideally also root traits) instead of single trait responses are needed for a more robust prediction of shifts in vegetation, ecosystem processes and related climate–vegetation feedbacks.