Anthropogenic disturbance modifies long- term changes of boreal mountain vegetation under contemporary climate warming

Aims: Accelerating high- latitude climate warming drives shrub expansion in open landscapes and alters species distributions and compositions of plant communi -ties. Simultaneously, various land use practices cause disturbance to the vegetation. However, not much documentation exists on how long- term intensive land use disturbance modifies high- latitude vegetation under climate warming. Here, we study how the composition of boreal mountain plant communities has changed during three decades in response to heavy land use disturbance, related to ski resort construction and management, and how these changes compare to those observed in adjacent less disturbed communities. Location: Iso- Syöte, Finland. Methods: We resurveyed vegetation along four elevational gradients (240– 426 m a.s.l.) on a boreal mountain in 2013– 14.

Particularly rare are the studies that consider both the entire community composition (canopy, shrub, field and ground layers) and cover over a sufficiently long time period (decades). These aspects are of high importance in order to better understand the often slow temporal dynamics of high-latitude plant communities with a generally high share of ground-layer species (Callaghan et al., 2004;Maliniemi et al., 2018Maliniemi et al., , 2019. One widely distributed form of local land use that often causes severe disturbance in the vegetation, is the construction and management of tourist attractions, such as ski resorts (Roux-Fouillet et al., 2011;Holtmeier & Broll, 2018;Bacchiocchi et al., 2019). In forested areas, this typically means a reduced or completely removed tree layer while the understorey vegetation is retained to varying degrees. The clearing of the tree layer strongly alters the microclimate (De Frenne et al., 2019), allowing solar radiation and wind flow near the soil surface (Holtmeier & Broll, 2005). Ski piste management and the use of artificial snow may affect snow cover duration, chemical properties of the soil, freezing conditions during the cold season and thermal conditions throughout the year (Rixen et al., 2004;Roux-Fouillet et al., 2011;Bacchiocchi et al., 2019). The managed snow cover and mechanical damage, e.g., due to grooming, has also been shown to cause dieback of dwarf shrubs (Stoeckli & Rixen, 2000).
Given the relatively slow process rates and long recovery time of high-latitude vegetation, these kinds of disturbances with relatively short intervals and large areal extent can destabilise the system and shift the compositional trajectories of vegetation (Turner et al., 1993). Therefore, by strongly altering site conditions, disturbances due to ski resort infrastructure may also be expected to considerably influence climate-driven vegetation changes, such as shrubification.
Shrub expansion is one of the most evident responses to climate warming in open high-latitude landscapes (Myers-Smith et al., 2011;Scharnagl et al., 2019). However, improved light conditions and enhanced wind dispersal of seeds due to tree clearance, longer-lasting burial under snow, or certain levels of soil disturbance may alone facilitate the establishment of good competitors, such as taller shrubs (Myers-Smith et al., 2011;Tilman, 1988). Importantly, warmer conditions can further accelerate the growth response of shrubs to the increased light (De Frenne et al., 2015). Thus, the combined effects of warming and certain levels of disturbance can be expected to enhance shrub abundance, unless the disturbance is too severe and begins to limit shrub growth (Verma et al., 2020). Increasing shrub cover is further associated with changes in local abiotic conditions (e.g., hydrology and nutrient cycling) but also to changes in understorey species composition, e.g., reduced cover of small species and species richness by plant-plant interactions (Pajunen, Oksanen & Virtanen, 2011;Wallace & Baltzer, 2020).
Late-successional forest communities in European boreal mountains are often relatively species-poor and dominated by Norway spruce (Picea abies), dwarf shrubs (Vaccinium spp.) and a moss layer with a few dominant species (Hylocomium splendens, Pleurozium schreberi). However, it is well known that natural or anthropogenic disturbances (e.g., fire or various silvicultural practices) can, at least temporarily, increase local species richness in boreal forests (Peltzer et al., 2000;Kumar et al., 2018). In arctic-alpine landscapes, in turn, disturbances due to ski resort management have been shown to have no or negative effects on species richness (Bacchiocci et al., 2019).
In undisturbed boreal mountain landscapes, tree and shrub cover thins out towards higher elevations giving way to open slopes and mountain tops that can maintain a tundra-like vegetation. Although it can be generally expected that species shift towards higher elevations in response to a warmer climate (Lenoir et al., 2008;Steinbauer et al., 2018) and that the shift is stronger the lower species were situated historically (Rumpf et al., 2018), it has been found possible that disturbances modify these range shifts. According to Lenoir et al. (2010), disturbance due to habitat modification and the consequent release in competition could also cause a downward shift of species that are filling only part of their potential distributional range along the elevational gradient. Agreeing with this, Guo et al. (2018) found that loss of forest cover caused downward shifts in species' elevational distributions Therefore, along with the climatic effects, land use effects on vegetation are important to consider in management actions and in projections of future vegetation. suggesting that land use disturbances may complicate recognition and interpretation of climate-driven vegetation changes.
Here, we study how heavy land use disturbance, due to the construction and management of a ski resort, influences long-term vegetation changes along elevational gradients and compare the observed changes to those occurring along only slightly disturbed gradients. In 2013-14, we resurveyed plant communities along four elevational gradients on a boreal mountain in northern Finland that were originally surveyed in 1980-81. Each gradient started from a closed forest, changed gradually into more open slopes, and reached the nearly treeless top. Right after the original survey, half of the slopes and the whole top were partly cleared of trees to accommodate ski resort infrastructure and were thus subjected to a relatively heavy disturbance that continued until the resurvey. In contrast, the other half of the slopes was located outside the resort and was only slightly disturbed. This unique setting allowed us to analyse longterm impacts of different disturbance levels on the boreal mountain plant communities within a similar environmental context and under a similar impact of macroclimatic warming. More specifically, we analysed temporal changes in the cover of different plant groups, species richness and species' elevational range means in relation to heavy vs slight disturbance and elevation.
We hypothesised that heavy disturbance has a strong influence on the long-term vegetation changes that are expected under a warmer climate. Under slight disturbance, we predicted to observe changes that are frequently reported as a response to warmer conditions, such as advancement of taller shrubs on the originally open upper slopes and upward movement of low-elevational species. In contrast, we predicted that heavy disturbance and the tree cover removal result in an increase of taller shrubs along the whole elevational gradient and cause both upward and downward movements of species in the long term. We further predicted that heavy disturbance reduces the cover of dwarf shrubs and ground layer species (bryophytes and lichens) but enhances species richness due to species introductions.

| Study area and vegetation resurvey
The study area, Iso-Syöte ( The NE vegetation gradient was situated in the immediate vicinity of the ski routes and lifts, while the SE gradient was situated next to paved roads, numerous cabins and a hotel. Although the specific disturbance types are partially different between these gradients, they generally represent heavy mechanical and physical land use disturbance. For most parts of the gradients, trees were cut or thinned out and forest understorey vegetation was typically disturbed by infrastructure maintenance or building yet was not completely removed (Appendix S1). No artificial revegetation or seeding was done during the study period. Based on the field observations, no prominent or directional variation in disturbance intensity was observed along the slopes are regarded as "heavily disturbed" gradients. In an apparent contrast, the SW and NW slopes (250-395 m a.s.l.) remained nearly unbuilt over time (Figure 1c, d, Appendix S1). Only a few hiking trails crossing parts of the slopes were established during the study period. The vegetation on these slopes was protected from any kind of management since 1910 and for the most part also prior to this (Metsähallitus, 2004). A visual inspection on site during the resurvey confirmed the absence of land use disturbance. Thus, these gradients are regarded as "slightly disturbed." Sample plots at the heavily disturbed top are treated separately in the analyses. Thus, the heavily and slightly disturbed slopes have comparable elevational range.

MALINIEMI ANd VIRTANEN
The vegetation plots in this study can be regarded as quasipermanent (Kapfer et al., 2017) and thus, the comparisons of vegetation between the surveys are affected by a relocation error. This may cause pseudo-turnover in species composition, which adds random error to the observations of temporal change in vegetation.
The presence of random error may increase the risk of not detecting real temporal changes yet is unlikely a source of a systematic error (Kopecký & Macek, 2015;Kapfer et al., 2017). In our case, considering the relatively detailed information on gradient and plot locations, we estimated that the relocation error for the 1 m × 1 m plots is generally <20 m. Due to the larger plot size, the estimates of shrub and canopy layers and of species richness are less affected by the relocation error.
Original and resurveyed species data were taxonomically harmonised before the analyses, i.e. synonyms were identified, and certain bryophyte and lichen species were treated collectively (Appendix S2). Because the studied vegetation is composed of relatively few species, most of which are common and identifiable in the field, we find it unlikely that the observed vegetation changes would reflect systematic error between the observers. To analyse temporal changes, species were grouped into the following morphological groups: shrubs (including both tall and low-growing shrubs), dwarf shrubs, graminoids, herbs, bryophytes, and lichens. Absolute cover values (%) were calculated for each group by summing the cover values of species in the group. In order to detect elevational shifts of biogeographical species groups, vascular plants were classified into four different general biogeographic categories: temperate-boreal, boreal, boreal-hemiarctic or hemiarctic after Engelskjøn, Skifte and Benum (1995). A class could be assigned for 93% (67 out of 72 species) of the vascular plant species (Appendix S2).

| Environmental data
E-OBS raster data (10 km × 10 km; Cornes et al., 2018) were used to calculate temporal changes in mean annual temperatures, annual precipitation and annual thermal sums (summed from daily temperatures exceeding the +5°C threshold during the growing season) in the Iso-Syöte area. Data on annual reindeer numbers in the herding district containing Iso-Syöte were obtained from the Natural Resources Institute Finland and converted into grazing pressure (reindeer/km 2 ).
Temporal differences in climate and grazing pressure variables were calculated using 15-year averages prior to each survey.

| Statistical analyses
To have an overall view on the long-term compositional changes in communities under slight and heavy disturbance, we used nonmetric multidimensional scaling (NMDS) ordination, based on Bray-Curtis distances and plot-level data (n = 43) on species' covers, from the R package vegan (Oksanen et al., 2019). Stable and adequate stress solution (stress > 0.2) was reached using three dimensions (k = 3). Correlation vectors of elevation and shrub cover were fitted to the ordination to find out if compositional changes correlate with elevation or the expected increase of shrubs. The goodness-of-fit of the fitted vectors was estimated using 999 permutations. Due to repeated sampling, permutations were not allowed within paired plots.
We used generalised linear mixed-effects models (GLMM) to analyse the effects of time, disturbance level and elevation on vegetation changes. The first set of models was built for the slope data that are comparable in elevational gradients between different disturbance levels (240-395 m a.s.l.). We built separate models for seven cover (%) variables (trees, shrubs, dwarf shrubs, graminoids, herbs, bryophytes and lichens) and four species richness variables (total, vascular plant, bryophyte and lichen richness). Time (original survey vs resurvey), disturbance level (slight vs heavy), elevation (m a.s.l.) and their interactions were added as fixed factors in the models. We were particularly interested in whether the response variables had changed over time (T) and whether these changes were dependent either on the disturbance level (T × D) or elevation (T × E) or on both the disturbance level and elevation (T × D × E). Plot identity was added as a random factor to control for the repeated sampling. Aspect could not be used as a random factor because half of the potential combinations were lacking (i.e. neither level of disturbance was found in all aspects). However, none of the disturbance levels were entirely directed to either north or south. This reduces the possible influence of the aspect together with the long summer days and low sun angles, which have been shown to even out differences between the aspects to some extent in the study area (Winkler et al., 2016). Model fits and normality of residuals were examined using diagnostic plots. There was no substantial spatial autocorrelation in the residuals of any model (Appendix S3). All cover variables were modelled with a Gaussian distribution. Four variables (shrub, graminoid, herb and lichen cover) were log(x + 1)-transformed to improve the model fit and the predictions of these models were back- (original survey vs resurvey) was added as a fixed factor. Plot identity was added as a random factor to account for repeated samplings.
All the models were run using the R package glmmTMB (Magnusson et al., 2020).
We used the slope data to test whether temporal shifts in species' elevational range means were influenced by the level MALINIEMI ANd VIRTANEN of disturbance and species' mean elevational position during the original survey (e.g., whether low-elevational forest species shifted their range means more than species with originally higher elevational position). We first calculated weighted elevational range mean for each species (including vascular plants, bryophytes and lichens), on both slightly and heavily disturbed slopes during both surveys. Species' abundances (percent cover) were used as weights. We then built a generalised linear mixed-effects model to explain species' elevational means during the resurvey by their elevational means during the original survey, disturbance (slight vs heavy) and their interaction. Species identity was added as a structural term in the random-effects part of the model. Only species that were present at least three times during both surveys were included in this analysis. Model fit and normality of residuals were examined using diagnostic plots. Additionally, we estimated temporal shifts in the weighted mean elevations of biogeographic plant groups under both disturbance levels using t tests for weighted means from the R package weights (Pasek, 2018). All statistical analyses were performed using the R software (R Core Team, 2019).

| RE SULTS
Mean annual temperature, as a 15-year average prior to each survey, rose by 1.5°C during the study period (from 0.08 to 1.61°C), while annual precipitation increased only marginally (c. 610-658 mm, Appendix S4a, b). Annual thermal sums increased 16% between the surveys (from 878 to 1,044°Cd, Appendix S4c). Grazing pressure remained relatively low and constant (c. 1.5 reindeer/ km 2 , Appendix S4d). The mean ± SE soil pH in the resurvey was 3.67 ± 0.04 on slightly and 4.24 ± 0.11 on heavily disturbed gradient.
The NMDS ordination suggested that the original slope communities (which were later subjected to different disturbance levels) had relatively similar composition, while the original composition of top communities differed from that of slope communities (Figure 2).
However, under long-term exposure to different disturbance levels, the composition of slope communities became dissimilar in relation to their original compositions. Expectedly, the NMDS suggested that the shift was stronger under heavy disturbance. Moreover, the composition of top communities became more similar to the composition of slope communities.

Generalised linear mixed-effects model fits indicated that
the elevational patterns of all the plant group covers were originally rather similar between the slopes that were later subjected to either slight or heavy disturbance (Figure 3). However, these patterns changed prominently during the three decades. The expected decrease in canopy cover over time was prominent under heavy disturbance (Table 1; T × D), being strongest at lower elevations ( Figure 3a). Surprisingly, a similar temporal change in the elevational pattern was observed also under slight disturbance (Table 1; T × E, Figure 3a). The cover of shrubs clearly increased over time, but this response was affected by both the disturbance level and elevation: shrubs increased towards the top under slight disturbance, whereas the increase was more uniform along the heavily disturbed slopes (Table 1; T × D × E, Figure 3b). The covers of dwarf shrubs, graminoids and herbs changed over time in terms of their elevational distribution (Figure 3c-e, Table 1; T × E).

F I G U R E 2
Compositional changes over time. The non-metric multidimensional scaling (NMDS) ordination shows the composition of communities subjected to slight (SD) and heavy disturbance (HD) and the corresponding original communities. Ellipses are 1 standard deviation for group centroids. Temporal shifts in the community composition seem most pronounced under heavily disturbed slopes and at the top. The composition of top communities has shifted towards the composition of more forested slope communities over time. Temporal shifts in community compositions correlate with increased shrub coverage (r 2 = 0.28, p = 0.001) and elevation (r 2 = 0.44, p = 0.001; ordination plots with the third axis in Appendix S5). Stress = 0.15 | 7 of 13

MALINIEMI ANd VIRTANEN
Non-overlapping confidence bands indicate that dwarf shrubs also gained ground at lower slopes under slight disturbance (Figure 3c), whereas graminoid increase was strongest at the upper slopes under heavy disturbance (Figure 3d). Bryophyte cover decreased strongly at lower slopes, regardless of the disturbance level, leading to temporal changes in their elevational distribution (Table 1: T × E, Figure 3f). Lichen cover tended to decrease only slightly on the heavily disturbed upper slopes (Figure 3g), yet this change was not significant. Total species richness increased prominently under both disturbance levels over time but even more so under heavy disturbance (Table 1; T × D, Figure 4a). Here, the increase of total species richness was primarily due to the increased number of vascular plants (Table 1; T × D, Figure 4b). Bryophyte richness increased under both disturbance levels (Table 1: T, Figure 4c) and lichen richness only little on the slightly disturbed upper slopes (Table 1; T × D, Figure 4d). At the heavily disturbed top, canopy cover decreased over time, while shrub cover, total species richness, vascular plant richness and bryophyte richness increased (Table 2).
Species' elevational position during the original survey and the level of disturbance had significant effects on species' weighted mean elevation during the resurvey (X 2 = 19.8, p < 0.001 and X 2 = 9.5, p = 0.002, respectively; Figure 5). However, the magnitude of change in mean elevation was different between the disturbance levels (original mean elevation × disturbance level: X 2 = 5.5, p = 0.019). Species with an originally lower mean elevation shifted upwards over time regardless of the level of disturbance (confidence bands diverge from the 1:1 line that represents no temporal change in mean elevation in Figure 5). Many species with originally higher mean elevation shifted downwards under heavy disturbance, whereas such a trend was not detected on slightly disturbed slopes.
Additional analyses for biogeographical groups showed that the weighted mean elevation of temperate-boreal, hemiarctic-boreal and hemiarctic vascular plant groups shifted upwards under slight disturbance, whereas under heavy disturbance, the mean elevation of the hemiarctic group shifted downwards (Appendix S7).

| D ISCUSS I ON
We resurveyed the vegetation of a boreal mountain along four elevational gradients that had been originally surveyed before the construction of ski resort facilities in the 1980s. Between the surveys, climate warmed in the study area for over two decades, which was expected to affect the vegetation. Half of the gradients were taken for ski resort facilities, while the other half remained relatively F I G U R E 3 Fits of generalized linear mixed-effects models with 95% confidence bands for the cover (%) of (a) canopy, (b) shrubs, (c) dwarf shrubs, (d) graminoids, (e) herbs, (f) bryophytes and (g) lichens along the elevational gradient on slightly and heavily disturbed slopes during the original survey and the resurvey. Model summaries (including estimates, standard deviations and pseudo-R 2 values) are in Appendix S6 MALINIEMI ANd VIRTANEN intact, allowing comparison of the effects of heavy and slight land use disturbance on vegetation changes that had taken place during 34 years. We observed differences in shrub encroachment, species richness and species' elevational shifts between slightly and heavily disturbed gradients. These findings therefore generally support our hypothesis that the level of land use disturbance considerably influences the long-term vegetation changes in a boreal landscape.
In the following, we discuss the main findings to establish the role of  (Harsch et al., 2009). Even though we cannot completely rule out other drivers that may have had effect on these observed changes (e.g., natural succession), the encroachment of trees and shrubs seem to rapidly transform formerly semi-open or open boreal heaths towards more closed forest stands. In contrast, tree cover clearance and continued heavy disturbance tended to limit shrub expansion on the upper slopes and to promote a uniform, yet modest, increase along the elevational gradient. Thus, our results indicate that heavy land use disturbance may limit major expansion of shrubs on formerly open habitats to some degree but can allow their presence at low to moderate abundances under cleared or thinned tree canopies. Taken together, these findings strongly suggest that increases in tree canopies and shrub abundance could be more common in the boreal mountain uplands if there were no constraints due to disturbances.
Contrary to our expectation, the overall cover of dwarf shrubs, mainly composed of ericoids (Vaccinium spp.), was not notably affected by the disturbance level and thus, seemed to resist relatively heavy disturbance. Like taller shrubs, dwarf shrubs also may benefit from lighter and warmer conditions (Myers-Smith et al., 2011). Thus, under reduced canopy shade due to heavy disturbance, they could retain at least part of their abundance. According to Zeidler et al. (2016), dwarf shrubs can also benefit from managed snow cover and perform well on ski slopes, which may partly explain the observed results. The increase of dwarf shrubs on slightly disturbed lower slopes coincides with the reduced tree cover. The latter is likely due to a natural dieback of old-growth trees and canopy damage due to pronounced crown snow-loads typical for the study area (Gregow et al., 2008).
There was a notable overall decrease in bryophyte cover between the surveys that depended on elevation. Most pronounced decreases took place on the lower slopes under slight disturbance, coinciding with the increased dwarf shrub cover on slightly disturbed slopes and with increased shrub cover on heavily disturbed slopes. Decreased bryophyte cover may also be a general response to a warmer microclimate, solar radiation damage and evaporation stress due to decreased canopy cover (Busby et al., 1978). The TA B L E 1 Results from generalized linear mixed models testing the effects of time (original survey vs resurvey), disturbance (heavy vs slight), elevation (m a.s.l.) and their interactions on different plant group covers (%) and species richness. df = 1 other main ground layer element, lichen cover, tended to decrease towards the top under heavy disturbance, which can be expected because fruticose lichens, especially in dry conditions, are sensitive to disturbances and recover slowly (Klein & Schulski, 2009;Heggenes et al., 2017).

Elevation (E)
Species' elevational shifts, in terms of their elevational range means, were directed primarily upwards under slight disturbance, which was expected given the warmer conditions (Lenoir et al., 2008;Parolo & Rossi, 2008). Species-specific shifts were stronger in those species that had originally lower elevational position, as also shown by Rumpf et al. (2018). Such clear upward shifts of species and also certain biogeographical groups (Appendix S7) may be highlighted in boreal landscapes where low-lying and short elevational gradients along with the proximity to closed forests may allow relatively rapid movement of species compared to, for instance, alpine environments. However, as anticipated based on Lenoir et al. (2010), our results showed that heavy disturbance promoted several species to expand their elevational range means also downwards. This likely reflects an increased immigration of local species following more open and partially disturbed forest floor vegetation and may also explain the absence of upward shifts in biogeographical group means under heavy disturbance (Appendix S7). Taken together, these patterns strongly suggest that disturbance opens new space for species immigration and thus, contributes to more pronounced elevational shifts of individual species in boreal landscapes.
As expected, the ski resort as a heavy disturbance agent enhanced plant richness in the boreal mountain environment, promoting richness particularly at the lower slopes compared to slight disturbance. crucial, since year-round recreational activities will likely increase due to growing tourism (Saarinen, 2007).
The negative effects of ski resorts on mountain environments are relatively well known, being strongest under soil-grading and revegetation by invasive species (Tsuyuzaki, 1994B;urt & Rice, 2009;Roux-Fouillet et al., 2011). Thus, the realisation of a potentially positive effect is an intriguing idea. We find it possible that disturbances owing to mainly tree clearance, but preserving soils and field layer vegetation to some extent, may contribute, at least temporarily, to preservation of species requiring open habitats (see also Steinbauer et al., 2018;Chardon et al., 2019)  we emphasise that the disturbance effects depend on the context (e.g., elevational position) and intensity, and these factors, along with the multiple negative effects that human land use disturbance imposes on the environment must always be considered when developing conservation and management plans.

ACK N OWLED G EM ENTS
We thank Karoliina Huusko and Anu Skog for assisting in the resurvey. John-Arvid Grytnes gave insightful comments that improved the manuscript. Jouko Kumpula helped acquiring the reindeer data. We acknowledge the E-OBS dataset from the EU-FP6 project UERRA (http://www.uerra.eu) and the Copernicus Climate Change Service, and the data providers in the ECA&D project (https:// www.ecad.eu).

AUTH O R CO NTR I B UTI O N S
TM and RV conceived the ideas and designed methodology; TM collected and analysed the data and led the writing of the manuscript.
Both authors contributed critically to the drafts and gave final approval for publication.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data used in this study will be deposited in EVA (European Vegetation Archive).

F I G U R E 5
Generalized linear mixed-effects model fit with 95% confidence bands for species' weighted mean elevation in the resurvey in relation to their original weighted mean elevation and the disturbance level (SD = slightly disturbed slopes, HD = heavily disturbed slopes). The diagonal dashed line indicates a condition where species' mean elevation is the same for the original survey and the resurvey. Only species that occurred at least three times in both surveys were included in the model (all species-specific shifts in Appendix S8)