Reindeer grazing controls willows but has only minor effects on plant communities in Fennoscandian oroarctic mires

ABSTRACT

Shrubification of arctic tundra is a well-recognized phenomenon, and it can be particularly rapid in moist habitats. Reindeer grazing can inhibit shrubification, but grazing impacts on mire vegetation have been overlooked. We studied grazing effects on plant communities and Salix lapponum in oroarctic mires at the border of Finland and Norway. We compared plant community structure and S. lapponum abundance and traits between (1) grazed fens (Finland); (2) experimental exclosures (Finland), where reindeer have been kept out for 13 years; and (3) nongrazed fens (Norway). Grazing effect on shrubification was assessed using the Normalized Difference Vegetation Index (NDVI) and leaf area index (LAI). We did not find a uniform direction of vegetation change connected to the exclosure treatment, and grazing treatments were overlapping in multivariate ordination. Neither NDVI nor LAI indicated clear differences. Instead, significant results were revealed in total abundance of species groups and in S. lapponum traits. The cover of bryophytes was significantly lower under free grazing. Reindeer grazing reduced the abundance, height, and flowering and increased leaf N concentration of S. lapponum. We conclude that reindeer grazing controls willows and affects total abundance of important species groups, and plant community structure is resistant to grazing effects in oroarctic mires.

KEYWORDS:

• Herbivory
• NDVI
• plant communities
• Salix
• shrubification

Introduction

A major part of the circumpolar arctic tundra and northern boreal regions is grazed and trampled by either semi-domestic or wild reindeer (Rangifer tarandus) populations. Reindeer are known to modify tundra vegetation and soil processes, but ecosystem responses to grazing and trampling vary between habitats and vegetation types (see Bernes et al. [2015] and Christie et al. [2015] for recent reviews). In tundra heaths, reindeer grazing can increase nutrient cycling (Olofsson et al. 2001; Stark, Strömmer, and Tuomi 2002) and cause transitions between alternative vegetation states (van der Wal 2006). Lichen cover and biomass tend to decrease in grazed areas (van der Wal, Brooker, et al. 2001; den Herder, Kytöviita, and Niemelä 2003) because they are preferred winter forage for reindeer (Nieminen and Heiskari 1989) and are sensitive to trampling during snow-free seasons (Forbes 2010; J. Kumpula et al. 2014). In dry, nutrient-poor habitats, grazing results in transitions from a lichen-rich to a moss-dominated vegetation state (Manseau, Huot, and Crête 1996; van der Wal 2006). In more productive dry tundra habitats, grazing tends to increase primary production and often leads to a shift from moss- or dwarf shrub– to graminoid-dominated vegetation (Olofsson et al. 2001; Olofsson, Stark, and Oksanen 2004; Post and Pedersen 2008). The effects of reindeer grazing in wet tundra habitats are less studied. Though willows and sedges in peat-accumulating mires provide important forage in summer, little is known about grazing effects on mire vegetation.

Though reindeer grazing is an important form of land use, climate change is rapidly proceeding in Arctic regions (Overland et al. 2017), and the role of herbivory may change. Remote sensing observations have indicated “greening” (measured by the Normalized Difference Vegetation Index, NDVI) in central and eastern Siberia, Canadian tundra, and the north slope of Alaska since the early 1980s (Jia, Epstein, and Walker 2009; Epstein et al. 2012; Raynolds et al. 2013). It is important to note, however, that the greening pattern has not been uniform throughout the arctic tundra zone (Phoenix and Bjerke 2016). Observed greening links to increased growth, abundance, and range of deciduous shrubs, particularly tundra willows (Salix spp.; Forbes, Macias-Fauria, and Zetterberg 2010; Blok et al. 2011; Myers-Smith et al. 2011; Macias-Fauria et al. 2012). Alpine willows have expanded to higher elevations likely in response to warming (Myers‐Smith and Hik 2017), and recent evidence suggests climate change–induced expansion of evergreen ericoid dwarf shrubs in Scandinavian mountain regions (Vuorinen et al. 2017; Vowles, Gunnarsson, et al. 2017). Shrub expansion can result in several positive feedbacks to climate change, such as decreases in surface reflectance (albedo; Blok et al. 2011) and increases in winter soil temperatures (Myers‐Smith and Hik 2013), both of which can further enhance shrub growth. However, reindeer grazing can inhibit climate change–induced shrub expansion (Olofsson et al. 2009; Ravolainen et al. 2014) because deciduous shrubs are a preferred food source for reindeer during summer (Nieminen and Heiskari 1989; Klein 1990).

As recently noted by Martin et al. (2017) and Vowles, Lovehave, et al. (2017), research on shrub expansion and reindeer grazing impacts is limited geographically and in terms of vegetation types. Northern mires (bogs and fens) comprise large carbon reservoirs and play a major role in regulating global carbon cycle and climate (Yu 2012). The capacity of ungulate herbivores to alter carbon dynamics in arctic mires has recently been recognized (Falk et al. 2015; Stark and Ylänne 2015). Herbivores can mediate ecosystem responses to warming and alter carbon dynamics in arctic tundra and northern mires through their effects on vegetation and soil processes (Ward et al. 2007; Cahoon et al. 2012; Laiho, Penttilä, and Fritze 2017). Given that northern mires represents approximately one third (500 ± 100 gigatons) of estimated global soil carbon reservoirs (Yu 2012), grazing-induced changes in arctic and oroarctic mire vegetation could have significant implications for the global carbon budget (Stark and Ylänne 2015). Although mires are important feeding habitats for many herbivore species, studies of grazing effects on mire vegetation and particularly on carbon budget are still rare.

Here, we present results from a long-term (13-year) experiment on the impact of reindeer grazing on mire vegetation. The study sites locate in sedge- and willow-dominated fens across the border of Finland and Norway. The border has been fenced since the late 1950s to prevent reindeer entering winter pastures in Norway during summer. We compared plant community structure between (1) grazed fens (Finland); (2) experimental exclosures (Finland), where reindeer have been kept out for 13 years; and (3) nongrazed fens (Norway). In addition, we compare several growth characteristics and flowering of Salix lapponum, the main willow species in the study area and a key reindeer forage species in summer, according to local reindeer herders (Juha Magga, pers. comm.). Other Salix species (mainly S. glauca, S. phylicifolia) occur sparsely, whereas S. lapponum is found on all mire sites in the area and forming continuous thickets in most favorable sites. Results of an earlier short-term (4-year) study of the same experiment showed no significant effects on vegetation cover or on S. lapponum height (Kitti, Forbes, and Oksanen 2009). We hypothesize that after a longer period (13 years) of reindeer exclusion, the grazing impacts have dissipated and vegetation has already become more similar to that in nongrazed fens on the Norwegian side. We further expect that reindeer exclusion will have enhanced the growth and reproductive success of willows and that increased abundance and size of willows will be reflected as higher NDVI values. Alternatively, the legacy of reindeer grazing pressure could still be present in the mire vegetation, either as continued suppression of willow growth or in vegetation composition.

Materials and methods

The study area

The Jauristunturit study area is located in the Fennoscandian oroarctic tundra zone (Virtanen et al. 2016). The experiment was carried out at ten study sites near the border of Finland and Norway (68°49ʹ N, 23°49ʹ E, 450–510 m.a.s.l.) in a non-permafrost region (Figure 1). The mean annual air temperature in the closest (27 km southwest of the study area) meteorological observation station (Enontekiö, Näkkälä, 68°36ʹ N, 23°35ʹ E, 370 m.a.s.l.) was −1.3°C and annual precipitation was 456 mm during the years 2001–2012 (Finnish Meteorological Institute, open data). The mean temperature of the warmest month (July) was +13.4 °C and that of the coldest month (February) was −14.6 °C. Snow cover persists on average for 152 days. The climate is slightly colder than indicated by the meteorological station data, due to the ca. 100-m higher elevation. The landscape is a mosaic of Betula nana–Cladina-type heaths and fjells (Virtanen, Oksanen, and Razzhivin 1999) and peat-accumulating fens. One characteristic feature of the studied fens is the abundance of S. lapponum, a willow species subject to summer grazing by reindeer. Other dominant species of the fens are sedges, such as Carex rostrata, Carex rotundata, and Eriophorum angustifolium. Among bryophytes, Sphagnum riparium, Sphagnum lindbergii, Sphagnum teres, Paludella squarrosa, and Warnstorfia sarmentosa are characteristic.

In the mid-1950s, a 3-m-high reindeer fence was built along the Finnish–Norwegian border. In the study area, the Finnish side of the fence belongs to the Näkkälä reindeer herding district, and it is grazed mainly in summer. In early summer, reindeer migrate to higher hills and mires and feed on new shoots of willows and graminoids. Reindeer feeding on willows is selective browsing of young leaves, shoots, and flowers, whereas graminoids are typically grazed to the ground. In this article, we refer to both feeding patterns. The grazing pressure in summer is approximately 10–12 reindeer per square kilometer over the 150-km² district (Kitti, Forbes, and Oksanen 2009). From the late 1970s to the mid-1990s, the Finnish side was used also as a winter pasture (T. Kumpula 2006).

The Norwegian side is a traditional winter range for migratory herds of reindeer, and summer grazing (by Finnish reindeer) has been ultimately prohibited by the fence building in the late 1950s. In winter, reindeer feed mainly on lichens in dry habitats (Nieminen and Heiskari 1989), if they are available (Helle and Jaakkola 2008; J. Kumpula et al. 2014). Therefore, the studied fens in Norway are virtually nongrazed year-round. It is important to bear in mind, however, that the nongrazed Norwegian plots are not similarly fenced and some occasional minor grazing may have taken place. The reindeer fence along the national border prevents reindeer entering winter pasture during summer, and it is an experimental exclusion of reindeer summer grazing. This was a true manipulation of the grazing regime that completely divided reindeer culture in two parts of a wilderness area. The reindeer fence along the national border has been used as an experimental design in many studies of grazing impacts (e.g., Käykhö and Pellikka 1994; T. Kumpula 2006; Cohen et al. 2013), and it is widely recognized as uniquely valuable for long-term study of grazing effects. Presumably, the starting point in the 1950s was similar in terms of summer grazing in all studied fens, because of their close proximity (Figure 1). However, at least since the 1950s, reindeer have been kept out of the Norwegian side in summer.

Figure 1. Locations of study sites across the national border with reindeer fence between Finland and Norway. Reindeer grazing pressure is lower and missing during summer on the Norwegian side, where abundance of lichens on dry tundra gives a whitish tone in the false-color image. The Finnish side has a historical legacy of high grazing pressure throughout the year, indicated by the lack of lichen cover, and is now used as summer pasture imposing high grazing pressure on mire habitats.

Figure 2. Grazing treatments have overlapping ranges with the nongrazed Norwegian subplots that have less variation in a two-dimensional nonmetric multidimensional scaling ordination. (a) Data from exclosure subplots in 2002 (before the effect of exclosure) and grazed subplots in 2015 combined as one group (Finland, grazed). Finland: Exclosures consist of subplots in exclosures in 2015. Norway: Nongrazed area consists of subplots in Norwegian fens in 2002 and 2015. (b) Combined data from 2002 and 2015 show nearly complete overlap; the range of vegetation variation is similar between years. Combinations of different years and grazing regimes are shown in more detail in Appendix 4.

Experimental design and fieldwork in 2002 and 2006

An exclosure experiment was established in 2002. Kitti, Forbes, and Oksanen (2009) described the details of experimental design and only the main aspects are summarized here. During summer 2001 and winter 2002, ten separate fens were selected as study sites with local reindeer herders, five sites on both sides of the border fence. At each site in Finland, three exclosures (5 m × 5 m) were built, and three freely grazed plots were established and permanently marked with plastic pipes. Exclosures effectively prevent reindeer grazing and trampling, though smaller vertebrate herbivores, such as arctic hare (Lepus timidus), willow ptarmigan (Lagopus lagopus), and rodents, can freely enter. In Norway, three plots per site were established and marked permanently with plastic pipes. Thus, three grazing treatment levels are included: (1) free reindeer grazing and trampling that mainly takes place in summer (FI), (2) exclusion of reindeer grazing since 2002 (FI), and (3) no reindeer grazing since the 1950s (NO). In July 2002, abundance of all plant species was estimated using the point intercept method (100 points per 0.25-m² subplot) and S. lapponum height measurements were made in exclosures (FI) and in nongrazed plots (NO). The height of S. lapponum was recorded by measuring the height of uppermost shoots (maximum of five measurements per subplot). In 2006, only vascular plant species were estimated with the point intercept method and S. lapponum heights were measured again in all subplots, including freely grazed plots (FI). One site from both Finland and Norway was excluded from analysis by Kitti, Forbes, and Oksanen (2009) due to disparate vegetation type.

Repetition of vegetation sampling in 2015

We repeated the vegetation measurements described above in late July 2015. We did not manage to relocate two of the original Norwegian sites, and at most remaining sites it was not possible to exactly relocate each original subplot (0.25 m²). In such cases, we randomly selected three subplots from each 5 m × 5 m plot (n = 39). In all cases, we ensured that the same rules of subplot selection were employed, in terms of distance from the edge of the 5 m × 5 m plot and between subplots. The cover of each species was estimated in all subplots (n = 117) using the point intercept method (100 points per 0.25-m² subplot). All mosses and vascular plants were identified to species level. Hepatics in the order Jungermanniales were treated collectively as one unit, but Aneura pinguis and Marchantia polymorpha were identified to the species level. For Salix spp. and B. nana, we counted each touch with the pin, but for other vascular plant and bryophyte species, we recorded only one hit per pin. In the earlier study by Kitti, Forbes, and Oksanen (2009) bryophytes, ericoid dwarf shrubs, B. nana, and willows were counted similarly, but they counted multiple pin hits for the remainder of vascular plants that we counted by single hits. This discrepancy in counting pin hits in the different years potentially had a minor influence on the results; in our experience, multiple pin hits are not common and are often difficult to judge in a comparable manner. Species inside a subplot but not touched by the pin were given a value of 0.1 counts in all inventories. Hereafter we refer to point intercept data as “abundance.”

Growth and catkin production measurements in 2015

We measured the height of uppermost shoots of all separate willow shrubs (S. lapponum) and calculated their averages for each plot (n = 45), because we could not relocate earlier subplot measurements. This makes the measures of our hypothesized increase in willow height conservative, because the earlier measures depicted average maximum height at subplots, whereas we averaged over all willows in the plots. In addition, we collected ten to twenty willow leaves per plot from the midsection of the highest annual shoots. We pressed the leaves and used LAMINA (Bylesjö et al. 2008) to determine the mean leaf area (mm²). We also collected approximately ten annual shoots per plot. After drying at 70°C for 48 hours, we measured leaf nitrogen concentration (% of dry mass) from two leaf tissue subsamples (~0.01 g) per plot using LECO FP-528 Nitrogen/Protein Determinator.

We counted the number of fertile female plants and average number of flowering catkins (fruit bodies) per plant in each plot to study the reproductive success of willows. Because sex ratios in the genus Salix are often female biased (Crawford and Balfour 1983; Myers-Smith and Hik 2012) and we detected very few senescent male catkins, we counted only fertile female plants and catkins. Because we did not manage to relocate two of the original sites in Norway, we measured S. lapponum heights and collected leaves and catkins from two new sites that had general characteristics similar to those of the other sites.

Normalized Difference Vegetation Index and leaf area index

We compared NDVI and leaf area index (LAI) between grazing treatments to assess the impact of reindeer grazing on vegetation greenness. NDVI is a commonly used vegetation index that quantifies vegetation intensity from multispectral remote sensing imagery. NDVI is calculated between the near-infrared channel (which vegetation strongly reflects) and red channel (which vegetation absorbs). NDVI values were measured from WorldView-3 very high resolution satellite images (July 2015) with a spatial resolution of 1.24 m in multispectral channels (DigitalGlobe Inc., Longmont, CO). NDVI values were averaged over 9-pixel grids inside the exclosures (FI), over the freely grazed plots (FI), and in the original nongrazed plots (NO) that we managed to relocate in the field. For new nongrazed plots, where willows were measured in 2015, we averaged NDVI values over 49-pixel grids to make our hypothesized pattern of higher NDVI values in nongrazed fens conservative.

LAI (one-sided leaf surface area/ground area) indicates the amount of leaf area in an ecosystem. LAI can be retrieved from satellite data sets, and it is often used in monitoring tundra vegetation along with NDVI (e.g., Walker et al. 2003; Blok et al. 2011). In the Jauristunturit study area, LAI was measured indirectly in situ with an LAI-2200 Plant Canopy Analyzer (LI-COR Inc., Lincoln, NE). The device estimates LAI from light transmittance measurements above and below the canopy at five zenith angles (7°, 22°, 38°, 52°, and 68°) made with a fisheye optical sensor (LI-COR Inc. 2010). The instrument applies a model of radiative transfer in vegetative canopies to calculate the leaf area index. LAI was measured in grazed plots and exclosures in July 2016. The LAI value is a difference in light levels above and below the canopy. We took one measurement above and five below the canopy at four points in each plot (i.e., five above–below pairs per point). The LAI value for each plot is an average of twenty values (i.e., twenty above–below pairs).

Statistical analysis

Nomenclature follows Hämet-Ahti et al. (1998) for vascular plants and Laaka-Lindberg, Anttila, and Syrjänen (2009) for bryophytes. In cases where all individual pin hits were not identified to the species level or there were inconsistencies in species identification between years, we grouped plants to the genus level. Mosses of the family Mniaceae were combined into one group, and liverworts (Hepaticae) were treated as one group (except Aneura pinguis and Marchantia polymorpha). Due to difficulties in identification of sterile and hybrid Carex species in the field, C. limosa, C. magellanica, and C. rariflora were treated collectively (C. limosa coll.), and C. aquatilis ×bigelowii was combined with C. aquatilis.

Data availability and quality differed between the separate tests (Appendix 1). Vegetation patterns and the main directions of change of community composition were analyzed with nonmetric multidimensional scaling (NMDS) using combined 2002 and 2015 subplot-level species data sets. We excluded the 2006 data from NMDS because bryophytes were not surveyed that year. In order to standardize data between different inventories, we applied a Wisconsin double standardization transformation on vegetation data sets—that is, adjusting each species abundance to its maximum value and then adjusting all species values to proportions of plot total—and arcsine square root transformation of these values. Transformation was made separately for 2002 and 2015 vegetation data sets, which emphasizes changes in community composition and evens out possible differences between data sets and in personal perceptions of pin interceptions by different workers. NMDS was performed using the Sørensen (Bray-Curtis) distance measure. The effect of exclosure treatment on plant community structure was tested against free grazing treatment using a two-factorial permutation multivariate analysis of variance (PERMANOVA; Anderson 2001) using site as a second factor; that is, for five different mires. Because our interest was in grazing treatment, rather than random effects of mires, only grazing effects are reported. This was done for all species together and for bryophytes and vascular plants separately. We used the program PC-ORD 7.04 (McCune and Mefford 2016) for NDMS and PERMANOVA.

We tested differences between grazed sites, exclosures, and nongrazed sites in the abundance of plant functional and taxonomic groups using the 2015 point frame data. We divided the vegetation data into the following groups: all bryophytes, Sphagnum mosses, non-sphagnaceous mosses, hepatics, Salix spp., forbs, Carex spp., and Eriophorum spp. and averaged their point intercepts for each plot. On a species level, we tested the effect of grazing treatment levels on B. nana. We used a linear mixed effects model, in which a three-level grazing treatment was a fixed factor (grazed, exclosure, nongrazed) and site was included as a random factor. Approximate normality of the model residuals was assessed by graphical inspection of histograms and q–q plots, and a square root or logarithm transformation was applied when normality assumptions were violated (square root transformation for hepatics and Eriophorum spp; log transformation for Salix spp. and B. nana). Analyses were performed in R 3.4.3 (R Core Team 2017), and we used the “nlme” package’s (Pinheiro et al. 2017) lme() function for linear mixed effects models. The lsmeans package (Lenth 2016) was used to obtain t ratios and Tukey-adjusted p values for pairwise comparisons of the three grazing treatment levels when treatment was significant in the model.

The same model was used to test the abundance, height, leaf nitrogen concentration, and reproductive success of S. lapponum, as well as NDVI and LAI. We further compared S. lapponum height and Salix spp. (S. glauca, S. lapponum, and S. phylicifolia) total abundance over time (2002, 2006, and 2015). For these comparisons we included year in the model as a fixed factor and the interaction of year and grazing treatment level. Square root transformation was used for Salix abundance and number of female plants and log transformation was used for catkin number to approximate the normal distribution.

Results

Results in ordinations

The two-dimensional NMDS ordination had a high correlation between ordination distances and distances in the original n-dimensional space (R² = 0.482). Common minerotrophic moss species, such as Straminergon stramineum and Sphagnum teres together with the focal willow species Salix lapponum, characterized the bulk of the plots in a wide central area of the ordination. Certain mainly site-specific features added dispersion of plots in the ordination. Radiating toward different corners of the ordination, Sphagnum riparium, C. rostrata, Eriophorum angustifolium, and Aulacomnium palustre showed the most distinct patterns among species to separate groups of plots from the center of the ordination (see Appendices 2 and 3 for further details). All grazing treatments were overlapping (Figure 2), and we did not find any uniform directions of change in vegetation connected to the exclosure treatment; that is, the exclosure plots did not move toward nongrazed plots as expected (Appendix 4). In particular, grazed subplots (FI) and exclosures (FI) had widely overlapping ranges. The nongrazed Norwegian subplots had less variation in the NMDS than the Finnish exclosures and grazed subplots that were together spread more widely in the ordination. This separation of some grazed and exclosure plots apart from nongrazed Norwegian plots mainly reflected differences in shrub abundance. In the ordination, dwarf shrubs (Andromeda polifolia, B. nana, Empetrum nigrum ssp. hermaphroditum, and Vaccinium spp.) were more common in the Finnish sites that were separated from the nongrazed subplots in Norway, where abundance of Salix lapponum was high (Appendix 3). In addition, mosses characteristic of bog hummocks (Aulacomniun palustre, Dicranum spp., and Polytrichum spp.) and certain sedges (C. rotundata and Eriophorum angustifolium) were almost totally confined to some Finnish subplots.

The PERMANOVA analysis for Finland 2015 data with exclosures and grazed plots showed significant differences in plant community structure between treatments (F = 1.7747, p = .043), although NMDS failed to separate exclosures and grazed subplots. The analysis showed significant effects of grazing specifically on vascular plant communities (F = 2.4242, p = .02), whereas grazing did not affect bryophyte community structure (F = 1.1358, p = .319).

Abundance estimates of plant groups, NDVI, and LAI

We used the 2015 point intercept data to compare the abundances of plant groups between grazing treatment levels (Table 1, Figure 3). We found a significant effect of treatment on total abundance of bryophytes (p = .026), which was greatest in nongrazed fens. The effect was particularly pronounced in the case of hepatics (p = .004). They were significantly more abundant in nongrazed plots than in grazed plots (Tukey’s test, p = .038) and in the exclosures compared to grazed plots (Tukey’s test, p = .021). No significant differences were found between exclosures and nongrazed plots. On average, the abundances of Sphagnum mosses and non-sphagnaceous mosses were also greatest in nongrazed fens in Norway, though the effect of treatment was not significant.

Table 1. Analysis of variance results of treatment effects (Finland grazed, Finland exclosures, Norway nongrazed) from linear mixed effect models for abundances of different plant functional and taxonomic groups. Abundance was estimated in July 2015 using a point intercept method.

	df	F	p
Bryophytes	29	4.157	.026
Sphagnum mosses	29	0.137	.872
Non-sphagnaceous mosses	29	1.404	.262
Hepaticae	29	6.729	.004
Salix spp.	29	2.684	.085
Betula nana	29	1.820	.180
Forbs	29	0.026	.974
Carex spp.	29	1.590	.221
Eriophorum spp.	29	4.179	.025

The average abundance of Salix spp. was greatest in nongrazed plots, and it was doubled in the exclosures compared to the grazed plots. However, there was much variation, and this effect of treatment was not quite significant (p = .085). Grazing treatment level showed no significant effects on abundances of B. nana, forbs, or Carex spp. However, grazing treatment level did have a significant effect on the abundance of cottongrasses (mainly E. angustifolium; p = .025). On average, cottongrasses were least abundant in nongrazed fens in Norway.

No significant differences in NDVI and LAI were found. On average, however, NDVI was highest in nongrazed fens (Figure 3j). In Finnish sites, both mean NDVI and LAI (Figure 3k) were higher in exclosures compared to grazed plots, although the differences were not significant.

Figure 3. Abundances of different plant functional and taxonomic groups (a)–(i) in grazed plots, exclosures (Finland), and nongrazed (Norway) plots in July 2015. The abundance of plant groups has been estimated by a point intercept method. (j) Normalized Difference Vegetation Index in July 2015 and (k) leaf area index measured in July 2016. Leaf area index was measured only in exclosures and grazed plots (Finland).

Figure 4. Changes in mean (±SE) (a) Salix lapponum height and (b) Salix spp. abundance over the study period (2002, 2006, and 2015).

Figure 5. Comparisons of Salix lapponum (a) height (cm), (b) leaf area (mm²), (c) leaf nitrogen concentration (% of dry mass), (d) number of fertile female plants, and (e) number of flowering catkins, and between grazed sites and exclosures (Finland) and nongrazed sites (Norway) measured in 2015.

Abundance, size, and reproduction of Salix

As expected, both grazing treatment level (p < .001) and study year (p < .001) had significant effects on S. lapponum height (Table 2). A significant Year × Grazing interaction (p = .019) was also found. From 2002 to 2006 both grazed and exclosure plots in Finland had lower willow height and abundance than nongrazed fens in Norway (Figure 4). In 2015, the abundance of Salix spp. in the exclosures had reached the values of nongrazed plots, and grazed plots still had supressed abundance. Willow height showed a similar response, but the nongrazed treatment willows appeared also to have gained more height on average.

Table 2. Analysis of variance results from linear effect models for S. lapponum height and Salix spp. abundance over time (2002, 2006, and 2015).

						Interaction
		Treatment		Year		(Treatment × Year)
	df	F	p	F	p	F	p
S. lapponum height	79	16.031	<.001	23.361	<.001	4.177	.019
Salix spp. abundance	90	3.740	.028	0.302	.583	0.211	.810

In 2015, willows (S. lapponum) were significantly taller in nongrazed plots (NO) compared to grazed plots (FI; Tukey’s test, p = .0014; Table 3). They were also significantly higher in exclosures versus grazed plots (Tukey’s test, p < .0001), but there was no significant difference between exclosures and nongrazed plots. We did not find significant differences in the leaf area (cm²) of willows (p > .05; Figure 5) between treatments. However, there was a general significant effect of grazing treatment on leaf nitrogen concentration (p = .048). On average, leaf nitrogen concentration was lowest in nongrazed fens.

Table 3. Analysis of variance results of treatment effects (Finland grazed, Finland exclosures, Norway nongrazed) from linear mixed effect models for Salix lapponum height, leaf area, nitrogen concentration, and flowering traits measured in July 2015.

	df	F	p
Height (cm)	33	26.209	<.001
Leaf area (mm^2)	32	0.598	.556
Leaf nitrogen (% of dry mass)	32	3.345	.048
Flowering female plants (n)	33	10.253	<.001
Female catkins (n)	33	9.900	<.001

We found significant differences in the number of fertile female plants (p < .0001) and the number (p = .0004) and weight (p = .003) of female catkins between treatments. Fertile female plants were significantly more abundant in nongrazed plots in Norway compared to grazed plots in Finland (Tukey’s test, p = .0114). They also produced significantly more catkins in nongrazed fens than in grazed fens (Tukey’s test, p = .0166).

Discussion

Different types of wetlands form an integral part of northern boreal and arctic landscapes. However, changes in ecosystem structures, such as main vegetation characters, could alter the functioning of these important ecosystems. Reindeer grazing is one major factor to regulate northern plant communities and ecosystem balance. We examined the effect of reindeer grazing on sedge- and willow-dominated fens, a typical oroarctic wetland type in Finnish Lapland. By comparison between areas with historical differences in grazing and with experimental exclosures, we found minor effects of grazing on plant community structure as a whole but clear impacts on some functional plant groups and on the growth and reproduction of S. lapponum.

Grazing effects on mire plant communities

In the NMDS ordination, the nongrazed subplots in Norway had less variation in plant community structure, whereas they were overlapping with the Finnish grazed and exclosure subplots. Some grazed and exclosure subplots in Finland were, however, separated in the ordination from the nongrazed Norwegian plots, and a part of this separation may have resulted from more intensive grazing effects at these sites, and other site-specific factors may have contributed as well. In fact, grazing impact is not supported by the observed changes in these sites during the 13 years of reindeer exclusion, because the experimental plots did not move in any consistent way toward the nongrazed subplots in the ordination. Instead of grazing-controlled patterns, the ordination revealed gradients in mire vegetation related to variation in groundwater influence (minerotrophy) and wetness; that is, common gradients repeatedly found in mire vegetation analyses. Nevertheless, the permutation test of distance matrix between exclosures and freely grazed plots did reveal a significant effect. These results indicate that in oroarctic wetlands, reindeer have only limited influence on plant community structure and that site-specific hydrological factors are more important. Instead of overarching impact on community structure as a whole, we found that reindeer grazing effects were confined to certain plant groups.

Bryophytes, particularly hepatics, were significantly more abundant in nongrazed fens in Norway than in grazed fens in Finland. Mosses constitute a less-preferred food source for reindeer (Danell et al. 1994), but Sphagnum mosses and hepatics are particularly sensitive to trampling (Studlar 1980; Arnesen 1999). Presumably, trampling has negatively affected bryophyte cover in grazed fens in our study area, but redistribution of nitrogen via feces and urine can also have a negative impact on moss layer (van der Wal et al. 2004). Indeed, in wet high-arctic tundra habitats, grazing and trampling have been observed to reduce moss layer depth (van der Wal, van Lieshout, and Loonen 2001), total bryophyte cover (van der Wal et al. 2007), and biomass (Falk et al. 2015; Mosbacher et al. 2019). In dry tundra habitats, exclusion of reindeer grazing often leads to an increase in the cover of both deciduous and evergreen shrubs at the expense of mosses (Pajunen, Oksanen, and Virtanen 2011; Vowles, Gunnarsson et al. 2017).In our study sites, the pattern of bryophytes differed from those found in previous studies, because they were most abundant in nongrazed fens, where willows have also increased in size and abundance.

The reduction in moss layer depth can increase soil temperature and support graminoid growth in arctic moss-dominated environments (van der Wal, van Lieshout, and Loonen 2001; van der Wal and Brooker 2004). Rhizomatous graminoids, such as sedges (Carex spp.) and cottongrasses (Eriophorum spp.), are tolerant to trampling (Arnesen 1999; Forbes, Ebersole, and Strandberg 2001) and may benefit from reindeer grazing (Henry 1998; Olofsson et al. 2001; Eskelinen and Oksanen 2006; Ravolainen et al. 2011). In our study area, reindeer grazing has maintained sedge-dominated vegetation, and the abundances of Salix spp. and B. nana were lowest in grazed areas, on average. Cottongrasses (mainly E. angustifolium) were less abundant in nongrazed fens in Norway compared to all plots in Finland, but we did not find significant differences in total abundance of sedges (Carex spp.) between treatments. At the individual species level, C. rotundata and E. angustifolium were almost totally confined to Finnish sites, whereas C. rostrata was most abundant in nongrazed plots in Norway. Sedges are eaten by reindeer and some of these differences thus may be caused by reindeer grazing. This remains highly speculative, however, because it is possible that some other differences between the areas affect abundance of these species. Falk et al. (2015) showed that a short-term (3-year) exclusion of muskox (Ovibos moschatus) decreased the abundance of Eriophorum scheuchzeri tillers in a high-arctic mire. Reindeer feed highly selectively in summer (Nieminen and Heiskari 1989), and plant species may respond differently to grazing when some species are selected over others.

Grazing effects on Salix lapponum

Though reindeer grazing effects on mire plant community structure were subtle, grazing evidently affects the growth and reproductive success of S. lapponum. Many previous studies have shown that reindeer grazing reduces the growth of tundra willows (e.g., den Herder, Virtanen, and Roininen 2004, 2008; Pajunen 2009), and 10 to 13 years of herbivore exclusion is sufficient to result in a significant increase in tundra shrub cover (Pajunen, Virtanen, and Roininen 2012). As expected, reindeer herbivory in summer suppresses the growth of S. lapponum in oroarctic mires as well. Results from a short-term study (Kitti, Forbes, and Oksanen 2009) of the same sites showed no significant effect of 4-year exclusion of reindeer grazing on S. lapponum height, but after 13 years, willows had recovered inside the experimental exclosures. In addition to growth response, willows produced less catkins in grazed fens compared to exclosures and nongrazed fens, which reduces the capacity of willows to colonize new areas. Reindeer use developing willow catkins as a food source (Klein 1990), but the lower number of catkins in grazed fen areas likely relates to decreased growth and indicates generally poor fitness of willows. Lower height and number of catkins in grazed fens indicates that S. lapponum responds to grazing by allocating resources to vegetative growth at the expense of sexual reproduction. Brookshire et al. (2002) observed that willows (S. boothii, S. geyeriana) subjected to grazing by wild ungulates and sheep did not flower until 2 years after herbivore exclusion. In our study area, reproductive success was studied only in 2015; therefore, we can only conclude that the 13-year time span was sufficient for revival of sexual reproduction of S. lapponum after reindeer exclusion.

Ungulate herbivores also affect ecosystem N cycling though several processes (Singer and Schoenecker 2003). Reindeer deposit feces and urine, which may accelerate N mineralization rates and increase N availability for plants (Olofsson et al. 2001; Stark, Strömmer, and Tuomi 2002). Increased N availability can be reflected in higher N concentration in plant leaves (Ritchie, Tilman, and Knops 1998). Indeed, we found that S. lapponum growing in the grazed fens in Finland had higher leaf N concentrations compared to willows growing in nongrazed fens in Norway. This suggests that reindeer grazing could increase the nutritional value of willows, which may lead to selection of previously grazed plants. Indeed, reindeer feeding is elsewhere observed to be more intense on small, young willow sprouts compared to old willows (den Herder, Virtanen, and Roininen 2004). This pattern in leaf N content may also reflect resource allocation between compensatory growth and chemical defense. Plant species have two alternative strategies against herbivory: (1) a well-developed regrowth capacity and a poor defense mechanism and (2) a well-developed defense mechanism and a poor regrowth capacity (van der Meijden, Wijn, and Verkaar 1988). Skarpe and Wal (2002) found that simulated browsing increased leaf N content but had no effect on the concentration of phenolic compounds in the leaves of deciduous dwarf shrub Salix polaris, indicating that under herbivory, S. polaris allocates recourses to growth rather than defense. In the absence of grazing, willows may in turn invest more on defense compounds (Dormann and Skarpe 2002), resulting in lowered palatability, and allocate nutrients in reproduction rather than vegetative growth. In our study, we did not measure defense compounds and therefore it remains speculative whether reindeer affect the chemical defense of S. lapponum. However, in every case, elevated N concentration in S. lapponum leaves was still observed 13 years after reindeer exclusion, indicating that grazing has a long-lasting impact on N balance of S. lapponum.

Shrubification of oroarctic mires

Expansion of deciduous tall shrubs in arctic tundra landscapes is a well-documented phenomenon (e.g., Myers-Smith et al. 2011, 2015; Vowles and Björk 2018). Reindeer grazing can inhibit climate change–induced shrub expansion (Olofsson et al. 2009; Ravolainen et al. 2014) and hence act to maintain resilience in ecosystems (Kaarlejärvi, Hoset, and Olofsson 2015). In the Finnish side of our study area, reindeer grazing evidently controls the growth and sexual reproduction of willows, and reindeer can be regarded as a bioengineer that regulates shrub expansion. In the absence of this control, the height of willows increased in nongrazed fens in Norway during the study period. It is likely that climate change has accentuated the growth of willows, because climate records show significant warming in northern Lapland for the past few decades (Kivinen et al. 2017). Soil moisture is another key factor affecting shrub growth in tundra habitats (Myers-Smith et al. 2015). Due to lack of temperature-induced moisture limitation, climate change–induced expansion of deciduous shrubs could be particularly rapid in moist and wet habitats, such as riparian sites (Ravolainen et al. 2011; Ackerman et al. 2017). Indeed, no evidence of larger-scale expansion of tall deciduous shrubs in Fennoscandian forest–tundra habitats exists, despite climate warming (Vuorinen et al. 2017; Maliniemi et al. 2018).

In this study, we combined remote sensing data with field measurements of Salix growth and abundance. The difference in S. lapponum height and abundance between the exclosures and grazed fens was clear and was also visible in some drone images taken in 2016 and in digital elevation models derived from orthomosaics (Timo Kumpula, unpublished; Appendix 5). However, we did not find significant differences in NDVI values derived from very high-resolution WorldView-3 images of exclosures and grazed areas. Despite the evident changes in S. lapponum growth and abundance in exclosures, it appears that even 13 years is not long enough to see significant differences in NDVI between grazed and nongrazed fens. Another possible explanation is that we could not detect differences in NDVI because sedge-dominated grazed fen areas have dense green vegetation as well. However, the results of NDVI were confirmed by the observation that LAI values did not show significant differences between exclosures and grazed areas. The increased height and number of flowering catkins suggest that 13 years of reindeer exclusion was enough for recovery of individual willows from grazing, though actual greening was not realized as yet. Furthermore, because we did not find significant difference between Finnish and Norwegian plots, it remains uncertain whether satellite-derived NDVI values of fens may indicate higher shrub cover in Norway, where reindeer summer grazing has been prohibited since the late 1950s. Cohen et al. (2013) compared temporal variation in NDVI values between Finnish summer pastures and Norwegian winter pastures, and one of their test areas overlaps with our study sites. They found that when the snowmelt started in the test area, NDVI increased faster on the Norwegian side, indicating denser shrub vegetation. However, later during the growing season, NDVI values were higher in the Finnish summer pasture. They included all vegetation types in their analysis, and higher NDVI values in summer pasture were probably due to higher lichen cover in the Norwegian winter pasture. In our study, we compared only fens, and the lack of inference is possibly due to the low number of sites. A different approach of mapping different types of mire areas across the border fence for NDVI could help to unravel how this most commonly used index performs in observing shrubification in northern mires.

Reindeer grazing can affect also ecosystem carbon dynamics and modulate ecosystem responses to climate change in boreal forest and arctic tundra habitats (e.g., Väisänen et al. 2014; Köster et al. 2018; Ylänne et al. 2018). However, studies focusing on grazing impacts on carbon dynamics particularly in mire ecosystems are still few and were mainly conducted in high-arctic mires with continuous permafrost (Falk et al. 2014, 2015; Mosbacher et al. 2019). It is difficult to draw conclusions regarding grazing impacts on carbon dynamics, because processes and mechanisms affecting ecosystem carbon dynamics are complex and largely dependent on habitat type. For example, in shrub and graminoid tundra in Greenland, an ecosystem CO₂ sink function became stronger after exclusion of muskox (O. moschatus) and followed the increase of tall shrubs (Cahoon et al. 2012), whereas in sedge- and grass-dominated mires, grazing exclusion resulted in a decrease in the CO₂ sink function (Falk et al. 2015). Reindeer can also increase methane emissions from mires via fecal deposition (Laiho, Penttilä, and Fritze 2017). The impacts of reindeer grazing on carbon dynamics in willow-dominated mires remain uncertain. However, considering the special role of bryophytes (especially Sphagnum) in carbon dynamics of northern mires, our observation of reduced abundance of bryophytes in freely grazed mires indicates that reindeer grazing might suppress formation of peat and, hence, carbon sink. In line with this expectation, bryophytes specifically increased the carbon pool in response to muskox exclusion in arctic mires (Mosbacher et al. 2019).

Conclusion

Fennoscandian oroarctic mire plant communities as a whole are resistant to reindeer grazing effects that could only be revealed with detailed individual species and species group abundance ca. 55 years after cessation of grazing. Reindeer grazing particularly affects the abundance of bryophytes and controls the abundance, size, and flowering of S. lapponum. Individual willows are able to tolerate moderately intensive reindeer grazing, however. In grazed fens, willows are still present and able to recover after grazing removal. In addition, nitrogen concentration of S. lapponum was elevated in grazed fens after 13 years of reindeer exclusion. Despite the recovery of willows inside the exclosures, 13 years of grazing exclusion was not enough to result in changes observable by NDVI and LAI, and greening in terms of increased vegetation cover was not realized in exclosures. These results indicate that early phases of vegetation changes potentially leading to shrubification of oroarctic and arctic mires may take place undetected even by the most accurate satellite-derived NDVI.

Acknowledgments

We thank Heidi Kitti and Jari Oksanen for providing data from 2002 and 2006. We thank Karoliina Heiskanen, Pasi Korpelainen, Andrew Martin, Franziska Rosaroja, and Veikko Tahvanainen for their assistance in fieldwork. Henni Ylänne gave valuable advices for the selection of statistical analyses. Juha Magga and Samuli Näkkälä are thanked for invaluable advice and help with many practical issues during fieldwork. Comments from the editor and three anonymous referees significantly improved the article.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary material

Supplemental data for this article can be accessed on the publisher’s website.