Herbivore Effects on Ecosystem Process Rates in a Low-Productive System

Mammalian herbivores shape the structure and function of many nutrient-limited or low-productive terrestrial ecosystems through modification of plant communities and plant–soil feedbacks. In the tundra biome, mammalian herbivores may both accelerate and decelerate plant biomass growth, microbial activity and nutrient cycling, that is, ecosystem process rates. Selective foraging and associated declines of palatable species are known to be major drivers of plant–soil feedbacks. However, declines in dominant plants of low palatability often linked with high herbivore densities may also modify ecosystem process rates, yet have received little attention. We present data from an island experiment with a 10-year vole density manipulation, to test the hypothesis that herbivores accelerate process rates by decreasing the relative abundance of poorly palatable plants to palatable ones. We measured plant species abundances and community composition, nitrogen contents of green plant tissues and multiple soil and litter variables under high and low vole density. Corroborating our hypothesis, periodic high vole density increased ecosystem process rates in low-productive tundra. High vole density was associated with both increasing relative abundance of palatable forbs over unpalatable evergreen dwarf shrubs and higher plant N content both at species and at community level. Changes in plant community composition, in turn, explained variation in microbial activity in litter and soil inorganic nutrient availability. We propose a new conceptual model with two distinct vole–plant–soil feedback pathways. Voles may drive local plant–soil feedbacks that either increase or decrease ecosystem process rates, in turn promoting heterogeneity in vegetation and soils across tundra landscapes.


INTRODUCTION
Mammalian herbivores exert strong control on vegetation biomass and species composition in many low-productive and nutrient-limited ecosystems, such as the northern tundra biome (Oksanen and others 1981;Aunapuu and others 2008). In parallel with changes in vegetation and leaf chemistry, herbivory may both accelerate and decelerate ecosystem process rates, i.e., rates of primary production, litter decomposition and nutrient cycling in tundra (Olofsson and others 2004b;Brå then and others 2007), boreal (Pastor and Naiman 1992) and grassland systems (Sirotnak and Huntly 2000;others 2004, 2009). However, it remains an open question if herbivoredriven acceleration of process rates can occur in nutrient-limited communities, and if so, under what conditions (Jefferies and others 1994;Augustine and McNaughton 1998;Bardgett and Wardle 2003).
Herbivore-induced change in process rates in nutrient-limited systems depends strongly on changes in plant and litter quality and litter input quantity (Bardgett and Wardle 2003;Cherif and Loreau 2013;Sitters and Olde Venterink 2015), at least at intermediate soil moisture levels (Schrama and others 2013). Direction and magnitude of process rate change are contingent on whether grazing-tolerant, palatable plants with high N content and low C:N ratio (e.g., forbs) will increase or decrease in abundance in comparison with nonpalatable, grazing-intolerant and N-poor plants with a high C:N ratio (for example, evergreen dwarf shrubs) (Pastor and Naiman 1992;Augustine and McNaughton 1998;Mulder 1999). This is because palatability and N content are linked with process rates such as plant growth rates and litter decomposability (Jefferies and others 1994;Cornelissen and others 2004). Plant-herbivore dynamics are traditionally thought to be driven by the effects of selective herbivory on palatable, fastgrowing plants (Grellmann 2002; Bardgett and Wardle 2003;Brå then and others 2007;Bakker and others 2009). More specifically, when plant growth is already limited by nutrient availability, selective feeding by relatively few herbivores reduces the abundance of fast-growing species further, which favors slow-growing species (Bardgett and Wardle 2003;Brå then and others 2007). Slow-growing species slow down process rates further by modifying plant-litter-soil feedbacks (Bardgett and Wardle 2003;Weintraub and Schimel 2005;Cornwell and others 2008), if the decreased N content at plant community level and in litter (cf. Garnier and others 2004) concomitantly decreases soil N availability and microbial activity (Cornwell and others 2008;Sitters and Olde Venterink 2015). In contrast, urine and feces deposition may enhance nutrient cycling others 2004, 2009) locally in low-productive systems (Stark and others 2015).
However, especially when present in high densities, herbivores may directly affect slowly growing non-palatable species (Hairston and others 1960). Changes in vegetation biomass and composition during population peaks of voles and lemmings, key herbivores across northern tundra systems (Ims and Fuglei 2005), represent a well-documented example of such effects (Ruffino and others 2015). Although voles feed selectively on fast-growing forbs and bilberry others 2013, 2014 ,  Table 1), the physical disturbances associated with high vole densities induce massive landscape-scale declines in evergreen dwarf shrubs others 2012, 2014;Ruffino and others 2015;Hoset and others 2017). In line with this, during the first three years of the long-term experimental study presented here, we documented strong suppression of dominant dwarf shrubs and an increase in forbs under high grey-sided vole (Myodes rufocanus) density (Dahlgren and others 2009). Declines in non-palatable species contradict patterns expected based on selective herbivory and classical models of herbivore-plant-soil interactions (Pastor and Naiman 1992; Bardgett and Wardle 2003, Table 1) and result from, for example, clearing paths during winter, shrub shoot clipping and (non-selective) grazing of unpalatable plants under extreme resource depletion (Jefferies and others 1994;Korslund and Steen 2006;Soininen and others 2013). Despite these strong and well-documented effects on vegetation, only a couple of studies have addressed the effects of voles on soil processes and plant-soil interactions in the arctic (Stark and Grellmann 2002) or elsewhere (Sirotnak and Huntly 2000;Bakker and others 2004).
Notwithstanding a consensus that modification of plant and litter quality and quantity are central drivers of herbivore effects on soil nutrient cycling and microbial activity, empirical findings often fail to corroborate the classical theory in low-productive systems (Sitters and Olde Venterink 2015;Stark and others 2015). Recent work has helped to account for these discrepancies, especially in lowproductive systems through consideration of stoichiometric mismatches between herbivores, plants and detritus, herbivore identity and soil physical properties (Bakker and others 2009;Cherif and Loreau 2013;Schrama and others 2013). However, we know little of how herbivore-induced increases or decreases in non-palatable plants affect plant and soil nutrient contents. Herbivore effects on non-palatable plants, especially under higher herbivore densities (vs. low desities, see Grellmann 2002;Stark and Grellmann 2002), are rarely explicitly considered as drivers of ecosystem process rates.
Here we test the complementary hypotheses that (1) cyclic vole populations with high-amplitude peaks accelerate ecosystem process rates in a nutrient-limited system by reducing the abundance of slow-growing species, and that (2) experimental reduction in vole density in a nutrient-limited system decelerates process rates as decreased disturbance favors slow-growing species as superior competitors. More specifically, we predict that 1. High vole density negatively affects dwarf shrubs and induces a shift in plant community composition toward increased abundance of forbs (compare Dahlgren and others 2009).  Dahlgren and others (2009) and Cornelissen and others (2004). 2 Soininen and others (2013). 3 Cornelissen and others (2004). 4 Kleyer and others (2008). 5 Olofsson and others (2004a) and Dahlgren and others (2009). 6 Hoset and others (2017).
Under experimentally reduced vole density, dwarf shrubs will become more abundant relative to forbs. 2. High vole density and induced shift in plant community composition result in higher plant N content in individual plants (hereafter species-N) and in community-weighted average (hereafter community-N, Garnier and others 2004), larger content of available soil inorganic N (hereafter soil N or NH 4 + -N content), and higher microbial activity in litter and soil in comparison with low vole density.

Study Area Characteristics and Habitat Description
The study was conducted on the islands of Lake Iešjá vri (68.16 km 2 ; 69°45¢N, 24°30¢E) in West Finnmark, Norway, between 2000 and 2013. The region has a weakly continental subarctic climate (Oksanen and Virtanen 1995) with average July and January mean temperatures of 12°C and -13°C and precipitation of 61 mm and 22 mm, respectively (2004-2014 averages from Suolovuopmi station, www.eklima.no). Areas surrounding the lake are characterized by mosaics of different types of dwarf shrub tundra and lowshrub wetland (Walker and others 2005) with patches of mountain birch woodland in depressions and in sites with warm microclimate. The most common vascular plant species in the dwarf shrub tundra are deciduous Betula nana, Vaccinium myrtillus, V. uliginosum and Salix spp., evergreen Empetrum nigrum ssp. hermaphroditum, V. vitis-idaea, Andromeda polifolia and V. microcarpum and small forbs Rubus chamaemorus, Cornus suecica L. (Chamaepericlymenum suecicum (L.) Ascherson and Graebner), Stellaria borealis and Linnea borealis. Graminoids such as Carex sp., Calamagrostis lapponica and Eriophorum angustifolium are present, but in very low abundances. The lake is adjacent to Joatka Research Area, where local vole and lemming population dynamics has been monitored since 1986. The main mammalian herbivores in the study area are grey-sided voles (Myodes rufocanus), Norwegian lemmings (Lemmus lemmus) and reindeer (Rangifer tarandus ssp. tarandus). Grey-sided voles exhibit regular 5-year cycles, whereas lemming dynamics follow a 10-year pseudoperiodic cycle with peaks roughly coinciding with vole population peaks (Ekerholm and others 2001). Reindeer migrate through the areas surrounding the lake during spring and autumn, but rarely visit the small islands (L. Oksanen, personal observation).
The experimental islands of Lake Iešjá vri are relatively small and isolated. Vegetation on the islands follows a zonal pattern from moist and lush vegetation dominated by willow thickets along the shoreline, to wind-exposed lichen heaths in the center up to a few meters above the lake water level ( Figure 1A). Above the willow thickets, each island harbors mesic dwarf shrub vegetation with distinct peat ridges running parallel to the shoreline ( Figure 1B) and with a vegetation community consisting mainly of evergreen and deciduous ericoid dwarf shrubs, dwarf birch and small forbs (Dahlgren and others 2009). Relative abundance of dwarf shrub species differs between exposed ridges and moister depressions (see arrows in Figure 1B). Based on trapping data from Joatka Research Area, this vegetation type constitutes an important habitat for grey-sided voles (Ekerholm and others 2001;Hoset and others, unpublished).

Study Design
We used two complementary sampling designs to test our predictions. First, we used a transplant experiment, where vegetation blocks from a large island with a moderate vole density were transplanted onto islands with manipulated vole densities that were either higher or lower than in the origin island. Second, we surveyed transects to analyze the vegetation composition in two adjacent microhabitats, ridges and depressions on the transplant islands, to assess to what extent the small-scale spatial heterogeneity of dwarf shrub tundra modulates the effects of herbivory on plant communities (Jefferies and others 1994).
The transplant experiment consists in total of nine islands of varying size, of which eight (0.02-4.2 ha, see Dahlgren and others 2009) were remote from the mainland and assigned as vole density manipulation islands. The remaining, largest island (105 ha) is the origin of transplanted vegetation blocks, but is excluded from further analysis due to abiotic and biotic differences in comparison with smaller islands and lack of replication. The experimental islands have low vole immigration rates and are effectively predator-free (Dahlgren and others 2009). On four of the islands, vole populations with cyclic fluctuations have been maintained since 1991 by introducing new individuals during springtime to mitigate effects of winter mortality on the populations (hereafter called ''high density''). The other four islands were kept nearly vole-free by recurring, mostly annual, trapping (hereafter ''low density''). One of the low-density islands was excluded from the current analysis due to divergent topography and ground moisture regime and a history of high vole density. The transplant experiment was started in July 2000 by excavating mesic heath vegetation blocks (70 by 70 cm, depth > 30 cm) from the origin island (see method details in Dahlgren and others 2009). Eight blocks (hereafter transplants) were transplanted onto each of the eight experimental islands along peat ridges, above the area prone to disturbance by ice scouring and avoiding wet depressions.

Herbivore Data
To monitor the success of the high-density treatment, vole live trapping was conducted on the high-density experimental islands during early summer (late June-early July) and/or late summer (late August-early September) in 2000-2003, 2007, 2010 and 2011; during other years voles were absent. Altogether, the vole density manipulation covered three peaks and two low-density phases, that is, 2.5 cycles ( Supplementary Figure 1). We used Ugglan live traps (Grahnab, Swe-den) in a grid with trapping stations at 10-m intervals across the islands, but excluding the most barren mid-ridge if present ( Figure 1A). The vole abundance estimates represent the number of uniquely trapped individuals, representative of the minimum number of individuals known alive (MNKA), corrected for island size to obtain density estimates per hectare (Supplementary Figure 1). All abundance estimates were calculated using data from late summer trapping.
Vole density fluctuations on the high-density islands were similar to dynamics on the origin island and the adjacent mainland (Dahlgren and others 2009), with peak-year densities 3-4 times higher than on the origin island. Vole densities remained very low on the low-density islands (Supplementary Figure 1). Vole and lemming populations peaked exceptionally high in 2010-2011 across northern Norway (Olofsson and others 2012) and also led to a spontaneous immigration of lemmings to high-density islands. Thus, our treatments successfully maintained high and low vole densities and mimicked the natural vole and lemming population dynamics in the area.

Change in Time: Transplant Vegetation Data
Transplant vegetation on high and low-density islands was monitored at the start of the experiment in 2000 and again in 2010 and 2011 by point intercept frequency method (Brå then and Hagberg 2004). We used 100 sampling pins (Ø 2 mm) distributed across a 50-by-50-cm Plexiglas plate centered onto each transplant plot. The pins had a fixed, initially randomly determined placing within the sampling plate. All hits of all vascular species were counted up to a maximum of 20 hits per species per pin. We refer to Dahlgren and others (2009) for vegetation changes during the first years of the transplant experiment (that is, 2000-2003).
To obtain plant species-specific biomass estimates for all years (that is, 2000, 2010 and 2011 for transplants, 2013 for transects, see below), we calibrated point frequency scores and plant biomass. For this purpose (and subsequent leaf N content analysis, see below), on each island the aboveground biomass of maximum four of the transplants was harvested in July-August 2011. The harvested vascular biomass of each transplant was sorted to species in the field, and later stems of woody dwarf shrubs (excl. V. myrtillus, whose green stems were pooled together with leaves to form shoot biomass) were further separated from the leaves. Samples were first air dried in room temperature and then dried in 50°C for 24 h, after which species-specific dry mass was measured. Correlation coefficients between point frequency scores and harvested biomass were calculated using weighted linear regression (Supplementary Table 1a-b; Brå then and Hagberg 2004). These coefficients were used to transform point frequency scores from all plots and all years into biomass estimates used in further analysis.

Variation in Space: Transect Vegetation Data
To gain a spatially broader view of the vole impacts on plant communities, we sampled vegetation plots along transects spanning through the entire mesic hummock tundra habitat on each island in 2013 (hereafter transects). We established transects along peat ridges, with length and number depending on the size and continuity of the habitat. The total length of transects per island was 20-38 m, but on some of the islands this was covered by one continuous transect, whereas on others we used two shorter transects. Along a transect, sample plots were placed every two meters, with one plot sampled on top of a ridge and one in the adjacent depression ( Figure 2B) to address the potentially confounding effect of abiotic differences between these two microhabitat types. Maximum distance between the paired plots was set to 2 m approximately perpendicular to the transect. Transect plot size was 40 by 40 cm, and we sampled vegetation by point intercept frequency method with nine pins ( Figure 2C, D), which is sufficient for capturing community composition attributes in dwarf shrub tundra (Brå then and Hagberg 2004). Transect plots that (1) were on top of transplant plots, (2) clearly disturbed by human activity (for example, a path, skidoo track), (3) vole activity marks (on low-density islands) or (4) had a stone larger than 10 by 10 cm were excluded and replaced by a plot within 1 m in either direction along the transect. The point frequency scores were transformed to biomass estimates as described above.

Plant Species Data for Analysis
For all plant community analyses, we used estimated biomasses of the six most common vascular plant species: evergreen dwarf shrubs E. nigrum and V. vitis-idaea, deciduous dwarf shrubs B. nana and V. myrtillus, and forbs R. chamaemorus and C. suecica. Their combined biomass encompassed roughly 93-95% of the vascular biomass on the transplant plots between 2000 and 2011, and 80% of the vascular biomass in the transect plots. We used a ratio of evergreen dwarf shrubs to forbs for analyzing vegetation change in transplant and transect datasets, because these growth forms are strong indicators of plant community composition and soil process rate changes in the study system (Chapin and others 1996;Freschet and others 2010). This ratio was log e -transformed to gain a normal distribution of values, after which zero value denotes an even ratio, negative values indicate dominance of forbs and positive values indicate dominance of evergreen dwarf shrubs. The log e -transformed values were used in further analyses (hereafter EHr).

Plant N Content
We analyzed the plant nitrogen (N) content (species-N, % of dry weight) using near infrared reflectance spectroscopy (NIRS). The samples were taken from a subset of the harvested transplant plot biomass (96 samples approximately equally distributed between species and islands). Stems and leaves for forbs, and leaves for evergreen and deciduous dwarf shrubs, were included in the samples. For forbs, the inclusion of stems may re-duce the recorded species-N values compared to measurements based on leaves only. Each sample was homogenized and ground to a fine powder using a ball mill (Mixer Mill, MM301; Retsch GmbH & Co. Haan, Germany), after which the dry powder was pressed into Ø 16-mm tablets by applying approximately 5 tons of pressure. Tablets were scanned with a FieldSpec 3 (ASD Inc., Boulder, Colorado) in 350-2500-nm range. Sample N contents were estimated based on an arctic-alpine model (Ancin and others, unpublished manuscript) which showed good accuracy in the predictions (R 2 validation ¼ 0:94, root mean standard error of the prediction (RMSEP) = 0.20% in the range 0.34-6.01% N DW). Estimated N values were included within the calibration range.
Vegetation community-N was used to evaluate effects of vole density at the community level and calculated following the method of Garnier and others (2004) where p i is the proportion of total green biomass contributed by species i (shoot biomass for forbs and V. myrtillus, leaf biomass for other dwarf shrubs), n equals the six most common species included in the analysis and N is the species-N described above. Species-N applied in the calculation of community-N was independent of vole density treatments, as the species-N did not show statistically significant differences between treatments for all species (see results). Our community-N values are therefore conservative, representing only responses in plant community composition to the vole density treatments and not including species-level chemical trait responses. For community-N calculations, we assumed that relative differences in species-N does not vary between years, using N content in plant individuals sampled in 2013 to represent species-N from all years.

Transplant Soil and Litter Sampling and Analysis
To test the hypotheses that high vole densities are associated with enhanced soil inorganic nutrient contents and microbial activity, intact cores containing both the litter and the soil organic layers were sampled from transplants in August 2012. Eight cores (Ø 3 cm, maximum 30 cm depth) per plot were collected along the rims of all transplants (two cores on each side) in intact vegetation outside the area sampled for vegetation data, but inside the transplant. Exact depth of each core was measured in the field. The top 5 cm of each core (top soil) was separated from the rest (bottom soil) and pooled to form one composite topsoil sample and one composite bottom soil sample per transplant plot. The composite samples were stored in plastic bags, kept cool (< 10°C) and frozen within 30 h of the sampling.
After thawing, litter was separated from the top soil sample and soil samples were homogenized and roots removed by sieving (mesh 3 mm), after which soil was subsampled for different analyses. Soil and litter moisture was analyzed (105°C, overnight) and soil and litter organic matter content by loss on ignition (475°C, 4 h). Values were then used to calculate soil and litter organic matter content (OMC %) and soil organic matter stock per area (kg/m 2 ). To analyze soil inorganic nutrient contents, a subsample of $ 3 g fresh soil was extracted with 50 ml of 0.5 M K 2 SO 4 . The contents of NH 4 -N and PO 4 -P in the extracts were analyzed following established colorimetric methods (SFS To analyze microbial activity in the litter and soil, microbial respiration (total CO 2 release) was measured in the laboratory by incubating soil and litter samples equivalent of 1 g OM in 120-ml glass vials. After a 48-h pre-incubation, the CO 2 concentrations in the headspace were analyzed twice during a 92-h incubation period in + 9°C (Stark and others 2015). Vials were aerated after the pre-incubation and between the CO 2 measurements. The soils were incubated in field moisture, which was 68.4% ± 0.53 (mean ± s.e.) of fresh soil and 59.14% ± 1.14 for litter and did not vary significantly by treatment.

Statistical Analyses
Due to the spatially and temporally nested structure of the different datasets, we used linear mixed effects modeling (Pinheiro and Bates 2009) with island as a random intercept for analyzing plant and soil responses to vole density treatment. Homogeneity of variances was confirmed visually by quantile-quantile plots and by plotting residuals against fitted values and response variable values (Zuur and others 2010). We selected models based on an information-theoretic approach using smallsample corrected Akaike's information criterion (AICc, Burnham and Anderson 2002). A priori determined candidate models were compared based on differences in their AICc values (DAICc), Akaike weights (Wt) and parameter estimates. In most cases, DAICc values were low, suggesting little support for any single model being best. Consequently, we applied a model averaging procedure that included all models in the candidate set to take into account model selection uncertainty (Burnham and Anderson 2002).
To link variation in plant community composition and plant N content with vole density, we applied two sets of linear mixed effects models (lme) to transplant and transect data separately. In the first set of analyses, two main gradients of plant community composition were used as response variables. These were identified as the first two axes of principal component analysis (PCA, based on the correlation matrix) summarizing the species-wise biomass estimates. Since all the species were present in almost all the plots, the sampled environmental gradient was short enough for the assumption of linear species responses to be reasonable (Legendre and Legendre 2012). Transplant-PC1 and transplant-PC2 were analyzed with only vole density in the fixed term. Fixed term of the global model for transect-PC1 and transect-PC2 consisted of interactive effects of vole density (low vs. high) and micro-topographic location (ridge vs. depression), to which model averaging was applied. In addition, to obtain comparable predictions for plant community responses to vole density between transplants and transects, the transect-PC1 and transect-PC2 scores were analyzed with vole density as the only fixed term.
The second set of mixed effects analyses included EHr, biomass estimates of the six most common species, species-N and community-N as response variables. First, we analyzed vole density and micro-topography effects on EHr, biomass estimates of the dominant species and community-N on transects. The global model included interactive effects of vole density (low vs. high) and micro-topographic location (ridge vs. depression). All models derived from the global model formed the candidate model set, and model averaging was applied to the entire candidate set. Then, we analyzed the temporal change in EHr, dominant species, species-N and community-N in response to vole density on transplants from year 2000 to 2010 and 2011. In the lme-model, we included year as random slope in addition to island as random intercept. The fixed term included the interactive effect of vole density and year. Because this interaction was central to our hypothesis, no model selection or model averaging was necessary. Model fit was assessed visually and if necessary, a log e -transformation of the response variable was applied.
Soil and litter variables from the transplants-microbial respiration, soil ln(NH 4 + -N) and PO 4 --P contents, soil organic matter (OM) content and soil bulk density-were response variables in separate linear mixed models. Soil OM stock was analyzed for the whole soil sample to avoid measurement errors due to compaction during sampling, but for all other response variables topsoil and bottom soil were analyzed separately. EHr was included as a fixed covariate in the analysis to model the vegetation-mediated vole density effects on soil variables. Thus, the full model consisted of the interactive effects of the categorical vole density variable and covariate EHr, to allow vole density and plant community to have nonadditive effects on soil and litter. Because the sample size was larger for positive values of EHr (see Figure 5), which might increase the risk of a type I error, we compared the results from the model with full dataset with separate linear mixed models for high and low-density islands with only EHr as the explanatory variable. As this did not impact our inference, we report only model results based on the full dataset. We base our inference jointly on (1) estimates derived from best models indicated by DAICc, (2) graphical representation of observed values and (3) model selection tables, where lowdensity treatment was set as the intercept. In addition, we explored the link between vegetation and soil nitrogen by applying a linear mixed model with ln(NH 4 + -N) as the response variable and community-N as the fixed term.
All statistical analyses were conducted using R statistical software (version 3.3.1). Linear mixed models were conducted with the nlme-package (Pinheiro and others 2016), model averaging and selection using the dredge function in package MuMIn (Barton 2016) and PCA analyses using the prcomp function in package stats. Transplant vegetation plot 95% confidence levels of the mean were obtained by stat_summary (fun.data = ''mean_cl_boot'') function in package ggplot2 (Wickham 2009).

RESULTS
Vole density manipulation had strong effects on the log e -transformed evergreen dwarf shrub:forb biomass ratio (EHr), biomass estimates of the six dominant plant species, species-and community-N as well as soil inorganic nutrients and microbial activity in the litter layer. Transplant plot results demonstrate a strong causal link between variation in vole density and in plant community composition over time.

Vole Density Effects on Plant Community Composition
Principal component analysis (PCA) and subsequent linear mixed models on PC1 and PC2 scores indicated divergence in vegetation communities between high and low vole densities in transplants and transects (Figure 2A-C, Supplementary Table S3). For both transplants and transects, small values of PC1 were associated with high abundance of forbs C. suecica and R. chamaemorus, whereas large values were associated with high abundance of evergreen dwarf shrubs E. nigrum and V. vitis-idaea ( Figure 2B, C). In both transplants and transects, PC1s indicated a marginal association with vole density treatment (Supplementary table S3), which was likely due to high degree of variation in vegetation community composition on islands with high vole density ( Figure 2B, C). Although PC1s seemed to represent a transition between forb and evergreen dwarf shrub abundance, PC2s were associated with vegetation patterns related to micro-topography, for example, exposure and moisture, and subsequent variation in abundance of deciduous dwarf shrubs B. nana and V. myrtillus (Supplementary Table S3).

EHr and Plant Species Responses to Vole Density
A high correlation between EHr and transplant-PC1 (Pearson r = 0.904) confirms that both approaches capture the same main gradient of vegetation composition. On the transplants, EHr was affected by 11 years of vole density manipulation. On lowdensity islands, EHr remained close to the initial positive values indicating strong evergreen dwarf shrub dominance. In contrast, on the high-density islands the index declined toward and below zero, indicating dominance of forbs (by 2011 t (149) = -4.074, p £ 0.001; Supplementary Table 4, Figure 3A). High vole density was also associated with higher within-and between-island variation in EHr, as with plant community composition PCA (details not shown). EHr responded to vole density in a similar way on the transects. Although mean EHr did not differ significantly between vole density treatments (z = 1.387, p = 0.166), large effect sizes indicate clear trends toward forb dominance in many plots on high-density islands (Supplementary Table 4, Figure 3B).
Species-specific responses varied consistently within growth forms: evergreen dwarf shrubs declined and forb biomass increased under high vole density on transplants, but with no effect of vole density on deciduous dwarf shrubs. On transects, evergreen dwarf shrubs and forbs were predominantly affected by micro-topography and less by voles. V. myrtillus was virtually absent from ridges, but present in depressions. In depressions, V. myrtillus and B.nana biomass under high vole density was less than half of their biomass compared to low-density islands (Table S5b, Figure S2f

Plant N Responses to Vole Density
The plant N content in the studied species (species-N, % dry weight) varied consistently between species belonging to different growth forms (Figure 4A). Although the expected N contents in green plant tissue would follow the order forbs > deciduous dwarf shrubs > evergreen dwarf shrubs (Reich and others 2007), we found Herbivore Effects on Ecosystem Process Rates the highest species-N contents in deciduous dwarf shrubs, intermediate in forbs and lowest in evergreen dwarf shrubs. This result may be due to the inclusion of green stems of forbs in the analyzed samples. Vole density effect on species-N was only statistically significant for C. suecica ( Figure 4A, Supplementary Table 2a). However, there was a tendency for higher species-N under high vole density for E. nigrum, B. nana and V. myrtillus ( Figure 4A, Supplementary Table 2a).
In transects ( Figure 4C), community-N was lower on the ridge communities than in depressions (z = 8.629, p < 0.001), and community-N responses to vole density varied between microhabitats (Supplementary Table 2c). In the ridge habitat, high vole density was associated with higher community-N compared to low-density islands (z = 3.595, p < 0.001). In depressions, vole density treatment did not have an effect on community-N (z = 1.329, p = 0.184).

Soil and Litter Responses to Vole Density and EHr
Decreasing EHr values were linked with increasing ammonium content in topsoil (t = -4.730, p £ 0.001) and bottom soil (t = -2.561, p = 0.015), phosphate content in bottom soil (t = -2.637, p = 0.013) and litter microbial respiration (t = -3.017, p = 0.005). These relationships were indicated in best model estimates (singlefactor main effect models, Supplementary Table 6a-b) and visual examination of observed values ( Figure 5).
Neither vole density nor EHr had impacts on soil OM stock, topsoil and bottom soil microbial respiration, topsoil bulk density or litter, topsoil and bottom soil OM content (Supplementary Table 6a  p = 0.031) showed a statistically significant negative response to high vole density.

DISCUSSION
Our experimental study revealed that high peak vole densities over multiple population cycles increased the relative abundance of forbs, plant N content, soil inorganic N content and microbial activity in the litter layer. Together, the findings support our hypothesis that small herbivores can increase ecosystem process rates in inherently lowproductive systems by reducing the abundance of slow-growing plants with low palatability (Figure 6, arrows a2-c2).

Vole Activity Effects on Plant Community Composition
Low vole density for over a decade led to an increase in the abundance of dominant E. nigrum and to a decline in forb R. chamaemorus (see also Gough and others 2012), indicative of vegetation changes under resource competition. However, decrease in vole density did not lead to significant changes in the ratio of evergreen dwarf shrubs and forbs, counter to our predictions. In contrast, our high vole density treatment reduced the biomass of evergreen dwarf shrubs and led to increasingly forb-dominated plant communities, in line with our prediction. Importantly, this shift persisted over multiple population cycles, including two 3-4-year low-density phases. Severity of biomass removal and time to recover seem key determinants in breaking the dominance of evergreen dwarf shrubs, especially the allelopathic E. nigrum (Aerts 2010), which suppresses growth of other plants (Brå then and others 2010). As voles and lemmings can decimate locally up to 100% of dwarf shrub vegetation during population peaks (Olofsson and others 2012;Hoset and others 2017), population lows of 3-4 years may be too short for E. nigrum to recover (Aerts 2010). Increased nutrient availability via herbivore urine and excreta (Bakker and others 2004) may have contributed to the shift in plant community composition by facilitating the growth of subordinate forbs under severe disturbance to E. nigrum (Aerts 2010), although this mechanism may be more pronounced in communities with low plant C:N ratio (Sitters and Olde Venterink 2015). Notably, effects of herbivores on non-palatable plants may require periodic high herbivore densities, as E. nigrum recovers well after 50%, but not after 100% aboveground biomass removal (Aerts 2010). This may indicate a threshold of vole impact intensity that would allow for a vegetation shift toward faster process rates, especially under nutrient-enriched conditions (Aerts 2010).
Although evergreen dwarf shrubs declined and forbs increased consistently with vole density, micro-topography modulated both plant community composition and vole density effects on vegetation. The expected strong negative effect of voles on deciduous dwarf shrubs, and especially on the preferred forage V. myrtillus, was evident in depressions. Also several previous studies have showed that cyclic vole populations can suppress the growth of deciduous dwarf shrubs and therefore counteract the effects of climate warming-induced shrubification (Olofsson and others 2004a;Kaarlejä rvi and others 2015). However, the way herbivory interacts with climatic warming to determine future species compositions may be contingent on herbivory effects on E. nigrum and its ability to act as a niche constructor (Brå then and others 2018), on summer warming effects of E. nigrum's recovery rates (Aerts 2010) and how these effects interact with landscape micro-topographic variability (Graae and others 2018). Our results indicate that while high-intensity vole grazing reduces both deciduous and evergreen dwarf shrubs across a micro-topographically heterogeneous landscape, deciduous and evergreen dwarf shrubs experience strongest reductions in different parts of the micro-topographic gradient. Given the ecosystem-level consequences of shrubification on tundra functioning (Myers-Smith and others 2011), acknowledging this potential small-scale variation in herbivory control over vegetation changes is important for climate change predictions.

Plant Tissue N Content at Individual and Community Levels
Consistent with our prediction, high vole density was linked with elevated N levels in shoots of C. suecica and R. chamaemorus, and leaves of V. myrtillus, B. nana, V. vitis-idaea and E. nigrum, albeit significantly only for C. suecica, likely due to low sample size and high variation. Elevated species-N levels of deciduous and evergreen dwarf shrubs under high vole density contrast previous studies reporting negative, species-rather than functional group specific, or no effects of herbivores on plant nutrient content in slow systems (Eckstein and others 1998;Stark and Grellmann 2002;Zamin and Grogan 2013). Hence, the elevated species-N may be a response to long-term (as opposed to short term, Nordin and others 2004) increased soil N availability (Sirotnak and Huntly 2000;Strengbom and others 2003) observed in this study.
Consideration of plant N content as a community-aggregated parameter provides a link between individual eco-physiological and population-level responses to the ecosystem-level effects of herbivores (Eskelinen and others 2012) and constitutes a long-term feedback measure for N cycling and herbivore forage quality (Augustine and McNaughton 1998;Bardgett and Wardle 2003). We found evidence of herbivores reinforcing plantlitter-soil feedback for N cycling (Bardgett and Wardle 2003;Bakker and others 2009) with elevated community-N under high vole density on Figure 6. Conceptual model of herbivore effects on plant-soil feedbacks in a slow system. The model incorporates effects of herbivores on both preferred and non-preferred forage plants, and litter-mediated effects of plants on soil process rates and forage quality. Within these boundaries, we suggest two dynamic feedback loops, one driven by the effects of selective foraging (pathway 1) and the other driven by herbivore-induced decimation of non-preferred plants (pathway 2). Pathway 1 ''decreasing process rates'': Slowly growing plants with relatively high LDMC, SLA and high C:N ratio (Freschet and others 2010;Reich 2014), low palatability and low tolerance of mechanical disturbance (Augustine and McNaughton 1998), for example, evergreen dwarf shrubs, dominate low-productive tundra systems (Walker and others 2005). They escape low-intensity selective grazing and browsing targeting subdominant, palatable forage (Bardgett and Wardle 2003) and thus are favored by it (a1). These species produce slowly decomposable litter (b1), and hence, increase in slow species dominance leads to lower rates of nutrient cycling, microbial decomposition and plant productivity (c1; Chapin and others 1996;Reich 2014). Reduced availability of labile nutrients in soil further favors slow growth species (d1), and the reduced community-level forage quality further decreases herbivore activity (e1). Pathway 2 ''increasing process rates'': Decimation of slow-growing dwarf shrubs and increased availability of microsites (a2; Nystuen and others 2014) allow for faster growing forbs to increase in abundance. Higher leaf and whole plant N content, higher SLA and lower lignin content typical of forbs and grasses are linked with higher quantity and quality litter (b2) promoting faster soil process rates (c2; Cornelissen and others 2004;Cornwell and others 2008;Reich 2014). Alleviation of nutrient limitation further promotes abundance of faster growth species (d2), increased availability of higher-quality forage and consequently increased herbivore activity (e2; Augustine and McNaughton 1998). Inputs of feces and urine alleviate conditions for forbs and further accelerate nutrient cycling (compare Bakker and others 2004). ridges, and positive correlation between community-N and topsoil NH 4 -N contents, corroborating our predictions. Our results indicate that community-level trait responses are informative of changes in N cycling in slow herbivore-plant-soil systems.

Plant-Litter-Soil Dynamics Under High and Low Vole Density
Contents of available soil NH 4 -N, PO 4 -P and litter microbial respiration rate increased under high vole density. However, since increases in soil nutrient contents and litter respiration rate were associated with increased forb dominance rather than with high vole density per se, our results support the view that herbivores affect soil processes mainly via the plant-litter-soil pathway (Chapin and others 1996;Olofsson and Oksanen 2002;Cornelissen and others 2004). Consistent with this view, in a study where vole and reindeer grazing rather favoured plant species with slow growth rates (Grellmann 2002) the aboveground changes were accompanied by decelerated soil ecosystem process rates (Stark and Grellmann 2002). The increase in the content of available soil NH 4 -N in our study did, however, not lead to increased soil microbial respiration, which counters our initial predictions, but agrees with previous findings from low-productive tundra systems, where a grazer-induced increase in soil NH 4 -N availability is not always linked with enhanced CO 2 release from soil others 2002, 2015). This could indicate herbivore-induced C-limitation of microbes (Stark and Grellmann 2002;Cherif and Loreau 2013;Sitters and Olde Venterink 2015; but see Stark and others 2015).

Conceptual Model of Small Rodent-Plant-Soil Dynamics
The potential of voles and lemmings to decimate dominant dwarf shrubs is well documented (Olofsson and others 2004a, 2012; Kaarlejä rvi and others 2015; Hoset and others 2017). However, these effects have not been regarded as leading to longterm vegetation community changes in low-productive systems (Jefferies and others 1994;Augustine and McNaughton 1998) nor have they been empirically or theoretically linked with changes in plant-soil interactions. Small rodent population fluctuations, weakening of predator control over small rodents (Aunapuu and others 2008;Hoset and others 2017) and behaviors leading to devastation of non-preferred food plants (Jefferies and others 1994;Korslund and Steen 2006;Olofsson and others 2012) seem to constitute a context often overlooked in herbivore-plant-soil interactions.
Drawing from our findings, previous studies and the models of Pastor and Naiman (1992), Augustine and McNaughton (1998) and Bardgett and Wardle (2003), we illustrate two alternative pathways ( Figure 6) distinguishing the effects of selective herbivory (pathway 1) and decimation of nonpreferred food plants (pathway 2) on herbivoreplant-soil interactions. Low-intensity selective herbivory targeting subdominant, palatable forage drives decreasing process rates (Grellmann 2002;Stark and Grellmann 2002;Bardgett and Wardle 2003) predominantly via reducing litter quality and quantity (pathway 1). These retarding effects may offset the accelerating effects of feces and urine of resident herbivores on nutrient cycling (compare Bakker and others 2004). In contrast, whenever herbivory leads to reduction in slowly growing plants in the community in favor of faster growing palatable species, increasing litter quality and quantity act together with urine and excrement input toward increasing process rates (pathway 2).
However, spatial-temporal variation in rodent dynamics and in the environment governing plant growth (Brå then and others 2010; Hoset and others 2017) will constrain the ability of voles to drive increased process rates over entire tundra landscapes (compare Olofsson and others 2004b) after peak-year disturbance events. Shifts to a pathway 2 feedback are likely local, and contingent on the severity and frequency of disturbance (Aerts 2010), as well as resources facilitating vegetation change (Brå then and others 2010). Rapid increases in population density often lead to spillover of animals to un-productive secondary habitats (Hoset and others 2017), but this does not always lead to foraging of non-preferred plants others 2013, 2014). Finally, seasonal variation in vole diet (Soininen and others 2013) and limitation of vole population size through winter forage availability constrain decimation of dwarf shrubs, and especially V. myrtillus over large areas. Hence, vole populations may maintain complex mosaic landscapes with varying herbivore-plant-soil dynamics and successional stages. Our findings indicate that voles may support resilience of tundra landscapes under warming climate (Kaarlejä rvi and others 2015), provided that their population cycles continue to include high-amplitude peaks (compare Ims and others 2008; Korpela and others 2013).

CONCLUSIONS
In this paper, we have shown that high vole densities during peak years may increase ecosystem process rates in an inherently low-productive ecosystem by decreasing the abundance of slowgrowing plants of low palatability that would increase in abundance in the absence of voles. Based on our results and previous research, decimation of non-palatable plant species and selective foraging of palatable plants are both relevant features of herbivore-plant-soil interactions in low-productive tundra. These findings represent a novel insight into the ecosystem effects of population peaks of cyclic small rodents, which drive changes in plant communities across entire landscapes.