Impact of Multiple Ecological Stressors on a Sub-Arctic Ecosystem: No Interaction Between Extreme Winter Warming Events, Nitrogen Addition and Grazing

Climate change is one of many ongoing human-induced environmental changes, but few studies consider interactive effects between multiple anthropogenic disturbances. In coastal sub-arctic heathland, we quantified the impact of a factorial design simulating extreme winter warming (WW) events (7 days at 6–7°C) combined with episodic summer nitrogen (+N) depositions (5 kg N ha-1) on plant winter physiology, plant community composition and ecosystem CO2 fluxes of an Empetrum nigrum dominated heathland during 3 consecutive years in northern Norway. We expected that the +N would exacerbate any stress effects caused by the WW treatment. During WW events, ecosystem respiration doubled, leaf respiration declined (-58%), efficiency of Photosystem II (Fv/Fm) increased (between 26 and 88%), while cell membrane fatty acids showed strong compositional changes as a result of the warming and freezing. In particular, longer fatty acid chains increased as a result of WW events, and eicosadienoic acid (C20:2) was lower when plants were exposed to the combination of WW and +N. A larval outbreak of geometrid moths (Epirrita autumnata and Operophtera brumata) following the first WW led to a near-complete leaf defoliation of the dominant dwarf shrubs E. nigrum (-87%) and Vaccinium myrtillus (-81%) across all experimental plots. Leaf emergence timing, plant biomass or composition, NDVI and growing season ecosystem CO2 fluxes were unresponsive to the WW and +N treatments. The limited plant community response reflected the relative mild winter freezing temperatures (-6.6°C to -11.8°C) recorded after the WW events, and that the grazing pressure probably overshadowed any potential treatment effects. The grazing pressure and WW both induce damage to the evergreen shrubs and their combination should therefore be even stronger. In addition, +N could have exacerbated the impact of both extreme events, but the ecosystem responses did not support this. Therefore, our results indicate that these sub-arctic Empetrum-dominated ecosystems are highly resilient and that their responses may be limited to the event with the strongest impact.

+N would exacerbate any stress effects caused by the WW treatment. During WW events, ecosystem respiration doubled, leaf respiration declined (−58%), efficiency of Photosystem II (Fv/Fm) increased (between 26 and 88%), while cell membrane fatty acids showed strong compositional changes as a result of the warming and freezing. In particular, longer fatty acid chains increased as a result of WW events, and eicosadienoic acid (C20:2) was lower when plants were exposed to the combination of WW and +N. A larval outbreak of geometrid moths (Epirrita autumnata and Operophtera brumata) following the first WW led to a near-complete leaf defoliation of the dominant dwarf shrubs E. nigrum (−87%) and Vaccinium myrtillus (−81%) across all experimental plots. Leaf emergence timing, plant biomass or composition, NDVI and growing season ecosystem CO 2 fluxes were unresponsive to the WW and +N treatments. The limited plant community response reflected the relative mild winter freezing temperatures (−6.6 • C to −11.8 • C) recorded after the WW events, and that the grazing pressure probably overshadowed any potential treatment effects. The grazing pressure and WW both induce damage to the evergreen shrubs and their combination should therefore

INTRODUCTION
The Arctic is experiencing extreme weather events more frequently due to climate change, causing high mortality rates among species when events surpass survival thresholds (Liston and Hiemstra, 2011;Kivinen et al., 2017). Especially during the winter period more winter extreme events are expected, such as rain on snow, unseasonal warm periods, ground ice formation and loss of snow cover (Vikhamar-Schuler et al., 2016). Abrupt changes in winter snow cover and depth following mid-winter thaw events (e.g., due to temperature rise from −20 • C to +5 • C in just 24 h) negatively affect plant survival, detectable across landscapes through remote sensing, because of reduced snow insulation against low temperature extremes (Bokhorst et al., 2009. In addition, loss of snow cover during a midwinter melt can induce plant physiological activity (i.e., respiration and fluorescence), disrupting winter dormancy of plants, and often resulting in high mortality due to drought and frost stress (Ögren, 1996;Schaberg et al., 1996;Bokhorst et al., 2009Bokhorst et al., , 2010. These extreme winter warming events (WW) occur against a background of gradually increasing temperature and extreme events such as intense herbivory by seasonal population explosions of defoliating insects (Jepsen et al., 2008) and nitrogen (N) deposition events originating from the industrialized regions further south (Karlsson et al., 2013). A combination of stressors may enhance the individual effects (Crain et al., 2008;Phoenix et al., 2012;Liess et al., 2016) and therefore have large impacts on plant development, such as phenology (Bokhorst et al., 2008), but also the plant community composition and its role in the sub-arctic terrestrial carbon budget. As such, a combination of these stressors may, through changing the dominance of plant functional types, affect the carbon balance of these subarctic ecosystems with potential feedbacks to greenhouse gas induced climate warming (De Deyn et al., 2008). Moreover, these community responses may differ greatly from summer warming responses and we therefore, need to understand the species and community responses for future sub-arctic vegetation predictions.
Recent experimental field studies of WW events in the subarctic have shown great vulnerability of evergreen dwarf shrubs to such events while deciduous plants and cryptogams are less affected (Bjerke et al., , 2017aBokhorst et al., 2015), which may govern future sub-arctic vegetation changes (Bokhorst et al., 2015). Differences in vulnerability between plant types are partly determined by differences in exposure of overwintering tissue to freezing, and in part by physiological adaptations, such as winter dormancy and changes in the membrane fatty acid composition . Higher N availability can also affect the plants' cell and physiological characteristics associated with drought and frost susceptibility (Carroll et al., 1999;Schaberg et al., 2002). Nitrogen input from agricultural practices, fossil fuel burning, and biomass burning can reach high latitudes (Forsius et al., 2010;Karlsson et al., 2013) and make sub-arctic plants more vulnerable to the temperature variability of WW events (Power et al., 1998;Phoenix et al., 2012). However, it is unclear which plant types will be most vulnerable to the combination of high N availability and WW events and whether this will result in a vegetation regime shift.
To address the impact of the combined stressors of N input and WW events on the community composition of sub-arctic vegetation we initiated an experiment that combined these factors in sub-arctic Norway. The simulation of WW events was done during three consecutive winters in February (2014)(2015)(2016), while N additions were applied during the summer months. However, during the growing season directly following the first WW simulation (2014), the field sites were subject to intense grazing by caterpillars of the geometrid moths Epirrita autumnata and Operophtera brumata (Bjerke et al., 2017bPepi et al., 2017). Such grazing pressure greatly affect the growth response of dwarf shrubs, as these are targeted by the caterpillars that drop onto the ground cover once the birch trees have been defoliated (Lehtonen and Heikkinen, 1995;Malmström and Raffa, 2000;Jepsen et al., 2013;Karlsen et al., 2013). Grazing pressure can result in sub-arctic vegetation regime shifts where dominant dwarf shrubs are suppressed and herbs and grasses emerge (Tømmervik et al., 2004;Olofsson et al., 2013). This therefore resulted in a factorial field experiment where the impacts on sub-arctic heath vegetation of +N and WW simulations were compared but also included intense grazing by caterpillars across all treatments. We hypothesized that (1) summer N additions will negatively impact the winter plant physiological adaptations to frost and drought. Because N additions are known to exacerbate other stressors (Power et al., 1998;Phoenix et al., 2012) we expect N additions to result in an increase of the damage caused by WW events and grazing impact. The grazing and WW effects will be strongest on evergreen dwarf shrubs because this plant type appears most vulnerable to WW events (Bokhorst et al., 2015 while also being heavily grazed when caterpillars fall to the ground Karlsen et al., 2013). Therefore, we hypothesize that in response to these extreme events (2) the plant community may start to shift from an evergreen dwarf shrub dominated community to one with higher dominance of cryptogams, grasses, herbs and deciduous plants with potential greater turnover of carbon flux rates and limiting carbon sequestration of these sub-arctic ecosystems.

MATERIALS AND METHODS
The study site was located on the small island Håkøya, situated in the fiord Balsfjorden between the larger island Kvaløya and the mainland (Tromsø, Norway, 69.66 • N 18.78 • E, 30 m a.s.l.). The western part of Håkøya, where this study was located, has a mosaic of open deciduous woodland dominated by Betula pubescens Ehrh. and treeless heathland dominated by the dwarf shrub E. nigrum L. The climate is relatively mild for these latitudes due to the warm Norwegian current (which is a branch of the North Atlantic current), resulting in mean summer and winter temperatures of 12 • C and −4 • C, respectively (Førland et al., 2009). Annual precipitation is ca. 1000 mm and the winter snowpack typically reach 60-80 cm depth. The experimental site was situated in an area with sparse distribution of birch trees (B. pubescens) and dominated by E. nigrum with a dense cover of the moss Pleurozium schreberi (Willd. ex Brid.) Mitt. Subdominant plant species included Vaccinium vitis-idaea L., V. uliginosum L., V. myrtillus L, Cornus suecica L., Avenella flexuosa (L.) Drejer, the moss Polytrichum commune Hedw. and Cladonia lichens.
The experiment consisted of 24 plots (1 m × 2 m) with four treatments replicated six times: control (C), N addition (+N), extreme winter warming events (WW), and WW with N addition (WW+N). Ammonium nitrate solutions (5 kg N ha −1 ) were applied by watering can (2 L volume) across each +N treatment plot three times during the growing season at monthly intervals. The N additions were at the lower limit of effect doses for most plants (Phoenix et al., 2012) and chosen to reflect realistic scenarios of N input resulting from airborne transport for these sub-arctic regions (Karlsson et al., 2013). The WW simulations were implemented by using four infrared heaters (800 W emitting at 3 µm; HS 2408, Kalglo Electronics Co., Bethlehem, PA, United States) that were suspended in parallel (65 cm apart) from wooden frames. This produced a thermal radiation flux of 270 W m −2 to the plots (at zero wind speed). Lamps were on between 9 and 16 February 2014, 13-20 February 2015, and 11-18 February 2016. The snow pack (60-80 cm deep) gradually melted out during 3 days after which the lamps remained on for another 4 days. Lamps were adjusted in height above the surface to achieve leaf temperatures of ca. 5 • C. Leaf temperatures were monitored at least twice per day by obtaining a reading of a thermocouple attached to the underside of a plant leaf in each plot. Lamps were turned off after 7 days and removed from the frames. The experimental plots were left exposed to ambient conditions and build-up of a new snowpack for the remainder of the winter. We marked out an additional 12 quadrats (1 m × 1 m) which were treated as control (n = 6) and nitrogen additions (n = 6) from which we could sample winter plant tissues and measure ecosystem gas fluxes without disturbing the main C and +N experimental plots. Temperature at canopy height was monitored throughout the year by temperature loggers (Hobo UA-001-08, Onset Computer Corporation, MA, United States) recording at hourly intervals in four control plots and four plots exposed to extreme WW events. Loggers were shielded from direct sunlight by a white dome cover.

Leaf Phenology and Vegetation Composition
Vegetative bud development was monitored during the subsequent growing season (early June onward). For this, 10 randomly selected shoots of the dwarf shrubs E. nigrum, V. vitis-idaea, and V. myrtillus were tagged in each plot and surveyed every week or more frequently depending on the speed of development. Due to the presence of large caterpillar numbers of geometrid moths, a large proportion (see results) of tagged shoots were grazed and new unaffected shoots were selected in spring 2014. Vegetative bud development was recorded by noting when the bud had burst and the first leaf had fully expanded (Bokhorst et al., 2008). During the 2015 WW event we quantified bud elongation of V. myrtillus, as it represents a first step in bud development (Bokhorst et al., 2010). The abundance and cover of plant species was quantified using the point intercept method in a quadrant (30 cm × 30 cm) in the middle of each plot during peak plant biomass (August). A total of 121 point counts at 2.5 cm intervals were made of the vegetation in each square by counting the number of times a vertical pin touched plant parts (Jonasson, 1988). Cryptogams were counted as present or absent, while vascular plants could be hit more than once by each vertical pin. For E. nigrum only shoots were counted rather than every leaf hit to avoid overrepresentation due to the high number of tightly packed, small needle-like leaves. Point intercept counts of vascular plants were converted to biomass using regression formulas (Supplementary Table S1) (Jonasson, 1988;Bokhorst et al., 2011). Species cover was quantified from point count survey based on presence or absence at each point. Changes in plant biomass and cover were compared from 1 year to the next starting at August 2013.

Normalized Difference Vegetation Index (NDVI)
As a measure of vegetation activity ("health") across the experimental plots we quantified NDVI by using two different handheld proximate sensor. We used a Maxmax-modified Canon camera (LDP LLC, Carlstadt, NJ, United States) where an infrared sensor replaced the normal sensor (the blue channel records the visible light and the red channel the near infrared). In addition, we used a GreenSeeker (Trimble Navigation Ltd., Sunnyvale, CA, United States) which generates its own radiation for NDVI measurements, while the MaxMax camera is a passive instrument, using reflected and incident radiation (Sakamoto et al., 2012). Plot pictures were taken during peak vegetation biomass, the second week of August, each year (2013August, each year ( -2016. During 2015 and 2016, the two types of NDVI measurements were also done at weekly intervals to monitor changes in plotlevel greenness during the growing season (May-August).

CO 2 Fluxes
Respiration rates of V. vitis-idaea (individual leaves) and shoots (2 cm) of E. nigrum (2014 and 2015) and the moss P. schreberi (2015 only) were quantified on an Infrared Gas Analyser (IRGA) (GFS-3000, Walz GmbH, Effeltrich, Germany). Measurements were done in closed cuvettes in complete darkness, at 7000 ppm H 2 O, with 380-ppm base level of CO 2 and temperature of the measuring head kept at 5 • C. We used the mean respiration rates of nine measurements taken at 15 s intervals for each sample. Samples were collected during the last day of warming, after the maximum exposure period, and 3 days after the lamps had been turned off.
Ecosystem CO 2 fluxes were measured during the WW events and in the following growing seasons (2014-2016) from May till August. Measurements were made by placing a clear chamber (20 cm × 20 cm × 20 cm) made from poly-methyl methacrylate over the vegetation and monitoring the rate of change in headspace CO 2 concentration, across nine measurements at 10 s intervals, using an IRGA (EGM-4 PP Systems, Amesbury, MA, United States). To minimize internal chamber air exchange with the external environment, plastic skirts (20 cm wide) were attached to a square frame -onto which the chamber could be attached -and weighed down with chains (and snow in winter). The square frame was slotted onto four metal pins that were fixed in the plots to ensure that CO 2 fluxes were measured at the same spot in each plot. Snow was removed 2 h before CO 2 measurements from the additional C and +N plots, to allow any build-up of CO 2 underneath the snow layer to diffuse out (Grogan and Jonasson, 2006).

Plant Physiological Measurements
As a measure of winter plant physiological activity we quantified potential activity of PSII using a mini-PAM (Walz, Effeltrich, Germany) for E. nigrum and V. vitis-idaea using leaf clips, in the experimental warmed plots and from underneath the snow as controls during early morning (06:00-06:30) when the sun was below the horizon and PAR levels were zero. Measurements were done during the last 2 days of warming when plants had the longest exposure to the warming, to quantify potential release from winter dormancy, and in addition 4 days after the lamps had been turned off to quantify the response of PSII following freezing. One leaf (V. vitis-idaea) or shoot (E. nigrum) was measured in each experimental plot during each measuring day.
To quantify changes in the fatty acid composition of the cell membranes as a result of the warming and freezing, we collected leaf samples (n = 5), while plants had been exposed to WW for 4 days and 4 days after the WW treatment had been turned off. Samples from control plots were retrieved during either of these sampling events. All samples were brought back to the laboratory (within 1-2 h) and frozen at -20 • C, freezedried, ground, and analyzed for fatty acids following the direct methylation procedure (Browse et al., 1986). Samples of 5-20 mg were dissolved in 1 mL methanolic hydrochloric acid (HCl) (1M) and an internal standard (heptadecanoic acid, C17:0) was added to a glass tube. The solution was heated to 80-100 • C for 1 h, and after cooling, 0.4 mL hexane and 1 mL of 0.9% sodium chloride (NaCl) were added to each sample. The fatty acid methyl esters (FAMES) were extracted into the hexane phase by vigorous shaking. The tubes were centrifuged for 10 min to separate the phases completely, and a sample was then taken directly from the hexane phase. Samples were stored at -20 • C until gas chromatography (GC) analyses, according to Maehre et al. (2013). The instrument used was an Agilent 6890N equipped with a flame ionization detector (FID) (Agilent Technologies Inc., Santa Clara, CA, United States) and a CP7419 capillary column (50 m × 250 µm × 0.25 µm nominal, Varian Inc., Middelburg, Netherlands). The fatty acids were identified by comparing against the commercial fatty acid standards PUFA 1, 2, and 3 (Sigma-Aldrich Chemicals Co., St. Louis, MO, United States) and the GLC standards 80, 411, and 412 (NuChec Prep. Inc., Elysian, MN, United States). The amount of each fatty acid was calculated by comparing peak area with the known amount of an internal standard (C17:0). The proportions of the single fatty acids were used in further analysis. We calculated the unsaturation to saturation ratio (U/S ratio) as the ratio between the total proportion of all unsaturated fatty acids and the total proportion for all saturated fatty acids (van Dooremalen et al., 2011). The detected fatty acid composition differed greatly between plant species, which limited direct comparisons of specific fatty acids between study species. However, changes in fatty acid composition could be quantified in response to the treatments.

Calculations and Statistical Analyses
We used repeated measures ANOVA to detect differences in leaf phenology, NDVI, and CO 2 fluxes between treatments (WW and +N) during the growing seasons of 2014-2016. To identify changes in plant biomass and cryptogam cover across the measurement years (2013-2016) we used a linear mixed effects model with treatments (WW vs. C and -N vs. +N) and years as fixed factors and block as a random factor. P-values were obtained by likelihood ratio tests of the full model with the effect in question against the model without the effect in question. A Principal Component Analyses (PCA) was used to summarize changes in vegetation composition of vascular plant and cryptogam cover between treatments, and the first two PCs were analyzed in the same way as the plant biomass. Treatment effects on fatty acid concentrations and U/S ratio were tested using a factorial ANOVA with Treatment and +N as fixed factors. We used Tukey HSD tests at P = 0.05 to identify differences in means between WW and C (with and without +N) whenever the interaction effect was significant. A visual inspection did not show any patterns in the residuals. All statistical analyses were carried out using R 3.3.0 (R Core Team, 2015).

Temperature Effects of Winter Warming Treatment
Canopy temperature increased to 6.3-7.2 • C in the WW plots while canopy temperatures underneath the snow were around freezing (−0.1-0.1 • C) in the control plots (Figure 1). Minimum temperatures were somewhat lower in WW compared to C during 2014 (−11.8 • C vs. −9.7 • C) and 2016 (−6.6 • C vs. −1.0 • C) but did not differ during 2015 (−9.1 • C). The number of freeze thaw cycles increased following WW by 67% and 57% during 2014 and 2015, respectively, and from 1 in C to 29 in WW FIGURE 1 | (A,B) Show canopy temperatures measured during 2014 and 2015 respectively in control plots (C) and those exposed to extreme winter warming (WW). (C,D) Show the monthly number of freeze-thaw cycles during 2014 and 2015 respectively. Mean daily canopy temperatures and monthly number of freeze thaw cycles during the winters of 2014 and 2015. Data points and bars are means of n = 4 with SE as error bars. * indicate significant (P < 0.05) differences between treatments. Canopy temperatures were measured following the extreme winter warming event of 2016 but not during the winter period before therefore, these are not shown.
during 2016 (Figures 1C,D). Leaf and moss temperatures were on average 5 and 10 • C, respectively, during the WW treatments.

Winter Physiology and Ecosystem Responses
Efficiency of Photosystem II (Fv/Fm) of E. nigrum increased (P < 0.001) by 77 and 80% in the WW (0.59) and WW+N (0.59) treatments, respectively, as compared to C (0.33) during 2014, indicating a release from winter dormancy. Fv/Fm of V. vitis-idaea increased (P < 0.01) by 26 and 33% in the WW (0.61) and WW+N (0.64) treatment, respectively, compared to C (0.48) during 2014. Fv/Fm no longer differed between the treatments after the warming lamps was turned off. Fv/Fm was unaffected by WW during 2015, and values across treatments were 0.68, 0.65, and 0.69 for E. nigrum, V. vitis-idaea, and P. schreberi, respectively. Winter leaf respiration decreased by 58% for E. nigrum (F 3,19 = 3.8, P = 0.027) and with 43% for P. schreberi (F 3,19 = 6.4, P = 0.005) in WW compared to C during 2015 but did not differ after lamps were turned off, and there were no differences between WW and WW+N. Vaccinium vitis-idaea leaf respiration was unaffected by the treatments (F 3,19 = 2.0, P = 0.155).
Bud length of V. myrtillus increased by 12% during the 2015 warming event, representing an increase of 0.14 mm and 0.17 mm in WW and WW+N, respectively.

Membrane Fatty Acids
The ratio of unsaturated to saturated fatty acids was not affected by the treatments for any of the study species during 2014 or 2015 (data not shown). The concentration of various membrane fatty acids declined for E. nigrum, V. vitis-idaea, and P. schreberi in response to WW during 2014 (Supplementary Table S2 and Figure 2). However, during the second WW event (2015), the response varied between treatments and species and not all fatty acids were found during both years. Five affected fatty acids were consistently lower (from 3 to 100%) during the WW of 2014 in E. nigrum compared to control plots and following the WW (Figure 2A). However, C13:0 was increased (263%) during WW compared to control plots. The fatty acids C16:1n-7, C18:1n-7, and C18:3n-3 were all lower following the 2015 WW compared to control plots, while C8:0, C10:0, and C12:0 were not found in detectable concentrations. C13:1 had highest concentrations (45% higher) following WW compared to control samples (Figures 2B,C). Nitrogen increased C18:1 n-12 by more than three times during 2014 while C12:0 was reduced by 69% during 2015. C18:0 and C20:2 (eicosadienoic acid) responded to WW in combination with +N; for C18:0 there were no Tukey post hoc differences, while C20:2 concentrations were consistently higher during WW without N compared to the other treatments.  Affected fatty acids (n = 5) were consistently lower (from 67 to 100%) during the WW of 2014 in V. vitisidaea compared to frozen samples from control plots and following the WW (Figure 2D). During 2015, C18:2 n-6 declined by 38% following the WW event compared to control samples, while C21:0 and C20:2 increased by 219 and 118%, respectively, following the WW event (Figures 2E,F). The remaining near-significant differences (P > 0.05) in fatty acid concentrations were found during and after the WW event. Nitrogen did not affect fatty acid concentrations during 2014. During 2015, C21:0 and C20:2 were both reduced (58 and 54%, respectively) when +N was applied. C14:0 had a significant WW × N interaction, but Tukey post hoc testing did not indicate significant treatment differences (Supplementary Table S2).
There were three affected fatty acids (C15:0, C18:1 n-12, and C18:3 n-6) responding to the interaction between WW and +N during 2014 in P. schreberi but all in a different way (Figure 2G). Concentrations of C15:0 were approximately 10 times higher in C +N compared to C without N but not compared to any of the other treatments. Concentrations of C18:1 n-12 were highest in C (without +N) compared to all other treatments, while concentrations of C18:3 n-6 were highest for WW+N compared to all other treatments. During the winter of 2015, the majority of the affected fatty acids (n = 16) were consistently higher (9-65%) in control plots, while only C18:1 n-12 was higher during and after the WW event (54 and 81%, respectively) (Figures 2H,I).

Ecosystem Responses in the Growing Season
Leaf emergence timing was not affected by the treatments during any of the years, although there were differences in the percentage of emerged leaves between treatments across the measuring periods (Supplementary Figure S1 and Table 1). Empetrum nigrum had higher total percentage of emerged leaves in WW (65%) compared to control (33%) at the final measuring date in 2014. There were no treatment effects on V. myrtillus leaf emergence during any of the years. While V. vitis-idaea had highest percentage (40%) of total fully developed leaves in the +N treatment at the final measuring dates in 2016.
Plant biomass and cryptogam cover were not affected by WW, +N or their interaction ( Table 2). Instead, biomass of most vascular plant species declined following the second WW events in all plots (including C) due to the caterpillar outbreak. Especially the dwarf shrubs were targeted by the caterpillars resulting in a large proportion of grazed shoots of E. nigrum (81%), V. myrtillus (87%), and V. vitis-idaea (17%) during 2014. There were no significant differences in the number of grazed shoots between treatments. Avenella flexuosa and C. suecica biomass were unaffected by the treatments and did not appear to be affected by caterpillar grazing. Plant biomass increased again during the last (2016) growing season (Figure 3). Moss cover (P. schreberi and P. commune) reached its peak following the second WW event (summer 2015) but declined again in 2016 (Figure 4). The lichen Cladonia uncialis and the liverwort Ptilidium ciliare had a stable cover between years and were unaffected by the treatments (Table 2 and Figure 4). The PCA analysis indicated that the vegetation community did not change between treatments or years ( Table 2).
Plot NDVI was not consistently affected by WW, +N or their interaction. During 2016 WW plots were 5% lower NDVI values compared to C in spring while later on in the season there was no difference (Table 1 and Supplementary Figure S2). NDVI plot values were on average 38% lower during August 2015 compared to starting values ( Table 1). CO 2 fluxes increased with time during the growing seasons ( Table 1), but there were no differences in ecosystem ER, NPP, and GPP between the treatments during the growing seasons (Supplementary Figure S3), with the exception of some lower NPP values under +N treatment in 2016.

DISCUSSION
We anticipated a vegetation regime shift from evergreen dwarf shrubs (heath) to one of cryptogams, grasses, and deciduous plants (meadow) as a result of WW with this being exacerbated by the N treatment. Although evergreen dwarf shrubs were reduced in biomass following the grazing and WW events, the communities did not significantly differ in species composition between treatments and years, indicating a strong resistance to change in these E. nigrum dominated ecosystems. Furthermore, there was no strong evidence for any interaction between the multiple extreme events, indicating that multiple stresses do not necessarily lead to accumulation of biological stress responses (Thompson et al., 2008;Jackson et al., 2016).
We expected that +N would exacerbate the effects of WW and grazing pressure but no such responses were observed. We expected that N additions would promote plant growth, but simultaneously increase vulnerability to freezing stress due to changes in cell and physiological characteristics associated with drought and frost susceptibility (MacGillivray et al., 1995;Carroll et al., 1999;Schaberg et al., 2002). However, the added N-levels (5 kg ha −1 ) did not lead to measurable responses in the control plots suggesting that these coastal E. nigrum vegetation types should be able to withstand N depositions events (Forsius et al., 2010;Karlsson et al., 2013). In addition, mid-winter bud swelling, Fv/Fm and leaf respiration, indicators of a potential release from winter dormancy responded similarly to WW with or without +N. The moss P. schreberi did show changes in membrane fatty acids to +N, but not consistently in combination with WW. The fatty acid C20:2 increased in WW-N compared to WW+N for both E. nigrum and V. vitis-idaea, indicating that soil N availability may affect cell membrane adaptations to freezing stress. However, as we did not find consistent phenology and plant biomass responses to WW+N in the following growing seasons, it appears that the role of N on the plants' winter physiology appears limited in this study. Nitrogen additions have been shown to increase insect damage to Calluna vulgaris communities (Power et al., 1998) but this was not observed in our study. Overall, the combination of stress events and low dose N deposition does not appear to seriously affect these coastal sub-arctic E. nigrum communities.
In support of our second hypothesis, the evergreen dwarf shrub E. nigrum declined in response to the grazing by geometrid caterpillars while the herb C. suecica, not often grazed upon, increased over the study period. Empetrum nigrum is frequently targeted by caterpillars when their main food plant birch trees have been defoliated , which explains the decline of this dwarf shrub over time. The increase of C. suecica was also found elsewhere in northern Fennoscandia as a combined effect of grazing pressure, increased precipitation and deposition of nitrogen through precipitation and caterpillar fecal matter (Tømmervik et al., 2004;Karlsen et al., 2013). However, we did not find any such response from N additions in our study, which indicates that the applied levels were insufficient to elicit a measurable plant growth response or that caterpillar fecal N input was much higher. In addition, moss cover was highest during 2015, reflecting the opening of the vascular plant cover canopy, as also observed following similar extreme events in physiognomically comparable open birch woodland in a more continental upland region of northern Scandinavia (Bokhorst et al., 2015). However, the evergreen species, such as E. nigrum and V. vitis-idaea, quickly recovered their biomass to starting levels in this study, indicating that this coastal E. nigrum dominated vegetation were fairly resistant to the combination of WW events, N additions and invertebrate grazing. The ecosystem CO 2 fluxes measured during the growing season also did not show differences between treatment plots indicating that the plant community shifts were too limited to impact this ecosystems carbon balance (Luyssaert et al., 2007;De Deyn et al., 2008), although the 2016 values were notably higher than during the 2 previous years (Supplementary Figure S3).
The muted response from the plant community, and NDVI values to the extreme WW events, may result from the intensive grazing impact overall, including the control plots, which could have diluted the effect of the imposed winter stress. The relative importance of such events, and their interaction, in shaping subarctic plant communities remains a complex interplay depending in large part on intensity, timing and duration (Olofsson et al., 2013). Furthermore, sub-dominant species cover was not affected enough to affect plot communities dominated by E. nigrum and P. schreberi, while high shoot defoliation of dwarf shrubs recorded during 2014 could have been partly compensated later on in the season . In addition, studies on C. vulgaris-dominated ecosystems show that plant community responses to multiple environmental drivers can be affected by successional stage and disturbance regimes across Europe (Kröel-Dulay et al., 2015). As such, the high dominance of E. nigrum in our field site may reflect a successional stage that is more resistant than stages where other species comprise a larger proportion of the community or these coastal habitats have plants with greater physiological adaptability to cope with multiple stressors than habitats with more stable climates. With respect to the WW events, more pronounced plant responses have been reported when temperatures were lower (Bokhorst et al., 2009;Bjerke et al., 2017b), which is an important factor for frost droughtinduced plant mortality (Kreyling et al., 2008;Bokhorst et al., 2010Bokhorst et al., , 2018. The relative mild freezing temperatures may have allowed for winter physiological adaptations by E. nigrum and V. vitis-idaea, such as the greater proportion of longer chain fatty acids (Dalmannsdóttir et al., 2001;Schaberg et al., 2002). The N addition levels were chosen at the lowest concentration effectresponse by most plant communities (Phoenix et al., 2012) as a realistic scenario, while higher N levels for a longer duration (>3 years) may result in additive biological responses (Power et al., 1998;Phoenix et al., 2012). Similarly, grazing intensity on plants is not consistent during each herbivore outbreak (Olofsson et al., 2013) indicating that sub-arctic plant community responses to multiple extreme events may depend predominantly on the event with the strongest intensity.

CONCLUSION
The findings indicate that coastal sub-arctic dwarf shrub plant communities and growing season CO 2 fluxes appear largely unaffected to a combination of extreme WW events, grazing and nitrogen additions in summer. This response may in part be the result of inherent physiological adaptations, such as changing the membrane fatty acids (Dalmannsdóttir et al., 2001), during winter and the mild winter freezing temperatures . The grazing impact appears to have overshadowed all plant responses to the WW and +N treatments indicating that multiple extreme events, that in theory can enhance each other's effects (Power et al., 1998;Phoenix et al., 2012), do not necessarily increase biological stress responses.