Responses of low Arctic tundra plant species to experimental manipulations: Differences between abiotic and biotic factors and between short- and long-term effects

ABSTRACT Climate change in arctic tundra is projected to affect air temperature, snow depth, soil fertility, and caribou herbivory, which may alter plant community composition by shifting niche space to favor particular species’ life history strategies. We report responses of a Canadian mesic birch hummock tundra plant community to a range of manipulative experiments (greenhouse warming, fertilization, snow fence, and caribou exclosure treatments). Aboveground biomass of each plant species was measured in the same permanent 1 m2 areas using the point frame method in 2005, 2011, and 2017. Although the greenhouse treatment had few effects on individual species, total vascular plant community biomass was enhanced between 2011 and 2017. Furthermore, species’ biomass across all control plots was stable from 2005 to 2011 but increased significantly from 2011 to 2017, with air temperatures also warmer over that same period. Species responded to high-level nitrogen and high-level nitrogen and phosphorus combined additions, with deciduous shrubs and graminoids increasing and evergreen shrubs decreasing. The snow fences and caribou exclosures had little effect on species biomass. Although vegetation greening trends have been reported in arctic environments that are primarily influenced by maritime climate, our study is one of the first to provide plot-based evidence of recent plant biomass increases in the low Arctic’s continental interior.

Warming in arctic tundra will affect ecosystem processes, such as carbon (C) and nutrient cycling, productivity and decomposition, and the consequent changes in greenhouse gas fluxes and surface albedo are likely to cause further climate change (Sturm et al. 2005;Dorrepaal 2007;Wilson and Humphreys 2010;. Plant community structure (i.e., species richness, composition, relative abundances, and interactions) is important in mediating these processes and feedbacks (Yang et al. 2011;Cantarel, Bloor, and Soussana 2013). According to fundamental niche theory, species can coexist alongside each other if they differ in the extents to which their growth is limited by the various growthdetermining environmental factors (Dorrepaal 2007), resulting in a unique multidimensional niche space for each species (Hutchinson 1957). In arctic tundra, the temporal ) and spatial availability of resources is highly variable, creating multiple niches into which different tundra species have evolved and can coexist McKane et al. 2002). Determining the particular suite of environmental factors and their relative importance for growth of each species is therefore critical to understanding and predicting the responses of tundra communities to climate change.
Shrubs (deciduous and evergreen) are woody perennial species that can live for centuries and are often the dominant vascular plant species across the low Arctic tundra biome (Bliss and Matveyeva 1992;CAVM Team 2003). Shrubs possess multiple traits that make them the ideal functional group for quantifying tundra vegetation responses to long-term environmental changes (Myers-Smith et al. 2015). Their long life span and predominantly vegetative propagation mean that they experience a larger amount of environmental change per generation (Havström, Callaghan, and Jonasson 1993). In addition, the tussock-forming graminoid species Eriophorum vaginatum L. dominates vast areas of the Arctic (Kummerow et al. 1988) and persists for decades to centuries, profoundly affecting tundra vegetation community dynamics (Wein and Bliss 1974;Mark et al. 1985). Therefore, the responses of individual shrub species and E. vaginatum to long-term environmental changes, while also facing concomitant changes in relative competition from co-existing neighboring plant species, will be a strong ultimate determinant of the future community structure of arctic vegetation (Havström, Callaghan, and Jonasson 1993).
The Arctic is warming two to three times more rapidly than lower latitudes (Intergovernmental Panel on Climate Change 2013) owing to climate amplification effects (Serreze and Francis 2006;Serreze et al. 2009;Cohen et al. 2014), and therefore one of the major current research themes in arctic tundra ecology is understanding and predicting the effects of ongoing air temperature rises on vegetation dynamics. The Arctic's particularly high warming rate is expected to exert dramatic environmental selection pressures on its biota (Post et al. 2009). Tundra plant species have been suggested as among the most responsive biological groupings because they operate close to their lower physiological temperature tolerance limits and thus are expected to exhibit significant biological responses to even relatively small temperature increases (Jonasson et al. 1996;Atkin et al. 2005;Stinziano and Way 2014).
The widespread expansion of arctic shrubs over the past thirty years based on remote sensing data, repeat photography, and vegetation surveys (Tape, Sturm, and Racine 2006;Pouliot, Latifovic, and Olthof 2009;Bhatt et al. 2010; has been consistently attributed to recent arctic air temperature increases. However, neither climate warming nor shrub increases have been observed at all arctic tundra monitoring sites . In particular, the primary greening trends have been reported for western Alaska, along the northern coast of Canada, and in northeastern Canada (Bhatt et al. 2010;Ju and Masek 2016;Myers-Smith et al. 2019). By contrast, there are locations in the central Canadian and Alaskan Arctic-the interior continental regions-where decadal warming and likewise vegetation greening are either absent or patchy (Ju and Masek 2016;Bonney, Danby, and Treitz 2018). Furthermore, several long-term monitoring field sites in those regions have observed no consistent trends in either surface air temperature (Toolik field station, interior Alaska, 1989Alaska, -2014Hobbie et al. 2017) or plant community composition (Wolf Creek drainage basin, Yukon, 1999Yukon, -2008Pieper et al. 2011). Indeed, a very recent study at our Daring Lake research site (Northwest Territories) indicated significant birch shrub expansion over the period between 2006 and 2016 but little evidence of a corresponding air temperature warming trend -2006Andruko, Danby, and Grogan 2020).
Factors other than direct impacts of climate warming, such as increases in soil nutrient availability (Shaver and Chapin 1980;Zamin, Bret-Harte, andGrogan 2014), winter precipitation (Wahren, Walker, andBret-Harte 2005;Björk and Molau 2007), and changes in grazing intensity (Olofsson et al. 2009) are strong drivers of tundra plant species composition . For example, soil nutrient availability, especially nitrogen (N), is widely recognized as a key factor limiting tundra plant growth, and its impacts on plant community may be more immediate and dramatic compared to those of other environmental factors, especially when high doses of nutrients are added (Zamin, Bret-Harte, and Grogan 2014). In addition, deepened snow depth caused by increased winter precipitation is hypothesized to enhance upward growth of deciduous shrub species, which are then able to trap more snow, thus providing a positive feedback between snow, winter soil temperature, microbial activity, nutrient release, and species composition (Sturm et al. 2001;Schimel, Bilbrough, and Welker 2004). Lastly, the topdown regulation of large herbivores can also alter plant community structure through selective browsing of graminoids and deciduous shrubs compared to evergreen shrubs (Grogan and Zamin 2018). However, only a few studies have simultaneously investigated whether the impacts of these individual factors on plant species' growth are consistent in direction and magnitude among tundra species over decadal timescales (Chapin and Shaver 1996;Grellmann 2002;Wahren, Walker, and Bret-Harte 2005).
Lack of information on this question greatly limits our ability to predict future vegetation feedbacks. If there are significant interspecific variations among species' responses to these individual environmental factors and these responses are time dependent (Magnuson 1990), as opposed to a common uniform positive or negative response to all factors as assumed by large-scale vegetation models, then different patterns of change in plant community structure and "greening" effects would be expected (Klanderud and Totland 2005;Pearson et al. 2013;Myers-Smith et al. 2015). High-resolution, long-term species monitoring data are a great resource to examine these important questions. Here, we reported a time series of plant community data for a single widespread low Arctic tundra vegetation type from multiple locations across the same landscape that provide a unique opportunity, and a necessary complement to the other methods (e.g., satellite remote sensing and photography), to make improved projections of vegetation feedbacks to future climate change.
In this study we examined how a suite of experimental manipulations (summer greenhouse, two levels of N and phosphorus [P] fertilization, snow fence, and exclusion of caribou grazing) that simulate various aspects of climate change affect shoot biomass of the principal tundra shrub species and one graminoid species over the short term (six years) and long term (twelve years). For the greenhouse data set in particular, we also investigated the treatment effects on shoot biomass of the whole vascular plant community to detect any consistent responses across the species group that would not be detected at the individual species level . We also measured biomass changes of lichen and moss at the functional group level as a supplementary to our main data sets because these nonvascular plants are important C stocks in many tundra ecosystems and are very sensitive (especially lichens) to abiotic and biotic environmental changes (e.g., changes in temperature, snow depth, and herbivore browsing; Heggberget, Gaare, and Ball 2002;Joly, Jandt, and Klein 2009). Each year over the full twelve years of the study, we manipulated (1) air temperature by installing plasticcovered greenhouses each summer; (2) soil nutrient availability by low-and high-level additions of N fertilizers, a high level of P fertilizer addition, and high levels of N + P fertilizers combined; (3) winter snow depth using snow fences to enhance snow accumulation; and (4) caribou herbivory using exclosures to prevent caribou foraging. Species biomass responses were measured using the nondestructive point frame method (Jonasson 1988). Finally, we analyzed the pooled data for each species from the unmanipulated control plots for each of the treatments (n = 25 in total across the landscape) to robustly investigate whether species' shoot biomass had changed over the twelveyear study period.
Our hypotheses were that (1) shoot biomass of each species would increase in response to the greenhouse manipulation over the six-and twelve-year periods, in association with enhanced soil nutrient availability; (2) shoot biomass of each species from the twenty-five pooled control plots would increase within twelve years, and those increases would correlate with a rising ambient air temperature trend; (3) shoot biomass of each species would show significant responses to the high-level nutrient addition treatments within six years, and some would be significantly affected by the low-level N fertilization within twelve years; (4) shoot biomass of species, especially the shrub species, would increase in response to the snow fence manipulation within twelve years, due to the positive feedback between snow, soil nutrient release, and plant uptake of nutrients; and (5) shoot biomass of deciduous shrub species and the graminoid species would increase in response to the exclusion of caribou foraging within twelve years.

Study site and experimental manipulations
This study was conducted in a mesic birch hummock tundra ecosystem near the Tundra Ecosystem Research Station at Daring Lake, Northwest Territories, Canada (64°52′ N, 111°33′ W). We focused on the six most abundant vascular plant species in this ecosystem type, which belong to three growth forms : two deciduous shrubs (Betula glandulosa Michx and Vaccinium uliginosum L.), three evergreen shrubs (Rhododendron subarcticum Harmaja [formerly Ledum decumbens (Aiton) Lodd. Ex Steud], Vaccinium vitisidaea L, and Andromeda polifolia L.), and one graminoid species (E. vaginatum). In addition to vascular species, there is a well-developed moss and lichen layer (see Nobrega and Grogan [2008] and Zamin, Bret-Harte, and Grogan [2014] for more details).
The four different field experimental manipulations were established in 2004 and have been maintained ever since (see Zamin [2013] and Christiansen [2016] for full details). All manipulations were set up on plots of representative and fairly homogenous birch hummock vegetation of similar gentle slope and aspect that were assigned randomly to control and treatment. Briefly, to enhance summer air temperatures, A-frame greenhouses (1.8 m × 4.7 m each, n = 10) were covered with transparent plastic during each growing season. Triangle vents were cut out of the tops to avoid extreme maximum temperatures and to reduce humidity within the greenhouses. These greenhouses warmed the summer mean diel air and soil (at 0-10 cm depth) temperatures by 2.1-2.4°C, based on manually collected one-time data in the 2006 growing season and multiyear continuous 4-hour measures from 2008 onwards (using replicate thermocouple probes connected to CR10X dataloggers, Campbell Scientifc, Logan, Utah; Zamin, Bret-Harte, and Grogan 2014). Note that soil moisture, which is an important factor in determining species' responses to environmental change (Walker et al. 2006), was not significantly altered by the greenhouses (Zamin, Bret-Harte, and Grogan 2014). To manipulate soil nutrient availabilities, annual low-level N (LN; 1 g N m −2 year −1 ; added as ammonium nitrate from 2004 to 2015 and as urea from 2016 onwards), high-level N (HN; 10 g N m −2 year −1 ), high-level P (HP; 5 g P m −2 year −1 ; added as 45 percent phosphorus pentoxide), and high-level N + high-level P (HNHP; 10 g N m −2 and 5 g P m −2 year −1 ) were added to plots (5 m × 7 m each; n = 5) for each nutrient treatment in July or August during each growing season. To manipulate winter snow depth, snow fences (15 m long, 1.2 m high each; n = 5) were set up that reduced wind speeds on their lee sides, leading to deepened snow patches extending out approximately 20 m from both sides of each fence. These snow fences generally increased the winter peak ambient snow cover from 0.3 m to 1.0 m (Christiansen 2016), and complete snowmelt within the snow fences was delayed by one to two weeks (Nobrega and Grogan 2008;. To exclude caribou browsing, ten large patches (19.8 m × 19.8 m) that were each approximately 200 m apart were assigned alternately to control or caribou exclosure treatment (n = 5). The caribou exclosures (1.2 m height) were constructed of metal "range" fencing (aperture 15 cm × 20 cm) to exclude caribou but not entry of other herbivores (e.g., hares, meadow voles, lemmings).

Shoot biomass estimation
Aboveground biomass of the six principal vascular plant species detailed above, together with all other co-occurring plant species (mainly one forb species [Rubus chamaemorus L.] and lichens and mosses), was measured after one year (2005), seven years (2011), and thirteen years (2017) of the experimental manipulations using the nondestructive point frame method (Jonasson 1988). This method records plant species' hit data and uses regressions to infer shoot biomass and has been widely applied by researchers studying tundra vegetation dynamics (Walker 1996;Schuur et al. 2007;Zamin 2013;Zamin, Bret-Harte, and Grogan 2014;Alatalo, Jägerbrand, and Molau 2015). Measurement biases introduced using this method are considered small compared to those associated with other ecological techniques (Goodall 1952;Jonasson 1988). For example, in contrast to the harvesting method, the point frame approach can be applied to exactly the same plot at multiple sampling times, therefore eliminating the confounding effect of spatial variation, which can be very large for tundra plant communities . Furthermore, compared with other nondestructive approaches, such as using regressions with height, diameter, or percentage ground cover to estimate biomass, point framing is considered more robust and consistent and therefore has a higher predictive power (Goodall 1952).
For the point frame method, a randomly selected area (1.0 m × 1.0 m) was marked out within each replicate experimental manipulation plot in 2005. Total number of hits for each species was recorded at each of 100 evenly spaced grid point intercepts in each of these permanent areas (referred to as subplots hereafter) during the growing seasons (between early July and late August) of 2005, 2011, and 2017. Note that the subplots were measured in a random order across all manipulations to avoid any temporal bias within the summer in the effects of the particular treatments on plant community structure. To calibrate the point frame hits data to shoot biomass, all aboveground biomass (i.e., down to the transition of the green-brown moss layer) was harvested from a 40 cm × 40 cm area within each of a total of sixty-four separate areas that were adjacent to the permanent point frame subplots described above from a range of the experimental manipulation replicate plots across the years 2009 (Zamin 2013), 2011 (Zamin, Bret-Harte, and Grogan 2014), and 2018 (Appendix S1). These sixty-four harvested areas were deliberately selected to include a large range of biomass for each species because the purpose was to generate robust hits versus shoot biomass calibrations for each species. See Appendix S1 in the Supporting Information for detailed information on these harvested areas.
Both simple linear regression and power regression models of hits/biomass data have high power to predict biomass from the number of intercept pin contacts (Jonasson 1988). The power regression model is preferred because the variance in biomass may increase as the number of pin contacts increases (Jonasson 1988). However, our analysis showed that the simple linear model yielded higher explained variance than the power model in most cases (Appendix S2), and therefore we chose the linear model for all data reported in this study. See Appendix S2 in the Supporting Information for the linear regression parameters. Separate calibration equations were used for leaf and stem tissues for each of the shrub species, whereas only one calibration equation was used for all aboveground biomass for the graminoid E. vaginatum. For the lichen functional group and the moss functional group, we adopted calibration equations directly from Zamin, Bret-Harte, and Grogan (2014).
Although the experimental manipulations were established over the summer of 2004, we consider the 2005 point frame data as representing initial biomass in all plots because the hits data indicated that there were no differences in any of the species' shoot biomass between experimental treatments and their corresponding control treatments in the 2005 data set (see Appendix S3 in the Supporting Information for statistical results). This overall homogeneity in plant community structure in the first year also supports the assumption that our initial choice of location for the permanent long-term monitoring subplots for point framing was unbiased for the individual treatment versus control comparisons.
We noted substantial pooling of water at the surfaces of one of the greenhouse control plots and one of the snow fence plots due to soil subsidence in 2017, probably associated with local melting of belowground ice in this palsa landscape. This happened sometime in the previous six-year period, creating an anoxic environment that led to foliar die off and thus biomass decreases (by 37-88 percent) for the two dominant species R. subarcticum and V. vitis-idaea but a four-to sevenfold increases in E. vaginatum. However, eliminating the data from these plots did not affect our overall conclusions, so they were included in all of the reported results and further statistical analyses below.
For the greenhouse data set in particular, we investigated effects at the whole vascular plant community level to detect any consistent responses across species that would not be detected at the individual species level. In addition, we pooled the data of each species from the twenty-five control plots for the various experimental manipulations that are widely distributed across the whole study valley to robustly determine whether species' shoot biomass changed with time and whether those changes (or lack of them) was associated with changes in ambient air temperature.

Soil nutrient availability estimation
We investigated soil ammonium (NH 4 -N) and phosphate (PO 4 -P) pools and fluxes to determine any potential fertilization, greenhouse, or snow fence treatment effects on soil nutrient availability. Specifically, we estimated soil NH 4 -N and PO 4 -P pools in the fertilization and their associated control plots during the 2016 growing season using the soil chemical extraction method (Gregorich and Carter 2007). Homogenized organic soil samples (0-7 cm depth) collected from intertussocks in each replicate sampling area were extracted in distilled water and NH 4 -N and PO 4 -P concentrations in the extracts were then determined colorimetrically using the indophenol, sulfanilamide (Mulvaney et al. 1996), and molybdate-ascorbic acid methods (Kuo et al. 1996). We also estimated intraseasonal changes in soil NH 4 -N and PO 4 -P fluxes in the greenhouses, the snow fences, and their associated control plots throughout the 2017 growing season using the in situ ion exchange membrane (IEM) incubation method (Giblin et al. 1994;Qian and Schoenau 2002;Gu and Grogan 2020). Full details of our development and specific application of the IEM method including tests and sensitivity analyses are described elsewhere (Gu and Grogan 2020). Briefly, two types of membranes, cation and anion exchange membranes, were used for accumulating the soil NH 4 -N and PO 4 -P ions, respectively. These membranes were placed in soil slits at a uniform depth (i.e., extending vertically from 2 to 7 cm underneath the soil surface) in each plot for a predetermined period and were then removed and replaced with new membranes in the same slits. Six rounds of incubations in the greenhouses and their control plots and five rounds of incubations in the snow fence and their control plots were conducted with incubation periods ranging from nine to seventeen days. The retrieved membranes were eluted with 2 M NaCl in 0.1 M HCl solution, and NH 4 -N and PO 4 -P concentrations from the elution solutions were determined using the methods mentioned above.

Air temperature trends analysis
We examined the pattern of late spring-early fall mean air temperatures at our study site over the time periods matching the phases of our vegetation monitoring measurements (i.e., 2000-2005, 2006-2011, and 2012-2017), using data from the local meteorological station. This station has been operated by Shawne Kokelj and others (Water Management and Monitoring Division of the Department of Environment and Natural Resources, Government of Northwest Territories) at the Tundra Ecosystem Research Station since 1996. Note that station data from 5 August to 30 September 2012 were not available due to equipment malfunction, and this gap was replaced with data collected from another location within our study area (i.e., from the two control plots for the greenhouse treatment) that were closely correlated to the climate station data (i.e., using a linear regression equation [slope = 0.88, R 2 = 0.91, p < .01] based on the July data sets).

Statistical analyses
We examined the effects of each different experimental manipulation treatment on aboveground biomass of each vascular plant species across the time series of three sampling years, using separate two-way mixed-model factorial analyses of variance (ANOVAs), with time as the within-subjects factor and the experimental manipulation as the between-subjects factor. We used separate Student's t tests to examine the impacts of each of the summer greenhouse, snow fence, and exclosure treatments on each species' aboveground biomass (as well as on lichen and moss functional groups) in each of the sampling years (2011 and 2017) and of the summer greenhouse and snow fence treatments on IEM nutrient fluxes for each of the individual incubation periods. We investigated changes in each plant species' aboveground biomass over time across the twenty-five control plots, changes in lichen biomass over time, and the overall effects of the summer greenhouse and snow fence treatments on IEM nutrient fluxes for the whole growing season, using separate one-way repeated measures ANOVAs. We examined the effect of varying N fertilizer levels (i.e., LN addition and HN addition) on each species' aboveground biomass in 2011 and in 2017 using separate one-way ANOVAs. We examined the impacts of high-level N and P additions on each species' aboveground biomass in 2011 and in 2017 using separate two-way ANOVAs, with high N addition and high P addition as main effects and an N × P interaction, and we used Student's t tests to determine the effects of either N or P alone when there were no significant interaction effects. Biomass data were natural log transformed for statistical analyses because the assumptions of constant variance and normality were violated, and a Wilcoxon test was performed when the assumptions for the t test were not met even after transformation. We used R software v3.2.3 (R Core Team 2015) for all statistical analyses. All statistically significant results (p < .05) and trends (p < .10) for all analyses are reported directly in the text.

Short-and long-term effects of the summer greenhouse treatment on plant species' aboveground biomass
Perhaps surprising, the shoot biomass of most vascular species significantly increased over twelve years at similar rates in both the greenhouse and the control plots (Tables 1 and 2, Figure 1, Appendix S4). These biomass increases generally happened during the second six-year period but not in the first six years of the study. A significant greenhouse effect was only observed for V. vitis-idaea, for which warming increased its shoot biomass in the period up to 2011 and then again up to 2017 (Figure 1d). By contrast, there was a positive interaction between the greenhouse effect and time for V. uliginosum due to an almost doubling of biomass in the second six-year period in the greenhouse plots and little response in the control plots (Figure 1b). Moreover, when examined at the whole vascular plant community level, the community biomass remained unchanged by 2011 (with a negligible 2.79 g/m 2 increase) but increased 69.15 g/m 2 (36 percent) in the greenhouses by 2017 (t = 1.99, p = .06; Figure 1g). In addition, total lichen community biomass decreased 24.81 g/m 2 (33 percent) in the greenhouses in 2011 (t = −1.82, p = .08; Figure  1h); and these declines continued up to 2017 (total decrease 64.57 g/m 2 [73 percent]; t = −4.40, p < .01; Figure 1h). Likewise, the greenhouses decreased the total moss

Effects of the summer greenhouse treatment on soil nutrient availability
Our IEM data showed that there were no long-term greenhouse effects on soil NH 4 -N fluxes, either for the whole growing season or for individual incubation periods within the growing season (Appendix S5). By contrast, the greenhouse treatment increased soil PO 4 -P fluxes by a factor of four (F 1,5 = 9.2, p < .05) for the whole growing season and by a factor of five in the fourth incubation period in particular (F 1,18 = 4.8, p < .05; Appendix S5).

Temporal changes in plant species' aboveground biomass of the unmanipulated control plots and in ambient air temperature
We used the complete data from the twenty-five control plots for the various experimental manipulations that are widely distributed across the whole study valley to more robustly determine whether and how shoot biomass of the investigated species changed with time. Shoot biomass of all vascular species in the ambient plots increased by 1.3-to 1.7-fold over twelve years without major shifts in relative abundance (Figure 2; Appendix S6). Moreover, these increases occurred only over the second six-year period, whereas biomass of these species was stable in the preceding period ( Figure 2; Appendix S6). For the nonvascular groups, lichen community biomass decreased in the first six-year period and then increased to the initial level during the second six-year period, whereas moss community biomass only showed a decrease in the second six-year period ( Figure 2; Appendix S6). Late spring and late summer air temperatures at our study site were generally significantly and substantially warmer over the period 2012-2017 compared to 2006-2011 and to the previous six years prior to our first sampling ( Figure 3). Specifically, mean monthly diel air temperatures were 2.9°C higher in May, 2.3°C higher in June, and 1.7°C higher in August in the period from 2012 to 2017 compared to the previous six-year period (Figure 3).

Effects of different N addition levels
Shoot biomass responses to the low and high levels of N fertilization addition differed among species and between short-versus long-term manipulations (Figure 4, Appendix S7). In the first six years, LN addition increased the biomass of V. uliginosum (t = 6.70, p < .01; Figure 4b) but had no effects on any of the other five species (Figures 4a, 4c-4f), whereas HN addition decreased the biomass of V. vitis-idaea (decreased by 78 percent, t = −7.47, p < .01; Figure 4d) and tended to decrease the biomass of R. subarcticum (decreased by 44 percent, t = −2.58, p = .06; Figure 4c) but had no effects on the other species (Figures 4a-4b, 4e-4f). Similar patterns were observed in 2017, after the full twelve years of manipulations (Figures 4g-4l): LN addition had no significant impacts, whereas the HN addition reduced biomass of R. subarcticum (by 74 percent, t = −3.93, p < .01; Figure 4i) and V. vitis-idaea (decreased by 88 percent, t = −11.33, p < .01; Figure 4j).

Effects of factorial high N and P additions
Species' biomass responses to HN and HP additions were also different among individual species and over the short versus long term (Figure 4, Appendix S7). In the first six years, a significant interactive effect between N and P additions was only observed for E. vaginatum, for which biomass increased 10.4 times in the HNHP plots (F 1,16 = 7.50, p = .01; Table 2, Figure 4f), but was unaffected by either  Monthly mean diel air temperature for late spring to early fall at Daring Lake over the three time periods of the point framing sampling (2000-2005, 2006-2011, and 2012-2017). Different labels (a), (b), and (c) indicate significant differences among time periods for each month at level of α = 0.05 (n = 6). Bolded lines indicate the median; boxes represent the 25th and 75th percentiles, and whiskers represent the 0th and 100th percentiles (excluding outliers). Black solid dots indicate outliers, which are defined as values more than 1.5 times the interquartile range either below the 25th quartile or above the 75th quartile. was negatively affected by both HN (decreased by 78 percent, t = −9.30, p < .01) and HP (decreased by 29 percent, t = −3.25, p = .01) additions, especially by N (Figure 4d). In contrast to these two evergreens, growth of the deciduous shrub V. uliginosum was more P limited. Its biomass was not affected by the HN addition but was stimulated by the HP addition by 7.4 fold (t = 6.08, p < .01; Figure 4b).
In the longer term (i.e., over twelve years), a negative interactive effect between N and P additions was observed for V. vitis-idaea (for which biomass decreased by 6.6 times; F 1,16 = 15.03, p < .01), primarily as a result of elevated N (Figure 4j), but not for the other species (Table 2,  . t Tests comparing the individual N and P treatments with the control indicated that HN reduced biomass of R. subarcticum (decreased by 74 percent, t = −3.58, p < .01) but tended to enhance biomass of B. glandulosa (increased twofold, t = 2.34, p = .05; Figure 4g). By contrast, no statistically significant P effects were observed (Figures 4g-4l). The HNHP treatment tended to increase graminoid species biomass very substantially (twelvefold), although there was very large variation in response among replicate plots (p = .09; Figure 4l).

Short-and long-term snow fence effects on plant species' aboveground biomass
No interactive effects between snow fence treatment and time were observed for any of the investigated species (Table 2). Species' shoot biomass generally remained unchanged in the first six years in both control and snow fence plots, except for a decrease for R. subarcticum in the control plots and an increase for E. vaginatum in the snow fence plots (Table 1). By contrast, the three dominant shrubs (R. subarcticum, V. vitis-idaea, and B. glandulosa) increased in the control but not in the snow fence plots in the second six years (Table 1).
The snow fence treatment generally had no effects on individual species in either sampling year, which may indicate the general lack of influence of changes in snow depth on soil nutrient availability as supported by our IEM NH 4 -N and PO 4 -P flux data (Appendix S5). In terms of responses of the nonvascular plant groups, the snow fences decreased lichen biomass by 21.43 g/m 2 (24 percent) and increased moss biomass by 68.61 g/m 2 (190 percent) in 2011 (t = −2.10, p = .07 and t = 1.95, p = .09, respectively; Appendix S8). Lichen biomass had decreased by 77.27 g/m 2 (49 percent; t = −2.97, p = .02; Appendix S8) in the snow-fenced plots in 2017, but there were no significant effects on moss biomass.

Short-and long-term caribou exclosure effects on plant species' aboveground biomass
Shoot biomass of all species remained unchanged in the first six years in both control and exclosure plots (Table  1). By contrast, the biomass of R. subarcticum, V. vitisidaea, and B. glandulosa in the control plots and of R. subarcticum, V. vitis-idaea, and V. uliginosum in the exclosure plots increased in the second six years ( Table  1). The exclosure treatment, however, did not affect shoot biomass of any of the six vascular species in either 2011 or 2017 (Table 2). No interactive effects between the exclosure treatment and time were observed for any of the investigated vascular species (Table 2). However, the exclosures tended to increase lichen community biomass (by 41.98 g/m 2 ; 50 percent) from 2005 to 2017 (t = 2.20, p = .09), whereas it remained unchanged in the control plots during the same period (Appendix S9).

Individual vascular species' biomass were generally unaffected by experimental greenhouse warming but did increase over time under ambient conditions, simultaneous with rising air temperature trends
Examining shoot biomass responses to both experimental summer warming using greenhouses and to natural ambient air temperature warming trends simultaneously within a single study allows us to better assess plant species' responses to warming, because the integration of these two approaches overcomes the methodological, spatial, and temporal limitations of each approach when adopted separately (Dunne et al. 2004;Pieper et al. 2011;Wolkovich et al. 2012). In contrast to our expectations, of any of the individual vascular species the greenhouse experiment did not enhance aboveground biomass of any of the vascular species (except for V. vitis-idaea) even after twelve years of manipulation. However, the consistent increases in biomass across the combination of all six vascular plant species in 2017 that led to a statistically significant greenhouse effect at the whole vascular plant community level indicate that this cumulative process is slowly occurring. On the other hand, we observed substantial biomass increases in all six species between 2011 and 2017 in the twenty-five control plots at our site, and these increases were coincident with warmer late spring and late summer air temperatures over that time period compared to the previous six years (Figures 2 and 3). This temperature increase (by 1.7°C-2.9°C) is within the higher range of what has been observed over the tundra biome between 1980 and 2010 (range = −1.5°C to 2.3°C; . Though there is a long history of significant results from experimental manipulations designed to mimic climate warming in tundra systems (e.g., Jonasson et al. 1996Jonasson et al. , 1999Wahren, Walker, and Bret-Harte 2005), few studies have shown significant responses at the level of plant growth and community structure in ambient plots that correlated to warmer air temperatures and possibly longer growing seasons.
The lack of a greenhouse treatment effect on individual species response may be caused by a mismatch in the seasonal timing when air temperatures were enhanced. The greenhouse plastic was installed from late June (or early July) to late August (or early September) each year, and this treatment time span only partly overlapped with the seasonal period of the natural ambient air temperature increases, which were primarily in May and June (Figure 3). Elevated ambient May-June air temperatures between 2012 and 2017 may have promoted species biomass by increasing soil temperature and thus the availabilities of soil nutrients, enhancing plant root capacities for nutrient uptake (Iversen et al. 2015), as well as by advancing the onset and therefore total length of the growing season. Previous studies of soil solution samples under ambient environmental conditions that were collected nine times over six months from early April to early September at our site indicate strong seasonal patterns in the soil biogeochemical properties Buckeridge et al. 2013). Dissolved inorganic N, inorganic P, organic N, and organic P were highest in May, declined in June, and then remained relatively stable throughout July and August Buckeridge et al. 2013). The trend of warmer May and June air temperatures since 2011 that we observed in this study may have enhanced these soluble soil nutrient peaks and thus enhanced spring to early summer plant growth. Indeed, other studies have shown that tundra vegetation is especially responsive to raised temperatures in early summer, in both aboveground and belowground growth (Macias-Fauria et al. 2012;Iversen et al. 2015).
Besides a mismatch in the seasonal timing when air temperatures were enhanced, the greenhouse treatment did not enhance mid-growing season soil nutrient availability to plants, at least for N. The absence of a greenhouse effect on soil available N in the data collected throughout the growing season in 2012 (Zamin, Bret-Harte, and Grogan 2014) and in our 2017 IEM data reported here supported this conclusion. However, the greenhouse treatment did increase overall mean summer soil PO 4 -P availability both in the earlier study (Zamin, Bret-Harte, and Grogan 2014) and in the IEM data set of this study. The extra PO 4 -P in the greenhouse plots may have been taken up primarily by the graminoid species E. vaginatum, which was much more responsive to P fertilization over the full study duration than the shrubs, suggesting that its growth was more strongly P limited and consistent with the conclusion of a previous study (Jonasson and Chapin 1991). Furthermore, it is noteworthy that although substantial variability precluded statistical significance, mean E. vaginatum biomass was more than doubled in the greenhouse plots compared to the controls by 2017 (p = .11), suggesting that the P limitation on its growth was alleviated by the greenhouse treatment. By contrast, and consistent with the general lack of P fertilization effects on the shrubs, the greenhouse-enhanced soil PO 4 -P did not stimulate growth of any of the other vascular species.
The absence of greenhouse effects in the 2011 data reported here conflicts with the results of another study conducted in the same treatment plots in the same year (Zamin, Bret-Harte, and Grogan 2014). The latter study used different 1 m 2 point frame areas that were directly adjacent to each of the permanent point frame areas used in our study. Zamin, Bret-Harte, and Grogan (2014) reported that the greenhouse treatment increased aboveground biomass of R. subarcticum (1.9 times) and of B. glandulosa (2.6 times). The disparities between the two studies holds even after the plant hits data in the earlier study (Zamin, Bret-Harte, and Grogan 2014) were recalibrated into biomass using the same calibration equations adopted in this study (data not shown). One possible cause may be subjective methodological issues related to the point frame method, but the fact that some species were not affected by the greenhouse treatment in either data set is at least a partial verification that the method was robust. We conclude that the most likely reason for the disparities is local spatial variation in vegetation composition and relative abundances of individual species within the treatment and control areas. Indeed, the relatively high species-level spatial variation compared to the whole vascular plant community responses to greenhouse reported above also supports this conclusion.
Lastly, the biomass of lichens declined with greenhouse treatment in both 2011 and 2017, consistent with several observations among tundra ecosystems (Chapin et al. 1995;Press et al. 1998;Walker et al. 2006;Hudson and Henry 2010). Increases in vascular plants and declines in lichens are well-recognized responses to experimental warming using open-topped chambers in the tundra biome and are attributed to increased height and density of higher plants resulting in decreased biomass of shadeintolerant lichens (Walker et al. 2006). In addition, in our greenhouse-based study, the exclusion of precipitation, which is the primary source of water and nutrients for lichens (Brodo, Sharnoff, and Sharnoff 2001), may partly explain the lichen declines we observed.

Tundra evergreen shrub species were more resistant to soil P increase than to soil N increase, whereas the dominant deciduous shrub (B. glandulosa) was more affected by the changes in soil N:P ratio
It is well recognized that large increases in soil nutrients due to high-level fertilizer additions are detrimental to tundra evergreen species but beneficial to deciduous species and graminoids because of the contrasting life history strategies between these plant growth forms (e.g., Chapin et al. 1995;Van Wijk et al. 2004;Grime 2006;Zamin, Bret-Harte, and Grogan 2014). Our results from the factorial N and P addition experiments are consistent with this conclusion. Two additional aspects of our N × P factorial results are novel: the relative importance of soil N versus soil P availability in determining evergreen shrubs' growth rates and the soil N:P ratio in determining the deciduous shrub growth responses.
In the first six years, the HN addition treatment led to a 233-fold increase in soil NH 4 -N:PO 4 -P ratio (caused by a 250-fold increase in soil salt-extractable NH 4 -N and a negligible change in soil salt-extractable PO 4 -P concentrations), whereas the HP addition led to a 3-fold decrease in soil NH 4 -N:PO 4 -P ratio in 2011 (caused by a 12-fold increase in soil salt-extractable NH 4 -N and a 35-fold increase in soil salt-extractable PO 4 -P concentrations), based on Zamin, Bret-Harte, and Grogan's (2014) study. This magnitude of change in available N and P in the soil solution of the HN addition versus the HP addition plots dramatically changed after another five years of manipulation. Specifically, the HN addition led to an 11-fold increase in soil NH 4 -N:PO 4 -P ratio (caused by a 68fold increase in soil distilled water-extractable NH 4 -N and a 7-fold increase in soil distilled water-extractable PO 4 -P), whereas HP addition led to a 214-fold decrease in soil NH 4 -N:PO 4 -P ratio in 2016 (caused by a 2-fold increase in soil distilled water-extractable NH 4 -N and a 366-fold increase in soil distilled waterextractable PO 4 -P concentrations; Gu and Grogan 2020). The stronger biomass declines in the dominant evergreen shrub species (R. subarcticum and V. vitisidaea) of the HN plots compared to the HP plots suggests that these species have a much lower tolerance for increased soil N availability than for increased soil P availability. On the contrary, biomass of the dominant deciduous shrub B. glandulosa seems to have been more affected by the magnitude of change in soil solution N:P ratio, supported by the fact that its biomass responded positively to the HN addition in 2016 that caused a much lessened imbalance in soil N:P ratio, whereas it was not responsive to the fertilization manipulations (i.e., HN addition in 2011 and HP addition in 2016) that caused a severe imbalance in soil N:P ratio, even though increased soil N or P availability individually is anticipated to increase deciduous species' biomass.

Experimentally deepened snow depth had negligible biomass impacts on vascular plants but reduced lichens
Our long-term snow fence treatment did not alter biomass for the shrub or graminoid species but substantially decreased biomass of lichens, which are more sensitive to environmental changes (Joly, Jandt, and Klein 2009). Although the deepened snow did not change soil nutrient availability, it led to wetter soils throughout November to July (Christiansen 2016) and shorter growing seasons (Nobrega and Grogan 2008;. These changes were detrimental to lichen growth, especially for moistureintolerant species such as Cetraria nivalis and Cetraria cucullata (Scott and Rouse 1995). Therefore, a heavier winter precipitation predicted under the current climate change scenarios (Wahren, Walker, and Bret-Harte 2005) may result in substantial decreases or even exclusion of lichens (Koerner 1980;Scott and Rouse 1995). In contrast to lichens, moss biomass was stimulated by the deepened snow cover, presumably due to the moderate increases in soil moisture, consistent with several observations from boreal systems (Rasmus, Lundell, and Saarinen 2011;Kreyling, Haei, and Laudon 2012). However, the statistical significance associated with this stimulation disappeared after another six years of manipulation, primarily due to a dramatic biomass decline in one snow fence plot that had subsided and flooded, causing a severely anaerobic environment due to the pooling of water that was clearly detrimental to moss growth.

Caribou exclosure had no biomass impacts on vascular plants but increased lichens over the decadal timescale
Our long-term caribou exclosure treatment did not alter biomass for the shrub and graminoid species but did significantly increase lichen biomass on the decadal timescale. Another study at Toolik Lake, Alaska, also found that caribou exclosures increased lichen biomass after seventeen years of manipulation (Gough et al. 2008). Deciduous shrubs are heavily selected by caribou during summer (White and Trudell 1980), whereas lichens are a dominant forage in spring, fall, and winter when the deciduous leaves are absent (Boertje 1984;Jefferies, Klein, and Shaver 1994;Saperstein 1996). The lack of biomass changes in B. glandulosa and an increase in lichens, therefore, may indicate negligible browsing of birch during the summer when its leaves are present but a stronger browsing of lichen when more caribous pass by during the nongrowing season. Indeed, our study valley lies close to the center of the Bathurst caribou (Rangifer tarandus groenlandicus) herd's movement between the spring calving and postcalving summer and winter ranges (Adamczewski et al. 2019). However, the Bathurst herd has experienced dramatic declines over the past several decades, from an estimated 470,000 in 1986 to 32,000 in 2009 and most recently to 8,200 in 2018 (Adamczewski et al. 2019). Birch shrub expansion within some landscape patches and birch encroachment across caribou trail pathways illustrated by repeat photography indicate evidence of this caribou herd decline at Daring Lake (Andruko, Danby, and Grogan 2020). Personal observations at our study site also confirmed that only small groups of caribou were regularly observed from 2004 to 2009 (Zamin 2013), and just a few caribou passed through our research valley each summer from 2015 to 2017 (Q. Gu, personal observation). The much longer time it takes for lichens to regenerate than for shrub foliage may also partly explain the different responses between lichens and shrubs since establishment of our exclosure.

Conclusion
Our study contributes three important conclusions to projections of future shifts in low Arctic tundra plant community structure in response to climate change. First, our results support the increased vegetation greening and shrub expansion that have been widely reported for arctic coastal and near-coastal environments (Bhatt et al. 2010;Ju and Masek 2016). Results from our site strongly suggest that similar responses are now becoming more evident in the continental interior of low Arctic Canada. The shoot biomass of each of the six vascular species that are widespread across the low Arctic biome and that represent three fundamentally different growth forms was unchanged between 2005 and 2011 but increased substantially between 2011 and 2017. Furthermore, these recent increases were correlated with simultaneous increases in late spring to early fall air temperatures at our site. Meanwhile, the greenhouses had a positive effect on aboveground biomass at the community scale but not at the individual species level due to high intraspecies spatial variability in vegetation composition. Experimentally enhanced snow depth or exclusion of caribou browsing had no impacts on vascular plant species biomass, although the latter did enhance lichen biomass. We conclude that our central continental low Arctic site is now undergoing the climate warming and associated vegetation changes that have been ongoing at other arctic locations for at least the past two decades. Secondly, our study highlights the relative importance of soil N versus P availability and of the soil solution N:P ratio in determining biomass responses of individual shrub species. We conclude that increases in soil nutrient conditions due to climate change may lead to species shifts within plant communities, which will be determined by the magnitude and the relative impacts of warming on soil N and P availabilities. Lastly, differences in sensitivity and direction among both the vascular and the non-vascular plant species (particularly lichens) in response to environmental changes indicate the importance of different niche space within the ecosystem in maintaining biodiversity in a rapidly changing world.