Decades of Recovery From Sheep Grazing Reveal No Effects on Plant Diversity Patterns Within Icelandic Tundra Landscapes

Tundra plant communities are often shaped by topography. Contrasting wind exposure, slopes of different inclination and landforms of different curvature affect habitat conditions and shape plant diversity patterns. The majority of tundra is also grazed by ungulates, which may alter topographically induced plant diversity patterns, but such effects may depend on the spatial scales of assessments. Here we ask whether topographically induced patterns of within (alpha) and between (beta) plant community diversity are different in contrasting grazing regimes. We studied plant communities within tundra landscapes that were located in the North and Northwest of Iceland. Half of the studied landscapes were grazed by sheep, whereas the other half was currently un-grazed and recovering for several decades (up to 60 years). Alpha and beta diversity were assessed on explicitly defined, nested spatial scales, which were determined by topographical units. Although we contrasted currently grazed vegetation to vegetation that witnessed several decades of grazing recovery, we found no statistically significant differences in plant diversity patterns. We relate these findings to the low resilience of our study system toward grazing disturbances, which has important implications for management practices in the tundra. Effects of topography on species richness were only found for specific spatial scales of analyses. Species rich topographical units were associated with relatively large biomass of plant growth forms that promote nutrient availability and potential plant productivity in the tundra, such as forbs. This suggests that biomass of such plant growth forms within habitats can be a useful proxy of potential plant productivity and may predict spatial patterns of plant species richness in tundra.


INTRODUCTION
Topographical structuring of different plant communities is particularly strong in alpine and arctic tundra ecosystems (Daubenmire, 1980;Evans et al., 1989;Ostendorf and Reynolds, 1998;Matsuura and Suzuki, 2013) and determines patterns of plant diversity (Jónsdóttir, 1984;Körner, 1995;Austrheim and Eriksson, 2001). The slope aspect for instance affects temperature and moisture regimes by either sheltering or exposing vegetation toward the main wind direction, with leeward slopes also accumulating larger amounts of snow than wind exposed slopes (Evans et al., 1989). Slope inclination further affects growing conditions with gentle slopes having higher snow accumulation, as well as higher amounts of water and nutrient influx as compared to steep sloping terrain (Ostendorf and Reynolds, 1998). The land surface curvature within such slopes can induce similar effects at smaller spatial scales, with concave depressions accumulating more snow, water and nutrients, as opposed to convex ridges (Florinsky and Kuryakova, 1996). Topography is therefore a main determinant of habitat conditions and potential plant productivity (further abbreviated as PPP) within tundra (Fisk et al., 1998). However, biotic interactions such as ungulate grazing may further modulate plant species diversity, but the role of grazing in shaping diversity in tundra ecosystems is still poorly understood.
Ungulate grazing influences species richness (Olff and Ritchie, 1998;Bakker et al., 2006;Kohyani et al., 2008;Bouahim et al., 2010) and relative abundance of plant species within communities (Augustine and McNaughton, 1998;Bråthen et al., 2007;Lezama et al., 2014) (further referred to as alpha diversity), but the direction of influence depends on intensity of grazing (Huston, 1979;Olff and Ritchie, 1998;Austrheim et al., 2008) as well as on the environmental growing conditions, which affect PPP (Proulx and Mazumder, 1998;Bakker et al., 2006;Lezama et al., 2014). Grazing reduces competitive exclusion of plant species under fertile growing conditions, promoting higher alpha diversity (Kaarlejärvi et al., 2017). Yet under nutrient poor conditions, even moderate grazing reduces alpha diversity (Proulx and Mazumder, 1998). Whilst grazing effects on alpha diversity are relatively well studied, there has been less focus on studying how free ranging ungulates modify the plant diversity component which resembles differences between plant communities within a landscape (further referred to as beta diversity).
Grazing can alter the relative abundance of graminoids, dicotyledonous forbs or woody plants by selectively grazing nutrient rich palatable species (Olofsson, 2006;Austrheim et al., 2008;Ravolainen et al., 2014), sometimes resulting in a more homogenous vegetation structure across habitats within a landscape (Bråthen et al., 2007;Lezama et al., 2014). Such homogenization effects across plant communities may lead to lower beta diversity within landscapes (Chaneton and Facelli, 1991;Olff and Ritchie, 1998;Ravolainen et al., 2010;Speed et al., 2013;Lezama et al., 2014). However, a reduction of beta diversity due to grazing is not always evident (e.g., Golodets et al., 2011), which is presumably due to dependencies on the spatial scales of assessments (Adler et al., 2001). Ultimately, patterns of both alpha and beta diversity result from the combined effects of apparent growing conditions and PPP, in relation to grazing intensity patterns within a landscape (Senft et al., 1987;Austrheim and Eriksson, 2001). In terms of tundra ecosystems, we still lack a profound understanding of how environmental conditions and grazing impact the diversity of plant communities, especially under the consideration of spatial scale dependencies. Grazing ungulates are present in almost all tundra areas throughout the northern hemisphere (Mulder, 1999;Van der Wal, 2006), being often managed as livestock or semi-domesticated herds (Barrio et al., 2016). However, we lack information on how grazers impact plant diversity in the tundra, and how topography may influence the impact at different scales (Ravolainen et al., 2010).
The aim of this study is to investigate grazing effects on alpha and beta diversity within tundra landscapes under explicit consideration of the spatial scales of assessment (Levin, 1992;Wiens, 1989;Barton et al., 2013). We used tundra landscapes in the North and Northwest of Iceland as our study system (Figure 1). Land use, including livestock grazing by sheep (Ovis aries L.), has strongly altered the natural vegetation of those landscapes since the island was settled around eleven hundred years ago (Lawson et al., 2007;Vickers et al., 2011;Brown et al., 2012;Eddudottir et al., 2016). Vegetation in tundra areas (at or above the Betula pubescens Ehrh. tree line) are nowadays strongly dominated by graminoids (grasses, sedges and rushes), but paleo records suggest that the same landscapes were previously dominated by deciduous shrubs (mainly Betula pubescens, Betula nana L. and thicket forming Salix species) and dicotyledonous forbs (Erlendsson et al., 2009;Streeter and Dugmore, 2014). Sheep grazing has maintained graminoid dominated vegetation and the relatively homogeneous appearance of most Icelandic landscapes today (Kristinsson, 1995;Thórhallsdóttir, 1996). Only a few small naturally protected areas escaped grazing in the Icelandic highland tundra. A study of one of these areas revealed that un-grazed vegetation had lower alpha diversity (species richness including bryophytes and lichens), whilst plant communities were stronger differentiated across topographic gradients compared to adjacent grazed areas (Jónsdóttir, 1984). However, the spatial coverage of this study is too small to generalize grazing impacts on plant diversity patterns of Icelandic tundra. Following Icelandic agricultural modernization during the 1940's, sheep grazing in many remote tundra valleys was terminated, creating the opportunity of vegetation recovery to un-grazed vegetation states. For the present study, we assessed plant diversity patterns of those valleys, which have been recovering from grazing for up to 60 years ( Figure 1A and Table 1) and compared them to similar valleys in close proximity that are still grazed. Within each valley, topography creates distinct growing conditions via contrasting slope aspects, elevations of differently inclined slopes and convex versus concave landforms, units that can be regarded as spatially nested (Figures 1B,C). We measured proxies of PPP in all these topographical units, using soil variables, such as soil pH, soil total carbon (C) and nitrogen (N) concentrations, as well as C:N ratios (Soil Survey Staff, 2011). High soil pH relates to a high proportion of bacteria in the soil microbial community, promoting nutrient availability and PPP in the tundra (Stark et al., 2012). Further, in Alaskan tundra, soil pH correlated positively with alpha diversity (Gough et al., 2000). We also assessed the biomass of plant growth forms (forbs, grasses, sedges and rushes, deciduous woody plants, evergreen woody plants, thicket forming shrubs) in each topographical unit. Such plant growth forms were previously shown to have different effects on PPP in the tundra, and to be preferentially grazed by ungulates   (Bråthen et al., 2007). Furthermore, these growth forms have different characteristics as niche constructors in tundra and may be important factors that drive plant diversity patterns . We expected different alpha diversity in recovering compared to grazed tundra valleys. In topographical units with high PPP (concave terrain, gentle slopes, slopes with leeward exposures), grazed valleys may have higher alpha diversity than recovering valleys, whilst alpha diversity may be lower in units with low PPP (convex terrain, steep slopes, slopes with exposures toward the general wind direction). We expected beta diversity to be generally lower in grazed compared to recovering valleys, due to homogenization effects. We also expected that proxies of PPP will confirm previous assumptions of differences amongst topographical units, with higher soil pH, total C and N concentrations, and lower C:N ratios in concave terrain, gentle slopes and slopes with leeward exposures, as compared their topographical counterparts. Topographical units of high PPP will be characterized by high biomass of productive plant growth forms such as forbs, grasses, sedges, rushes, or erect shrubs, whilst units of low PPP will have generally higher biomass of deciduous or evergreen dwarf shrubs (Bråthen et al., 2007). However, sheep grazing may reduce the biomass of plant growth forms in topographical units of high PPP (Bråthen et al., 2007).

Selection of Study Sites
We selected six valleys of contrasting sheep grazing regime in the Northwest and in the North of Iceland ( Figure 1A and Table 1). In general, three of the valleys were recovering from grazing (further referred to as recovering valleys), whilst the other three were presently still grazed (further referred to as grazed valleys). However, valleys had different land use history, which we found out by consulting the farmers who managed the respective valleys ( Table 1). The recovering valleys had recovering periods between 60 and 22 years. Grazed valleys were currently stocked by different numbers of sheep by the farmers, as indicated by the numbers of winter fed sheep ( Table 1). Sheep in Iceland are roaming freely throughout the whole summer season making it difficult to assess the grazing pressure within respective valley landscapes over time periods during the summer. We therefore conducted feces counts during our sampling (see Sampling collection below), providing us with an estimate of the different grazing intensities in the grazed valleys (Barrio et al., 2016). All valleys were of similar size, geographical orientation, and had comparable growing conditions ( Figure 1A), being situated north of 66 • N and within the low arctic subzone E according to the arctic bioclimatic zonation (CAVM Team, 2003). Long term (1949 to 2014) average monthly temperatures during the growing season (June to August) were 9.4 • C (min. 7.7 • C, max. 10.8 • C) in Northwest Iceland (weather station Bolungarvík), and 10.1 • C (min. 8.2 • C, max. 11.9 • C) in North Iceland (weather station Akureyri, 50 km south of the study sites) (Icelandic Meteorological Office 1 ). Average annual precipitation during the same measurement period was 841 mm (min. 590 mm, max. 1181 mm) in Northwest Iceland and 515 mm (min. 320 mm, max. 744 mm) in North Iceland (Icelandic Meteorological Office 2 ). Snow played an important role at our sites, with a continuous snow cover from October to Mid-June. All field sites were well outside the active volcanic zone in Iceland, with bedrocks of Tertiary basalts being more than 3.3 million years of age (Jóhannesson and Saemundsson, 2009). Therefore, in contrast with the active volcanic zone, the study sites were not heavily influenced by frequent deposits of volcanic ash or tephra, and typical soil types were classified as Brown Andosols (Arnalds, 2015). The valleys were glaciated during the last glacial maximum and became de-glaciated about 11 000 years ago (Norðodahl et al., 2008). The valley morphology was shaped by glacial erosion resulting in a typical U-shape. Steep valley slopes were covered with scree and solifluction lobes, which were important in shaping the smaller-scale topography.
The vegetation was generally classified as "low shrub tundra" (CAVM Team, 2003;Walker et al., 2005). The prevailing wind direction from the East and North East (Einarsson, 1976) lead 1 http://en.vedur.is/Medaltalstoflur-txt/Manadargildi.html 2 http://en.vedur.is/Medaltalstoflur-txt/Arsgildi.html to greater snow deposition on west facing than on east facing slopes (Evans et al., 1989). The concave morphology of valley slopes caused a vertical topo-sequence from xeric to mesic and moist conditions toward the valley bottom ( Figure 1B). However, small streams and alluvial fans, which ran down the valley slopes, caused a horizontally altering pattern of convex and concave landforms within this vertical sequence, leading to differences in growing conditions at even smaller scale. Domesticated sheep were the main vertebrate herbivores in our valleys. Other herbivores could occur occasionally such as flocks of migratory geese, resident ptarmigans (Lagopus muta Montin, own observations of droppings) or wood mice (Apodemus sylvaticus L.) (Unnsteinsdóttir and Hersteinsson, 2009).

Sampling Design
We aimed for a sampling design that enabled us to capture the vegetation differentiation according to the three spatially nested topographical units, i.e., topographical differentiation according to (i) the slope aspect (largest spatial scale), (ii) high and low elevations within slopes (intermediate spatial scale) and (iii) concave and convex landforms within different elevations (smallest spatial scale).
There were no available vegetation maps for our valleys and available digital data were too coarse to allow stratification by small scale landform differences. We therefore emphasized that all steps of the sampling design were as transparent as possible and based on clearly defined criteria (Mörsdorf et al., 2015). We used topographical maps within a geographical information system (esri ArcGIS 10.1), and marked a cross section along the rivers that run through the bottom of the valleys ( Figure 1B). To ensure a spread of sampling units throughout the valley, we further stratified the sampling to three zones which were defined by specific distances from the coast line: zone A (1-2 km from the coast), B (2-3 km from the coast) and C (3-4 km from the coast) ( Figure 1B). Within each zone and for both valley slopes, we defined transects which were running from the river at the valley bottom upwards the valley slopes. These transects were spaced at 100 m intervals along the river line and had to traverse a concave valley slope. Transects that crossed transitions to convex topography were discarded. We also used aerial photographs to discard transects that crossed boulder fields, most of which had very low vascular plant cover. To restrict sampling to the foot of the slopes (mild snow bed conditions) and the more inclined parts of the slope (mesic conditions), we noted the GPS coordinates of all remaining transects that intersected with a contour line of 40, 60, and 80 m elevation for zone A, B and C, respectively. The difference in elevations for each zone was due to a general uplift of the valley bottom from the coastline to the inner parts of the valleys. Those GPS coordinates built the sampling frame for the present study.
Two GPS coordinates were selected randomly from the sampling frame of each zone, one from either side of the valley. In the field we visited these coordinates and used a priori defined rules to guide us to sampling units of interest that are shaped on smaller spatial scales, i.e., convex and concave landforms: Arriving at the GPS location, we kept the altitude and moved horizontally toward the sea until we reached the approximate edge of an alluvial fan. Alluvial fans themselves were characterized by a convex surface topography, and passing the edge lead into concave terrain. The transition between both terrain forms was smooth in reality though, and we had to visually assess the approximate location of the border between both landforms, which entailed some subjective judgments (Supplementary Figure 1A). However, this subjective definition of the border between both landforms was only used to determine the location where we placed the center of a 30 m long measuring tape. To avoid a subjective placement of vegetation plots, whilst assuring that we sampled within convex and concave terrain, we stretched both ends of the tape into the convex and concave landform respectively (Supplementary Figure 1B). After that, we sampled the vegetation systematically, starting at both ends of the measuring tape (see section "Sampling Collection"). We repeated the same procedure at an elevation 60 m above the selected GPS coordinates to sample vegetation data at steeper (mesic) parts of the valley slopes.

Sampling Collection
The vegetation was analyzed across the concave and convex landforms along the 30 m measuring tape, using 40 × 40 cm plots. Beginning at both ends of the measuring tape, we sampled four plots within each landform at spatial intervals of three meters. To measure plant species abundance, we applied a refined version of the point intercept method (Jonasson, 1988), which is designed to sample vegetation over extensive spatial scales (Bråthen and Hagberg, 2004). We used a 40 × 40 cm metal frame with 5 metal pins of 2 mm diameter, one in each corner of the frame and one in its center. The frame was placed at the uphill side of the measuring tape and at each pin all hits through the vascular plant canopy were recorded and vascular plants were identified down to species level. To measure vascular plant species richness, we recorded all additional species within the plot, which were not hit by the pins. We analyzed 576 plots in total.
Soil samples were taken next to each vegetation plot: Approximately 50 g of fresh soil were excavated from the soil surface down to about five cm soil depth, which corresponded to the rooting zone at our study sites. The four soil samples of each landform (convex versus concave) were pooled and stored in cooled conditions until arrival in the lab (max. four days). Vegetation and soil sampling was conducted during the peak of growing season (11 Jul 2012 to 15 Aug 2012). In the lab, soil samples were air dried at ambient temperature, sieved using a two millimeter mesh size and homogenized with a mortar. We measured the soil pH after extraction in distilled water, using a soil to water ratio of 1:5 (method adapted for dried soil samples from Blakemore et al., 1987). In addition, we analyzed total C and total N concentration of all samples using a vario MAX cube CN analyzer 3 .
To assist estimates of current grazing pressure, we counted the number of herbivore droppings within a one-meter zone along the 30 m measuring tape ( Table 1).

Selection of Diversity Metrics and Plant Growth Form Classification
We selected alpha and beta diversity metrics that incorporated occurrence as well as the abundance information of plant species within and between communities. We used species richness to measure properties of alpha diversity in terms of species occurrences and Gini-Simpson index in terms of relative species abundances within a community ( Table 2). For beta diversity, we used dissimilarity indices that excluded information on joint species absences. We chose Jaccard dissimilarity to reflect community differentiation based on species occurrences. Representing community differentiation based on relative species abundances, we used a modified version of Gower's distance (Anderson et al., 2006). This "Modified Gower" distance (sensu Anderson et al., 2006) enabled us to weigh the change in abundance over orders of magnitude. By applying a prior logarithmic transformation on the raw abundances, where weighing is done according to the base of the logarithm (Anderson et al., 2006), the distance measure can be interpreted as an average change in orders of magnitude per species between two different plant communities. We chose to use a log base of two for this study (and further termed the distance measure "MG2" throughout this article), as this gives most weight to a change in relative species abundance. Using a log base of two weights a doubling in abundance of one species as much as a compositional change of one species. We used the R environment for all our data evaluations (version 3.1.3 4 ) and applied vegdist and decostand function of the vegan package to calculate Jaccard dissimilarities and MG2 distances. beta Modified Gower distance, using a log base of 2 = sum (w k (abs(x' 1kx' 2k )/sum (w k ) x': log 2 (x) + 1 x 1k : abundance of species k in community 1 x 2k : abundance of species k in community 2 w k = 0 when x 1k = x 2k = 0, otherwise w k = 1 Frontiers in Ecology and Evolution | www.frontiersin.org To evaluated differences of plant growth form biomass between topographical units and grazing regimes, we chose growth forms that were previously found to affect PPP, whilst being selectively grazed by ungulates (Table 3). This classification was according to findings of Bråthen et al. (2007). The plant growth form classification we used can also be relevant to predict plant diversity in the tundra . As another proxy of PPP, we additionally calculated the total vascular plant biomass (sum of all growth forms) for each landform.

Statistical Analyses
We analyzed our data by applying linear mixed effects models, using the nlme package in R. For analyses of species diversity, we regarded the spatial nestedness of topographical units and the data recordings within those units as different grain sizes ( Figure 1C). For analyses of the smallest grain size we aggregated all plant hits (or species numbers for species richness) of the four plots within each concave and convex landform. Accordingly, we aggregated all the plant data within each high and low elevation transect, representing an intermediate grain size. Finally, we aggregated all the plant data within east and west facing slopes within one zone, which was the largest grain size in our study.
Next, we converted all plant hits into biomass using weighted linear regression methods (Bråthen and Hagberg, 2004), which yielded the average biomass of each plant species per square meter (grams * m −2 ) for each grain size of assessment. The conversion was based on Ravolainen et al. (2010). Species in our Icelandic data that did not exist in their study were given the conversion factor of the morphologically most similar species (Supplementary Table 1).
Alpha diversity was assessed by setting species richness or Gini-Simpson index as response variable in the linear mixed effects models. We defined the effects of grazing regime, the topographical unit and their two-way interaction as fixed factors in our models. As none of the interactions were statistically significant (based on a 5% significance level), we reduced all models to include the grazing regime and topographical unit as additive fixed effects. Models thereby included the grazing regime (recovering from grazing versus grazed) and landform (convex versus concave), the grazing regime and elevation (high versus low) or the grazing regime and slope aspect (east versus west facing) as additive fixed effects, depending on the spatial grain size of the analyses (Supplementary Table 2). Beta diversity was calculated as the dissimilarity (Jaccard, MG2) between FIGURE 2 | Model estimates of diversity using a small spatial grain size. (A,B) Alpha diversity is presented given the influence of grazing regime (grazed versus recovering) and landform curvature (concave versus convex). (C,D) Beta diversity, which was calculated between the landform curvature units, is presented given the influence of the grazing regime. Predicted means of alpha and beta diversity are based on linear mixed effects models. " + " indicates marginally significant effects and error bars represent 95% confidence intervals of predicted means.
topographical units for the respective grain sizes of assessment ( Figure 1C). Models for beta diversity had therefore either Jaccard dissimilarity or MG2 distance as response variables and the grazing regime as a fixed factor. The random structure of all our models reflected the spatial hierarchy of our design on the respective spatial scale (Supplementary Table 2). We used the exact same modeling approach as for alpha diversity, to assess grazing regime and topography effects on soil variables, plant growth form biomass and total vascular plant biomass. As none of the interactions were statistically significant, we repeatedly reduced the models to only include the additive effects of grazing regime and topography (Supplementary Table 3). This procedure was followed for small and medium scale analyses only, as those were the only scales of analyses where we found effects on vascular plant species diversity.
We assessed the assumptions of all models in terms of constant residual variance and normality of residuals, and checked for outliers, using diagnostic plots. The response variables of all plant growth form models, and of the total vascular plant biomass model, were log e (x + v) transformed to fulfill model assumptions, with v being the smallest biomass value of the data set. Within the results section, we report statistically significant effect sizes based on a 5% significance level. Based on a 10% significance level, we annotate effects as "marginal."

Alpha and Beta Diversity Patterns of Contrasting Grazing Regimes and Topography
The species richness within communities was not different between grazing regimes, but the landform curvature displayed effects on species richness, with marginally lower values in convex The intercept represents average diversity within currently grazed landscapes and within concave landforms, high elevation or east facing slopes, depending on the spatial grain size of assessments.
than in concave landforms (Figure 2A and Table 4). For Gini-Simpson index, we found no difference between grazing regimes or landform units ( Figure 2B and Table 4). We found no differences in beta diversity between grazing regimes on small spatial scale (Table 4). Both, Jaccard dissimilarities ( Figure 2C) and MG2 distances ( Figure 2D) were similar in recovering and in grazed valleys. Analyses on intermediate spatial scale revealed no differences in species richness between grazing regimes, but species richness was lower at low compared to high elevations ( Figure 3A and Table 4). Using the Gini-Simpson index as a measure of alpha diversity, we found no difference between grazing regimes and elevations ( Figure 3B and Table 4). Also, beta diversity between high and low elevations was not different between grazing regimes (Table 4), with both Jaccard dissimilarities ( Figure 3C) and MG2 distances ( Figure 3D) being similar in recovering and grazed valleys.
Using the large spatial scale for analyses did not reveal any difference in species richness ( Figure 4A and Table 4) or Gini-Simpson index (Figure 4B and Table 4) between grazing regimes or contrasting slope aspects. For measurements of beta diversity between east and west facing slopes, we found no differences between grazing regimes (Table 4). Both Jaccard dissimilarities ( Figure 4C) and MG2 distances ( Figure 4D) were similar in recovering and grazed valleys.

Soil Variables and Plant Growth Form Biomass of Contrasting Grazing Regimes and Topography
Some soil variables showed differences between contrasting landform units but there was no difference between grazing regimes at this spatial scale ( Table 5). We found significantly lower soil pH and higher soil carbon and nitrogen concentrations within concave than within convex landforms. For intermediate spatial scale, we also found elevation differences in soil variables, but no differences between grazing regimes ( Table 6). Low elevations had lower values of soil pH and higher values of soil carbon and nitrogen concentrations than high elevations.
Deciduous and evergreen woody plants were the dominating growth forms within our study sites, but we found no significant biomass differences for any plant growth form between grazing FIGURE 3 | Model estimates of diversity using an intermediate spatial grain size. (A,B) Alpha diversity is presented given the influence of grazing regime and elevation. (C,D) Beta diversity, which was calculated between different elevations, is presented given the influence of the grazing regime. Predicted means of alpha and beta diversity are based on linear mixed effects models. "*" indicates statistically significant effects and error bars represent 95% confidence intervals of predicted means. Table 4), based on a small spatial grain size. However, the landform curvature had a significant effect on the biomass of forbs, which were more abundant within concave than convex landforms (Figure 5A and Supplementary Table 4). Table 4). Also the elevation had strong but contrasting effects on plant growth forms. We found larger biomass of deciduous woody plants and forbs at high compared to low elevations, whereas graminoid plants had larger biomass at low elevations ( Figure 5B and Supplementary Table 5). Analyses of plant growth form biomass on intermediate spatial scale also revealed marginally larger biomass of thicket forming Betula and Salix species in recovering compared to grazed valleys ( Figure 5C and Supplementary Table 5). Total vascular plant biomass was not different between elevations and grazing regimes on intermediate spatial scale (Supplementary Table 5).

DISCUSSION
We originally expected the grazing regime to affect alpha diversity, with potentially opposite effects in topographical units of high versus low PPP (Proulx and Mazumder, 1998;Bakker et al., 2006). However, for the three spatial scales we investigated, we did not find such interaction effects, or any effects of the grazing regime on alpha diversity. On all spatial scales of assessment, we did not find evidence that beta diversity differed between grazing regimes. Topography alone was the structuring element of plant diversity patterns at our sites, but these patterns depended on spatial scale. In line with these findings, we found almost no effects of the grazing regime on soil variables or on plant growth form biomass. Some soil variables and plant growth form biomasses were significantly different between topographical units.

Grazing History and the Topography as Drivers of Current Plant Diversity Patterns
We relate the lack of grazing effects on alpha and beta diversity to the low resilience of formerly grazed tundra vegetation after the cessation of grazing. It has been shown that ungulates can push vegetation into different stable states (Westoby et al., 1989;Laycock, 1991) and the same has been suggested for the arctic tundra ( Van der Wal, 2006). In general, two important characteristics of an ecosystem are relevant to evaluate its resilience to grazing impacts, which are i) the history of grazing and ii) the availability of resources (Milchunas et al., 1988;Cingolani et al., 2005). In Iceland, livestock grazing has been extensively practiced since the time of the Norse settlement 1100 years ago (Erlendsson et al., 2009). Sheep grazing is assumed to have maintained graminoid dominated vegetation in many Icelandic landscapes, preventing vegetation shifts back to states which were dominated by shrub and forb species (Kristinsson, 1995;Thórhallsdóttir, 1996). Grass dominance can also be expected because long grazing history usually selects for a subset of the plant species pool which is tolerant to grazing, preventing re-establishment of species that are less tolerant (Milchunas et al., 1988;Cingolani et al., 2005). It is important to notice that the vegetation in our valleys was dominated by woody plants (Supplementary Figure 2) even though their biomass is usually very low under high grazing pressure (Olofsson, 2006;Austrheim et al., 2008). The seeming discrepancy to tundra studies that showed graminoid dominated vegetation under grazing (Jónsdóttir, 1984;Eskelinen and Oksanen, 2006;Olofsson, 2006) is presumably because graminoids profit only under relatively high animal densities, including high defecation rates and trampling. However, we also have to consider the specific grazing histories of each valley that was grouped into the categories "grazed" and "recovering" in our analyses. Relatively little woody plant biomass was found in two currently grazed valleys (Skálavík and Ingjaldssandur). Personal communication with farmers informed that these valleys were intensively grazed in the past, leading to a persisting low biomass of woody growth forms. The currently most intensively grazed valley (Thorgeirsfjördur) had relatively high woody plant biomass (Supplementary Figure 2), which may indicate lower historical grazing intensity. On the other hand, the biomass of evergreen shrubs was highest in Thorgeirsfjörður, which may reflect the responses of these unplatable plants to the high current grazing pressure (Bråthen et al., 2007). Furthermore, we worked in an ecosystem of relatively scarce nutrient availability, which may cause recovery to un-grazed vegetation states to take a long time (Cingolani et al., 2005). Also vegetation analyses from tundra sites in the Icelandic highlands indicated that recovery from continuous grazing is slow: No difference in vegetation properties inside and outside an exclosure was detected after 4 years without grazing (Jónsdóttir et al., 2005), but grazing influenced vegetation could theoretically persist for decades (Laycock, 1991) or even centuries after grazing cessation (Ransijn et al., 2015). On average, we only found weak trends of higher woody plant biomass, represented by Salix and Betula shrubs, in recovering valleys of our study ( Figure 5C). One general explanation may lay in the dynamics of recruitment pulses of such shrubs. Those are generally temperature dependent, but operate with a time-lag even after years with good growing conditions in the tundra (Büntgen et al., 2015). Climate warming during the last decades has been shown to promote shrub growth throughout the Arctic (Myers-Smith et al., 2011;Elmendorf et al., 2012) and future increases of temperature may speed up the currently weak trends of shrub biomass increase that we found in recovering valleys here. Furthermore, higher snow fall during the winter time could promote shrub growth via enhancing nutrient availability in the immediate surrounding of shrubs (Sturm et al., 2005). Herbivores have been shown to counter-act shrub expansion (Christie et al., 2015;Bråthen et al., 2017). Yet, it is difficult to predict how interactions of climate driven shrub expansion and herbivory will affect plant diversity in future. Some effects of herbivory may indeed decrease plant diversity in the tundra, since shrub canopies can be important niche constructors and promote plant species richness in tundra (Bråthen and Lortie, 2015;Bråthen and Ravolainen, 2015). However, studies on the interactive effects of herbivory and and grazing regimes. "*" indicates statistically significant effects based on a 5% significance level. Error bars represent 95% confidence intervals of the predicted means. Note that y-axes have different dimensions for specific growth forms due to large differences in biomass.
shrub expansion on plant diversity require specifically designed experiments and coordinated long-term monitoring (Christie et al., 2015;Barrio et al., 2016).
Across grazing regimes, we found that current patterns of plant diversity in our valleys relate to topography at various spatial scales. Small scale contrasts in landform curvature and intermediate-scale contrasts in elevation were associated with differences in plant species richness, whilst contrasts of slope aspect at a large scale had no effect. Species richness within communities (alpha diversity) has generally been attributed to the conditions of PPP (e.g., Grime, 1973;Tilman, 1987). Contrasts in such conditions may have been more pronounced amongst landforms and elevations, than amongst slope aspects at our sites. However, our data exemplifies how the utilized surrogate measure of PPP affects interpretations of productivity-diversity patterns across scales.

Spatial Patterns of Plant Diversity and Their Relation to Soil Variables and Plant Growth Form Biomass
Our analyses of soil variables and plant growth form biomass were in accordance with the species diversity patterns we found. There were no effects of grazing recovery (except for a marginal increase of Betula and Salix thickets as outlined above), but soil variables and plant growth form biomass differed amongst landforms and elevations. The theoretical basis of how PPP affects plant diversity is often difficult to test in practice, since field studies have to rely on a variety of proxies for PPP, as we did here (Whittaker and Heegaard, 2003;Adler et al., 2011). The utilized proxy can thereby affect the outcome of study results. In our case, soil carbon and nitrogen concentrations indicated higher PPP in concave than in convex landforms. Therefore, species richness patterns between landform curvatures may indicate a positive relationship between PPP and species richness. However, our study also shows that productivity -diversity relationships could be interpreted in the opposite way when using the intermediate spatial scale of assessment (low and high elevation). Here we found higher richness at high, compared to low, elevations. The high units were being characterized by lower soil carbon and nitrogen concentrations.
A more congruent interpretation of productivity -diversity patterns can be based on the biomass differences of plant growth forms within the communities, which were shown to be as important as abiotic growing conditions in determining alpha diversity in the tundra . Across scales, the topographical units with highest species richness at our sites (concave landforms; high elevations) all had a larger forb biomass than their less diverse counterparts (convex landforms; low elevations). Forb dominated communities, have often high bacterial:fungal ratios (Sundqvist et al., 2011), causing fast rates of nutrient recycling (Eskelinen et al., 2009). Plant communities with relatively large biomass of forbs may thus, at times, promote higher inorganic nutrient availability than communities that are dominated by other growth forms, as for instance evergreen shrubs . In other words, forbs may have promoted species richness at our study sites. On intermediate scale, we also found that graminoids had more biomass at low than at high elevations. Although very large graminoid biomass may contribute to the exclusion of plant species due to both nutrient competition and shading, such negative plant interactions are only plausible in communities with biomass of more than 400 g * m −2 (Bråthen and Lortie, 2015). We did not encounter such saturated communities in our study ( Figure 5, and total biomass estimates of Supplementary Tables 4, 5) and therefore assume that facilitation through the abundance of forbs is the major mechanism that promoted plant species richness at our sites. Besides the promoting effects on nutrient cycling, forbs may also enhance plant species richness via being the plant group which had by far the largest species pool size at our sites (Table 3).
Altogether, productivity -diversity patterns in our study support the notion that productivity -diversity relationships depend on the spatial grain size of investigation (Mittelbach et al., 2001;Chase and Leibold, 2002;Whittaker and Heegaard, 2003). More specifically, we conclude that spatial patterns of vascular plant species richness at our sites are best explained by the biomass of different plant growth forms that dominate in specific topographical units. In tundra, plant growth form-specific biomass can be a useful proxy of PPP  and should perhaps be preferred to abiotic variables, such as soil nutrient concentrations. Topographical units that promote biomass of growth forms with a large actual species pool size (Zobel, 1997), such as forbs here, may be especially species rich in un-saturated tundra plant communities.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
IJ conceived of the research idea. MM, VR, IJ, TT, and NY designed the research. MM and IJ collected the data. MM, VR, and NY performed statistical analyses. MM, VR, and IJ with contributions of TT and NY wrote the manuscript. All authors read and approved the submitted version of this manuscript.

FUNDING
This study was funded by the Icelandic Research Fund (grant # 110235021-3) and the University of Iceland Research Fund (grants to IJ). The article processing charge was funded by the University of Freiburg via the funding program "Open Access Publishing."