Modeling dynamics of tundra plant communities on the Yamal Peninsula, Russia, in response to climate change and grazing pressure

Understanding the responses of the arctic tundra biome to a changing climate requires knowledge of the complex interactions among the climate, soils and biological system. This study investigates the individual and interaction effects of climate change and reindeer grazing across a variety of climate zones and soil texture types on tundra vegetation community dynamics using an arctic vegetation model that incorporates the reindeer diet, where grazing is a function of both foliar nitrogen concentration and reindeer forage preference. We found that grazing is important, in addition to the latitudinal climate gradient, in controlling tundra plant community composition, explaining about 13% of the total variance in model simulations for all arctic tundra subzones. The decrease in biomass of lichen, deciduous shrub and graminoid plant functional types caused by grazing is potentially dampened by climate warming. Moss biomass had a nonlinear response to increased grazing intensity, and such responses were stronger when warming was present. Our results suggest that evergreen shrubs may benefit from increased grazing intensity due to their low palatability, yet a growth rate sensitivity analysis suggests that changes in nutrient uptake rates may result in different shrub responses to grazing pressure. Heavy grazing caused plant communities to shift from shrub tundra toward moss, graminoid-dominated tundra in subzones C and D when evergreen shrub growth rates were decreased in the model. The response of moss, lichen and forbs to warming varied across the different subzones. Initial vegetation responses to climate change during transient warming are different from the long term equilibrium responses due to shifts in the controlling mechanisms (nutrient limitation versus competition) within tundra plant communities.


Introduction
The Arctic is a complex system with strong interconnectedness among system components , Post et al 2009. Climate and vegetation interactions have been studied extensively.
From a vegetation perspective, the effects of climate change on such a complex system may be challenging to predict, as shifts in tundra plant communities affect ecosystem processes including net primary production, nutrient cycling and trophic interactions (Walker et al 2006, Post et al 2009. Repeat photography has documented expansion of shrubs, particularly deciduous shrubs such as Salix and Alnus spp. over the past 50 years in northern Alaska, likely caused by a warming of the arctic climate (Tape et al 2006). Assessing the responses of tundra vegetation at the plant functional type level is challenging using remote sensing (which includes aerial photography) and has most often been accomplished with field experiments. Post and Pedersen (2008) revealed that graminoid-dominated tundra shifted to dwarf birch-dominated tundra in five years within enclosures with passive warming using open-top chambers (OTCs) in the inland area of Kangerlussuaq Fjord, West Greenland. In metaanalyses of tundra warming experiments, warming increased height and cover of deciduous shrubs and graminoids and decreased cover of mosses and lichens, however individual studies show that the responses differed among sites (Henry and Molau 1997, Arft et al 1999, Walker et al 2006. Grazing exists in most ecosystems and, with concomitant trampling, constitutes a regular disturbance regime to which some plants are more or less adapted (Milchunas et al 1988). The reindeer (Rangifer tarandus) is a ruminant of the family Cervidae, which is distributed circumpolarly and has long been a key component of high latitude ecosystems (Forbes and Kumpula 2009). Reindeer grazing and trampling effects on arctic plant species have been studied relatively extensively (Bråthen and Oksanen 2001, Moen and Danell 2003, Olofsson et al 2004, Bråthen et al 2007, Gough et al 2007, Pajunen et al 2008, Kitti et al 2009. However, grazing interactions with climate warming and the effects on plant community structure have received less attention (Post andPedersen 2008, Olofsson et al 2009). Post and Pedersen (2008) found that grazing mitigated long term warming effects on tundra plant communities, as was also suggested in a study where herbivores were found to inhibit shrub expansion (Olofsson et al 2009). These results are consistent with a recent modeling study which found that grazing caused tundra community biomass to decline, and that grazing interacted with warming to buffer plant responses (Yu et al 2009). In contrast, when erect willows (Salix lanata), that would otherwise constitute an important source of forage for reindeer in areas of heavy grazing, are already above the limit of browsing (≈180 cm), their growth increases significantly in response to decadal warming (Forbes et al 2010). Further investigation of the responses of the major tundra plant functional types (PFTs) is needed to improve understanding of warming and grazing interactions on tundra plant community biomass and composition.
Different plant species may respond individualistically to changes in temperature, nutrients, and herbivory, and these responses may differ over time (Hollister et al 2005). Short term experiments, however, are not able to provide an understanding of how tundra plant communities may reorganize over long term trajectories. Simulation modeling can be used to examine multiple factors collectively and provide insights to the functionality of arctic systems across a broad range of scenarios (Yu et al 2009). Our objective is to investigate the effects of the tundra climate gradient, climate warming, soil nutrients and grazing on tundra vegetation dynamics at the PFT level. Our research questions are: (1) how does climate warming affect PFTs across the five latitudinal subzones (Walker et al 2005) of the Arctic tundra biome (i.e. from polar desert to Low Arctic)? (2) How does grazing affect each PFT across this latitudinal gradient? (3) How do grazing and climate warming effects differ for PFTs along the latitudinal gradient, and how do these factors interact?

Study region
The Yamal Peninsula, in northwest Siberia, Russia, is a region with many environmental impacts including nomadic reindeer herding by indigenous Nenets, gas development and variable climatic change (Forbes and Kumpula 2009, Walker et al 2009b. Our study region stretches from southern Low Arctic tundra at Laborovaya on the Peninsula (67 • 42 N) to the polar desert at the Krenkel station (80 • 38 N) on Hayes Island, Franz Josef Land, and consists of study locations along a transect spanning the five arctic bioclimate subzones (figure 1).
Yamal soils consist of nutrient poor marine sands and nutrient-rich clays. The loamy sites tend to have greater soil organic nitrogen contents than the sandy sites. Comparable sandy and loamy sites were established in each location, and different vegetation cover is associated with each site (table 1). On mainland Yamal, we have study sites in subzones C, D and E located near Kharasavey, Vaskiny Dachi and Laborovaya respectively. For the loamy sites, dwarf-shrub, moss tundra is often dominant, while for sandy sites, prostrate dwarf-shrub, lichen tundra dominates. Study sites north of the Yamal include Ostrov Belyy (subzone B) and Krenkel station (subzone A). Ostrov Belyy is an island north of the Yamal Peninsula and includes two study sites: one on a mesic loamy landscape (OB-1) and the other on a drier sandy landscape (OB-2). The dominant plant communities at OB-1 include graminoid, prostrate dwarf-shrub, moss tundra on inner-circle areas and prostrate dwarf-shrub, crustose-lichen barrenson non-sorted circles. The communities at OB-2 are moss, prostrate dwarf-shrub tundra. These sites are used for summer and autumn pasture for migrating reindeer. Two study sites were established at Krenkel station on Hayes Island (subzone A) with one zonal sandy loam site (KR-1) and one drier sandy loam site (KR-2) in 2010. These sites are within the polar desert geobotanical subregion, and the vegetation is characterized by lichens, mosses, and cushion forbs at KR-1, and less moss but more cryptogamic crusts at KR-2.
More detailed site descriptions and sampling methods can be found in the data reports of the 2007-2010 Yamal expeditions (Walker et al 2009a(Walker et al , 2009c(Walker et al , 2011, and several published studies associated with the sites along the Yamal transect (Goetz et al 2011, Walker et al 2009b, Yu et al 2009, Kumpula et al 2011, Leibman et al 2011.

ArcVeg model
ArcVeg is a nutrient-based dynamic vegetation model (see Epstein et al 2000 for model description) that has been applied extensively for North American tundra (Epstein et al 2001(Epstein et al , 2004(Epstein et al , 2007. Arctic tundra vegetation in the model is grouped into 12 plant functional types including mosses (MOSS), lichen (LICH), forbs (FORB), tussock sedges (TUSS), nontussock sedges (SEDG), rushes (RUSH), grasses (GRAS), evergreen prostrate dwarf shrubs (<5 cm, EPDS), evergreen erect dwarf shrubs (15-40 cm, EEDS), deciduous prostrate dwarf shrubs (<5 cm, DPDS), deciduous erect dwarf shrubs (15-40 cm, DEDS) and low shrubs (40-200 cm, LOWS). The model essentially runs with nitrogen mass balance equations, moving nitrogen among three major pools (total soil nitrogen, plant available nitrogen and plant biomass nitrogen by PFT) (Epstein et al 2007). Climate and grazing control nitrogen fluxes among these pools (see figure 2). The twelve plant functional types compete for the limited nitrogen and use the acquired nitrogen to generate new biomass. Climate subzone, soil organic nitrogen levels and grazing regime dictate the simulated plant functional type composition and biomass for a particular tundra system. The model is parameterized for the five arctic bioclimatic subzones that range from the polar desert (subzone A) to the Low Arctic tundra at the southern extent of the tundra biome (subzone E) (Walker et al 2005). Soil organic nitrogen levels were taken from field data for each site, and extrapolated to account available nitrogen in active layer. Climate warming in the model alters the available nitrogen for plants, assuming greater mineralization rates with higher temperatures (Epstein et al 2000, Cornelissen et al 2007. Additionally, the length of the growing season in the models is affected by climate change.

Model adjustment for selective grazing
Grazing, as a major disturbance in ArcVeg, is characterized using two parameters: frequency of grazing and percentage of biomass eaten. For example, a grazing regime of (0.1, 25%) means that there is 0.1 probability of a grazing herd each year, and a maximum of 25% biomass is removed with grazing. Grazing is selective as herbivores generally prefer forage with a high content of nutrients and a low level of structural and chemical defenses (Hanley 1997). To study plant-herbivore interactions, diet selection is important for inclusion in model simulations. Reindeer are ruminants with a highly selective forage strategy to adapt to cold weather, short growing seasons and small plants. Grazing selectivity considered in ArcVeg is a function of both foliar nitrogen concentration and reindeer diet preference. The diet preference used in this letter is based on studies by the Reindeer Research Program in Fairbanks, Alaska, and reflects the selectivity of reindeer diet in general (Bucki et al 2004, Flenniken 2007. The data show a high  reliance on lichen with a seasonal spike in use of vascular plants during the growing season. Salix spp., Eriophorum spp., and Carex spp., make a big proportion of reindeer diet during the growing season (Podkorytov 1995, Kitti et al 2006, Forbes and Kumpula 2009).

Simulation design and model input data
Our model simulations were driven by bioclimatic subzones and field collected soil organic nitrogen (calculated based on N%, bulk density and active layer depth) (table 2). Input data include soil organic nitrogen (SON), climate subzone, grazing regimes ((0.1, 25%), (0.1, 50%), (0.5, 25%) and (0.5, 50%)) and climate warming manipulations (before warming, during warming (transient), and after warming (equilibrium)). All parameter combinations (including eleven sites-which span the five subzones and have different SON levels, four grazing regimes and three climate scenarios) were simulated with 20 replicate runs for each scenario. The three climate change scenarios included control/before warming, transient warming and equilibrium/after warming, and the four grazing regimes included: (0.1, 25%), (0.1, 50%), (0.5, 25%) and (0.5, 50%). The (0.1, 25%) scenario was assumed to be the control grazing regime or low grazing, and (0.5, 50%) was the highest grazing intensity, indicating reindeer herds will graze on the same site every two years, and each visit a maximum of 50% total biomass will be removed. These are reasonable grazing scenarios on the Yamal Peninsula, where the Nenets may set up camp, and herds will graze essentially around the camp for several days. The Nenets normally migrate across the peninsula every year with designated migration routes for each brigade (Stammler 2005), but there are also a large number of privately owned herds grazing on the same territories . Although there are no managed reindeer herds at Krenkel station on Hayes Island, for this study, we simulated grazing effects in subzone A for zonal vegetation dynamics comparison. Plant functional type biomass values were compared and evaluated across all sites on the Yamal Peninsula under the different grazing regimes with one subzone climate warming (essentially 2 • C). In total we simulated for 1500 years. At year 1000, the 2 • C warming was initiated and ramped linearly for 50 years; after that, the climate stayed at the warmer state for another 450 years. We calculated the mean biomass for each plant functional type for 100 years before warming (year 901-1000), for 100 years during the transient warming period (year 1001-1100), and for 100 years after warming (year 1401-1500).

Statistical analysis
Nonmetric multi-dimensional scaling (NMS) was used to assess simulated biomass of 12 plant functional types and their relationships to the environmental controls including subzone, soil organic nitrogen, grazing, and climate change. NMS was developed by Kruskal (1964) and is widely used for ecological applications (Urban et al 2002). It is based on ecological distance and essentially produces a low-dimension ordination space in which sample separation in this space reflects sample separation in a multi-dimensional 'species' space. The algorithm ranks correlation between ordination distance and ecological distance. NMS makes no assumptions about the nature of species response to underlying gradients, and it is well suited for non-normal data (McCune et al 2002). Any similarity or dissimilarity (distance) matrix (e.g. Jaccard, Bray-Curtis or Euclidean) can be used. A total of 132 ArcVeg simulated plant communities (as a multivariate of PFT biomass values) and their relationships to environmental variables (latitudinal climate gradient or five climate subzones, eleven soil organic nitrogen levels, three temporal climate scenarios and four grazing intensities) were analyzed. We used PC-Ord (version 5.0) for NMS with a Bray-Curtis distance matrix for this study. The NMS was run randomly with 250 runs with real and randomized data respectively. The final stress was 6.7 for a three-dimension solution, and the instability was 10 −5 , which suggests a good ordination with no real risk of drawing false inferences (McCune et al 2002).
The simulated biomass values of 12 PFTs were grouped into 6 growth forms (moss (MOSS), lichen (LICH), forbs (FORB), graminoids (SEDG, TUSS, RUSH and GRAS), evergreen shrubs (EPDS and EEDS), and deciduous shrubs (DPDS, DEDS and LOWS)) and were analyzed as dependent variables using multi-factor ANOVA (analysis of variance). The main effects were climate subzone, SON, warming, and grazing (all categorical variables except for SON). The interactions of the main effects were also included in the ANOVA. Least square means and Type III sum of squares were used to account for the unbalanced data. This analysis was performed in SAS version 10.0 for Windows (SAS institute Inc.).

Overall results in ordination
NMS ordination was applied to assess controls on simulated tundra plant community properties (figure 3). Each point in figure 3 represents one model scenario, and the colors indicate the five subzones. We used the biplot function in PC-Ord to show the direction and strength of correlations of environmental variables (figure 3(a)) within the ordination space, which are shown with the red arrows. In figure 3(a) the vertical axis is strongly correlated with the latitudinal temperature gradient (∼12 • C) and summer warmth (42% of total variance in model simulations for all arctic tundra subzones), the horizontal axis is strongly correlated with the grazing gradient (13% of total variance), and the rest of the variance is explained by SON, soil texture, warming and the interaction effects. Sites at Krenkel and Ostrov Belyy (subzones A and B) formed one cluster in the north of the ordination space and sites at Kharasavey, Vaskiny Dachi and Laboravaya (subzones C, D and E) clustered together probably because they have similar community properties which are determined by SON and SWI. The biplot in figure 3(b) illustrates the absolute value of PFT biomass. Most plant functional types, except evergreen shrubs, were negatively affected by grazing, especially lichen, forbs, rushes, grasses and deciduous prostrate dwarf shrubs ( figure 3(b)). Relative abundance of each PFT for each simulation was calculated and shown in figure 3(c). Relative abundances of mosses, lichen, forbs, rushes and grasses were greater to the north, in the polar desert (northern High Arctic), while relative abundances of tussock sedges, non-tussock sedges, evergreen and deciduous shrubs were the greatest in the southern subzones. Relative abundance of evergreen shrubs (figure 3(c)) correlated positively with grazing, suggesting minimal grazing impact on evergreen shrubs. Total biomass in figure 3(d) increased toward the low latitudes or south subzones (represented by the size of each dot, which increased toward south), but decreased toward the high grazing intensity.
Multi-factor ANOVA was conducted to test the individual and interaction effects and the results are presented in the supplemental materials (table A.1 available at stacks. iop.org/ERL/6/045505/mmedia).
Each individual effect (subzone, SON, warming and grazing) and interaction effects (subzone*SON, subzone*warming) on our grouped plant functional types were significant ( p < 0.01).

Grazing effects
Overall, grazing caused total plant community biomass to decrease, but PFTs were affected differently by grazing ( figure 4). Lichen, forb, graminoid, and deciduous shrub biomass declined in response to increased grazing frequency and percentage. Lichen and deciduous shrubs were affected the greatest by grazing. For mainland Yamal (subzones C, D and E) before warming, on average deciduous shrub biomass declined from 249 g m −2 when grazing was (0.1, 25%) to 198 g m −2 when grazing was (0.1, 50%), and to 111 g m −2 when grazing was (0.5, 25%) and 17 g m −2 when grazing was (0.5, 50%). Evergreen shrub biomass increased as grazing intensity increased, and the increase was most substantial (about 29 g m −2 ) when grazing was most intense (0.5, 50%). Moss biomass responded to grazing nonlinearly, increasing 1% (before warming) and 2% (after warming (equilibrium)) when grazing percentage increased from 25% to 50% and then declining (−11% and −13%) when grazing frequency increased from 0.1 to 0.5 (every ten years to every two years); this response was consistent under transient warming and equilibrium warming scenarios.

Temporal climate warming effects on tundra plant communities
Climate warming affected PFTs differently across our study sites (figure 5). Shrubs responded to warming with increased biomass at all study sites. The responses of moss, lichen and forbs to warming varied across the different subzones. In addition, some PFT responses during the transient warming period were different from their responses following the transient warming period.

Warming period effect (transient versus equilibrium warming)
. Different PFTs responded to transient warming and equilibrium warming differently. Evergreen and deciduous shrubs responded to both warming scenarios with continued positive responses, and these responses were consistent across both latitudinal and soil nutrient gradients, although the magnitude varied slightly. Evergreen and deciduous shrub biomass increased from 64 g m −2 and 93 g m −2 to 100 g m −2 (54% increase) and 132 g m −2 (42% increase) respectively during transient warming, and to 126 g m −2 (96% increase) and 188 g m −2 (102% increase) during equilibrium warming. Graminoids also had a positive response to warming at all sites but responded to transient warming with greater biomass increase (13 g m −2 or 26% increase) than equilibrium warming (10 g m −2 or 19% increase). Lichens responded to transient warming with only a 16% biomass increase and had a 3% biomass decrease during equilibrium warming when grazing was (0.1, 25%). Forbs responded with greater biomass increase during transient warming (15%) than during equilibrium warming (7%), although these responses differed across the latitudinal gradient and with different soil nutrient levels.

Interaction with latitudinal gradient and soil properties.
In general, warming caused most PFTs to increase biomass, but this response interacted with the latitudinal gradient and soil properties including soil texture and soil organic nitrogen levels. Multi-factor ANOVA showed that interaction effects between subzone and warming on most plant types were significant ( p < 0.01, table A.1 available at stacks.iop. org/ERL/6/045505/mmedia). For example, moss responded to warming positively in subzones A (136% increase) and B (171% increase), and negatively in subzone E (52% decrease). In subzones C and D, moss biomass increased during transient warming but decreased during the equilibrium warming period. Lichen responded to warming positively with increased biomass in subzone A, but the response to warming was negative in subzones B, C, D and E. Essentially, warming affected moss and lichen negatively where shrubs became largely dominant. Graminoids had a positive response to warming at all sites but responded during transient warming with a greater biomass increase than equilibrium warming in subzones A (40% versus 22% increase), B (28% versus 16% increase), and D (16% versus 7% increase), yet had a smaller biomass increase during transient warming than during equilibrium warming in subzones C (35% versus 41% increase) and E (11% versus 13% increase). The interaction effect of SON and warming significantly affect deciduous and evergreen shrubs ( p < 0.05, table A.1 available at stacks. iop.org/ERL/6/045505/mmedia). Across all sites, deciduous and evergreen shrubs responded to warming consistently, but with less biomass increase in sandy sites or low SON sites (85 and 33 g m −2 ) than in loamy sites or high SON sites (139 and 53 g m −2 ).

Warming and grazing interactions
Compared to climate change, grazing had a more substantial effect on plant communities ( figure 7(a)). We found that plant responses to grazing were different during the control climate, transient and equilibrium warming, indicating interactions of climate change with grazing. The interaction effects between subzone and grazing, and warming and grazing on deciduous shrubs and lichen were significant ( p < 0.01, table A.1 available at stacks.iop.org/ERL/6/045505/mmedia). Deciduous shrub, lichen, forb, and graminoid biomass decreased continually with increased grazing intensities. However, since warming may promote greater plant growth, such a decline due to grazing is mediated by climate warming, particularly for deciduous shrubs already above the browse limit for reindeer.

Discussion
According to the presented data and analyses, we found that grazing can be important in addition to the latitudinal temperature gradient (∼12 • C) in controlling tundra plant community properties. The NMS results of this study showed that grazing, as one of the controlling factors in addition to climate change (2 • C) and soil organic nitrogen, explained about 13% of total data variance in model simulations for all arctic tundra subzones, while latitude explained about 42%. This can have important implications for reindeer husbandry across the Arctic. The interactions of warming and grazing are potentially complicating our understanding of tundra vegetation dynamics and need to be taken into account in future research efforts.

Grazing effects
Grazing may cause total plant community biomass to decrease as indicated in our previous modeling study (Yu et al 2009). But, understanding how major plant functional types may respond, especially the types that compose the reindeer diet, is critical for projecting future changes in tundra vegetation cover affected by grazing and trampling. Our results suggest that as lichen and deciduous shrubs are preferred by reindeer, they can be the plant functional types that are most affected by grazing management patterns. This finding is in line with remote sensing data from along the Finnish-Norwegian border where a reindeer fence has been visible due to higher reindeer lichen coverage in Norway; the fence prevented lichen being grazed by reindeer in Norway in comparison to Finland (Kumpula 2006, Forbes andKumpula 2009). Study conducted in Brøggerhalvøya, Svalbard also found that lichen and vascular plants declined with increased grazing (Staaland et al 1993). Moss, which is not preferred as much by reindeer (Flenniken 2007), showed a nonlinear response to increased grazing in High Arctic. For example, the biomass in subzone C increased by 6% when grazing scenarios shifted from (0.1, 25%) to (0.1, 50%) and then declined by 12% and 19% respectively when grazing intensity continued to increase to the (0.5, 25%) and (0.5, 50%) scenarios. Increased grazing may favor moss growth by reducing species competition and changing the soil moisture regime, possibly causing the tundra plant community to shift from lichen-dominated tundra to moss-dominated (van der Wal 2006). But with continual increase in grazing, moss may also be foraged (Staaland et al 1993) and negatively affected by trampling.
Grazing interacts with soil nutrients, which can contribute to more complicated tundra plant responses to warming. For example, forb biomass generally decreased with increased grazing intensity, but there are exceptions in nutrient poor sites or sandy sites (KH-2, VD-2 and LA-2). In these sites, forbs responded positively when grazing increased from (0.1, 25%) to (0.1, 50%) and negatively with further grazing increases. This may due to the fact that in nutrient poor sites, grazing affected deciduous shrubs substantially and thus there may be nutrients available to other plant functional types such as forbs.
Evergreen shrubs in our current model simulations respond to increased grazing intensity with continual increase of biomass. This is due to the fact that they are least preferred by reindeer in the diet data used by ArcVeg. Moderate grazing has been shown to cause evergreen shrubs to increase albeit not always (Olofsson et al 2004, Bråthen et al 2007. The increase of evergreen shrubs in this model study may not represent the situation that we see in the field on the Yamal ; one reason may be that the susceptibility of vegetation to reindeer trampling has not been taken into account in the model. Additionally, the rapid cycling of nutrients in the model over time through grazing may not be accurately represented in the model as reindeer expedite decomposition and relocate nutrients through feces and urine (Olofsson et al 2001). While trampling during intense grazing can affect all taxa, and cause both deciduous and evergreen shrubs to be damaged , deciduous shrubs may have an advantage over some other plant types in that they have higher photosynthetic rates and higher reproductive capacity as mechanisms of resilience to disturbance. The more common prostrate deciduous shrubs on the Yamal (e.g., Salix polaris, S. arctica, S. nummularia) are fairly resilient to grazing (Walker et al 2009a(Walker et al , 2009c.

Resilience and palatability
We tested how the resilience of evergreen shrubs affected their response to grazing to help us understand how growth rates may affect the responses of evergreen shrubs to different grazing regimes. Evergreen shrubs are less palatable in comparison to deciduous shrubs due to their low nitrogen concentration and their poor digestibility, which contribute to their low percentage in reindeer diet throughout the year. Evergreen shrubs also have slower growth rates than deciduous shrubs in our model; however this did not prevent evergreen shrubs from becoming dominant with increased grazing intensity. Based on observations from the Yamal, deciduous shrubs still have greater abundance than evergreen shrubs on the Yamal, which is a heavily grazed system compared to other places in northern Russia and North America (Walker et al 2009a). We examined this phenomenon by altering evergreen shrub growth rates (nitrogen uptake: biomass ratio) to see how responses differed under the four grazing regimes.
The response of evergreen shrubs to grazing changed along a small gradient of evergreen shrub growth rates under different grazing scenarios.
Slower growth rates limited evergreen shrubs, and thus there was substantially less biomass when growth rates decreased for evergreen shrubs from 0.030 g N uptake g biomass −1 to 0.020 g N uptake g biomass −1 , and the response to increased grazing became negative in subzones D and E (figure 6). Moreover, as the growth rates of evergreen shrubs decreased, the negative effect of increased grazing pressure on deciduous shrubs became reduced.
Heavy grazing may cause plant communities to shift. When evergreen shrub growth rate was reduced to 0.020 g N uptake g biomass −1 , the proportional abundance of moss and graminoids increased with continual increase in grazing pressure. Figure 7(b) shows the plant community transition from shrub-dominated tundra toward moss, graminoiddominated tundra when grazing intensity increased. We tested this grazing effect with before warming biomass data for shrubs (both evergreen and deciduous shrubs) using a singlefactor ANOVA (table A.2 available at stacks.iop.org/ERL/ 6/045505/mmedia). When evergreen shrub nutrient uptake rate was 0.02 g N uptake g biomass −1 , heavy grazing caused subzones C and D to shift significantly from shrub-dominated to graminoid-moss-dominated tundra. Decrease in shrubs may be beneficial for understory plants, since there will likely be more light and nutrients available. Heavy trampling however may negatively affect moss growth and reduce the moss layer depth (Van der Wal and Brooker 2004), but this is not included in current model version. Altered or expedited nutrient cycling by reindeer digestion appears most advantageous for graminoids, according to field nutrient addition manipulations (Olofsson et al 2001. Tundra plant communities tend to be altered by grazing and trampling, and have been suggested to shift toward graminoid-dominated systems under these conditions (Zimov et al 1995). Transition toward graminoid-dominated tundra caused by heavy grazing has been observed on the Yamal at the local scale .

Warming effects
Shrub expansion in the Arctic, presumably caused by warming, has been captured by studies using repeat aerial photographs (Sturm et al 2001, Tape et al 2006. Other remotely sensed imagery such as from AVHRR has documented greening trends in the tundra region (Jia et al 2003, Goetz et al 2005, Bhatt et al 2010. Shrub dendrochronology has also been used in the Arctic (Forbes et al 2010, Hallinger et al 2010. Forbes et al Figure 6. Resilience test of evergreen shrubs (evergreen shrub growth rates included 0.030 g N uptake g biomass −1 , 0.025 g N uptake g biomass −1 , and 0.020 g N uptake g biomass −1 ) shown in the left column. Deciduous shrub biomass was compared without change of deciduous shrub growth rates (right column).

Figure 7.
Proportional abundance of the main functional groups across all sites. (a) Illustrates how warming affects the responses of PFTs to grazing (BW-before warming, TW-transient warming and EW-equilibrium warming). Proportional abundance was calculated when evergreen shrub growth rate was 0.030 g N uptake g biomass −1 . (b) Illustrates the change in dominant plant types caused by grazing when evergreen shrub growth rate changed from 0.030 g N uptake g biomass −1 to 0.020 g N uptake g biomass −1 .
(2010) explored the chronology of Salix lanata L. (sensu lato), an abundant erect deciduous shrub in the Yamal region and an important source of reindeer forage throughout the Low Arctic. They found a significant increase in shrub willow growth over the last six decades, and a strong correlation with the normalized difference vegetation index (NDVI), corroborating the 'greening' trend detected by remote sensing studies. Our model results also suggest that shrubs may respond to warming with a continual biomass increase, at least under a consistent grazing regime (figure 5). This may be highly applicable to Alaskan tundra where caribou grazing dominates, and grazing pressure is generally low. This study expands our understanding of shrub growth under different grazing regimes, including the Yamal Peninsula, Russia, where large herds of reindeer managed by the Nenets can graze and trample extensively at relatively high levels and, in many places, quite intensively.
Our model-simulated warming responses are consistent with field studies, such as experimental warming studies from the International Tundra Experiment (ITEX) using open-top chambers (OTCs), which can increase mean growing season air temperature by 1-3 • C. Our modeled warming was a 2 • C increase linearly ramped over a 50 year period, thus the rate of temperature increase was not as great as in the ITEX studies (Henry and Molau 1997), however we simulate for a longer period of time. Moss and lichen biomass measured in ITEXassociated studies declined in response to warming in Low Arctic sites, while deciduous shrubs and graminoids increased (Chapin et al 1995, Henry and Molau 1997, Walker et al 2006. This is consistent with our modeling results in the Low Arctic. Moreover, our results provide insights as to how these plants respond differently in soils of different organic nitrogen contents. Limitation of soil organic nitrogen may constrain plant responses to warming (Henry and Molau 1997). Additionally, graminoids and forbs responded with greater biomass increase on average during transient warming (26% and 14%) than during equilibrium warming (19% and 7%), and this may be due to shifts in controlling mechanisms from direct warming response to species competition for nutrients (Epstein et al 2000). Species competition, essentially for limited nutrients, may cause some PFTs to decline in biomass after transient warming. Such responses cannot be detected with short term experiments (Chapin et al 1995). Shifts in controlling mechanisms need to be considered to understand and predict how arctic systems will respond over the long term to changing climate.

Interaction effects
PFT responses to grazing varied under the three climate scenarios. As grazing may cause most PFTs to decline in biomass, more substantially at higher intensity grazing regimes, warming interacts with grazing and may contribute to more complicated plant responses. Lichen and deciduous shrubs were significantly affected by warming and grazing interactions ( p < 0.001, table A.1 available at stacks.iop. org/ERL/6/045505/mmedia). In general, grazing negated plant biomass increases in response to warming. Shrub biomass increased under warming scenarios alone, yet responses were buffered when grazing pressure increased (figure 7). Deciduous shrubs responded to warming with increased biomass, and the increases were profound at lower intensity grazing regimes. This is consistent with what has been found in northern Fennoscandia, where the abundance of the dominant shrub (Betula nana) has increased, and the increase has been more pronounced when grazing was absent from the study site (Olofsson et al 2009). The responses to grazing under different climate change conditions were different. Graminoid and forb biomass declined with increased grazing, and this decrease became greater under the equilibrium climate state due to enhanced shrub growth and thus competition for nutrients despite that warming can accelerate nitrogen mineralization and provide more nutrients for vegetation. Large scale surveys in regions of Finnmark and northern Norway with reindeer herds of high densities found that high grazing substantially reduced large dicotyledons and grasses in fertile sites (Bråthen et al 2007). Long term grazing, trampling and feces deposition may cause the plant community state to shift from woodydominated tundra to clonal rhizomatous graminoids . Although our model provides understanding of long term tundra vegetation dynamics, the current version does not include either trampling or rigorous nutrient recycling routine for reindeer waste and thus need further improvement and investigation to better understand plant community state shifts.

Conclusion
Simulation results in this study suggest that grazing can be important in addition to the latitudinal temperature gradient (12 • C) for tundra plant communities, explaining about 13% of the total variance. Plant functional types (PFTs) such as lichen, deciduous shrubs and graminoids responded to grazing with decreasing biomass. However, such a decline is potentially mediated by climate warming since generally warming promotes the growth of shrubs and graminoids, particularly when erect deciduous shrubs are above the browse line. Moss biomass had a nonlinear response to grazing, and such responses were stronger when warming was present. Our results suggest that evergreen shrubs may benefit from increased grazing intensity due to their low palatability, yet a growth rate sensitivity analysis suggests that changes in nutrient uptake rates may result in different shrub responses to grazing pressure. Heavy grazing may cause plant communities to shift from shrub tundra toward moss, graminoid-dominated tundra in subzones C and D. Further analyses with inclusion of trampling effects are strongly needed for better understanding of the interactions between shrub dominance and grazing. In response to climate warming alone, moss, lichen and forbs varied across the different subzones. Deciduous and evergreen shrubs responded positively to warming consistently across all subzones, but with less biomass increase in low SON sites than in high SON sites. Initial vegetation responses to climate change during transient warming are different from the long term equilibrium responses due to shifts in the controlling mechanisms (nutrient limitation versus competition) on tundra plant communities.