Arctic shrubification mediates the impacts of warming climate on changes to tundra vegetation

Climate change has been observed to expand distributions of woody plants in many areas of arctic and alpine environments—a phenomenon called shrubification. New spatial arrangements of shrubs cause further changes in vegetation via changing dynamics of biotic interactions. However, the mediating influence of shrubification is rarely acknowledged in predictions of tundra vegetation change. Here, we examine possible warming-induced landscape-level vegetation changes in a high-latitude environment using species distribution modelling (SDM), specifically concentrating on the impacts of shrubification on ambient vegetation. First, we produced estimates of current shrub and tree cover and forecasts of their expansion under climate change scenarios to be incorporated to SDMs of 116 vascular plants. Second, the predictions of vegetation change based on the models including only abiotic predictors and the models including abiotic, shrub and tree predictors were compared in a representative test area. Based on our model predictions, abundance of woody plants will expand, thus decreasing predicted species richness, amplifying species turnover and increasing the local extinction risk for ambient vegetation. However, the spatial variation demonstrated in our predictions highlights that tundra vegetation can be expected to show a wide variety of different responses to the combined effects of warming and shrubification, depending on the original plant species pool and environmental conditions. We conclude that realistic forecasts of the future require acknowledging the role of shrubification in warming-induced tundra vegetation change.


Introduction
Climate change alters vegetation by affecting species distributions and composition [1][2][3]. Additionally, the new spatial arrangement of species will further affect surrounding ambient vegetation via biotic interactions [4][5][6][7]. The mediating role of biotic interactions might be of high importance especially in arctic and alpine areas where the expansion of woody plants has been recorded and is predicted to continue [8][9][10][11][12]. This increase in shrub and tree abundance is also known as shrubification and shrub expansion, all referring to increased biomass, cover, height and/or volume of woody plants [13]. Shrubification could intensify both competition [e.g. 14,15] and facilitation between plant species [16,17]. Consequently, the changing dynamics of biotic interactions may amplify, reduce or even reverse the effects of a changing climate on biodiversity [15,[18][19][20]. However, the mediating effects of shrubification have been mainly studied in regard to their climate feedbacks, e.g. effects on atmospheric and soil CO 2 -content [21,22], albedo and evapotranspiration [10,23], and snow and permafrost [24,25].
Estimating the influence of shrubification is complicated by varied manifestation of biotic interactions between species [26,27] and across landscapes [28,29]. For example, at the transition border between boreal and arctic biomes, smaller arctic plants might be more sensitive to increases in shrub abundance that increases competition for light [14]. Varying relationships between shrubs and other plant functional groups (e.g. graminoids, forbs) due to the different ecological strategies have also been observed [14,30,31]. The spatially varying outcomes of biotic interactions are based on a hypothesis that the high productivity under benign conditions promotes competition between plant species, whereas facilitative impacts occur under environmentally harsh conditions [e.g. 29, [32][33][34][35][36][37]. Consequently, instead of the impacts following the isotherms, shrubification could cause spatially fluctuating changes in structure and composition of vegetation depending on the original species pool and environmental conditions.
To estimate these potentially multifaceted manifestations of the impacts of shrubification, species distribution modelling [SDM; 38,39], which can simultaneously cover multiple species and a spectrum of environmental conditions as well as enable predictions in space and time, is a convenient method [10,28]. However, it must be acknowledged that it is a correlative method, where the outcomes are based on statistical relationships among used variables and an assumption that the species are in equilibrium with their current environment. Therefore, interpretations of SDMs should be reinforced by comparisons with results of experimental studies and ecological theory [40].
Here, we examine potential changes in tundra vegetation as a result of warming, and the role of shrubification therein. This is done by comparing spatial predictions of vegetation changes in an oroarctic test area, located just above the current ecotone between boreal and tundra biomes. The predictions are based on two sets of SDMs: one disregarding and the second acknowledging shrubification. The used data comprise 2292 study plots (1 m 2 ) covering major environmental gradients, and thus, provide a basis for fineresolution, yet exhaustive, spatial examination of the impacts of warming and shrubification. Specifically, we firstly modelled and predicted the distributions of 130 vascular plant species, incorporating only abiotic predictors. Predictions were made under current and two future climatic conditions to assess the influence of warming alone on vegetation. Secondly, we produced spatial estimates of shrub and tree cover under current and future climatic conditions. Thirdly, the estimates of shrub and tree cover were incorporated to the abiotic SDMs of non-shrub species (n=116), followed with predictions of species distributions under current and future conditions. Finally, the vegetation changes based on the predictions in the test area were compared to assess the effects of shrubification. Parts of analyses were also repeated separately for different species groups: arctic and boreal species, and different functional groups.

Material and methods
Study area and data The data were gathered from an ecotone between oroarctic tundra and boreal mountain birch (Betula pubescens ssp. czerepanovii) forests in Finland and Norway (c. 69°N, 21°E; figure 1(a)) [29,41]. Identities and cover of all vascular plants were recorded from 1 m 2 plots (n=2292; figure 1(b)), and six abiotic variables concerning growing season, over-wintering conditions, moisture, solar radiation and soil nutrients were also derived for each plot [for rationale, see 42]. The plots are organised into groups of four so that each group forms a round study site with a radius of five meters. The sites were positioned so that they cover a wide spectrum of environmental conditions, see table S1 in appendix.
The climatic predictors in the dataset, growing degree days (GDD;°C, sum of the daily temperatures when air temperature >3°C), temperature of the coldest quarter (TCQ;°C Dec-Feb) and water balance (WAB; ratio of precipitation to evaporation), were derived from 50 m resolution data produced following [43]. In brief, the climatic data layers are based on weather station (n=942) measurements modelled and predicted using variables representing topography (solar radiation, topographic position index) and land-surface (proximity to sea, lake cover). Other abiotic factors in the data were topographic wetness index (TWI; [44] and solar radiation (RAD; MJ cm -2 d -1 ), derived in ArcGIS [45] utilising a 5 m-resolution digital elevation model, and soil nutrients (CALC), interpolated from the bedrock map of Finland [46] as the proportion of nutrient-rich calcareous bedrock. To create a variable to represent shrub cover, the percentage cover values of (dwarf) shrub species creating dense growths and reaching a height of at least 15 cm (resulted in 20 shrub species; see table S2 and [13]) were summed in each 1 m 2 plot. For tree cover, the tree canopy cover (%) of a site was used. The WAB variable was excluded from the subsequent analyses due to its high correlation (r s =−0.81) with GDD and low variable importance in SDMs as based on preliminary analyses.
To produce the environmental data for projections, the abiotic variables were also generated for a test area (8×9 km) located in the middle of the study area (figures 1(c) and S1). The test area was chosen so that it represents a variety of environmental conditions, yet so that the prevailing conditions are covered by the SDM training data, even in the future scenarios [47]. The variables were produced for 1 m 2 plots spaced every 25 m (i.e. the sub area was divided into 25×25 m squares, with a 1 m 2 plot in the centre of each square for which the data was produced). To generate the spatial estimates of shrub and tree covers for the test area, generalised boosted method (GBM from R package gbm; [48]), with five abiotic variables (GDD, TCQ, TWI, RAD, CALC) was used. To account for the effects of a changing climate, GDD and TCQ values were produced for the test area under a Representative Concentration Pathway scenario 4.5 (RCP 4.5) for years 2050 and 2070 [49]. Shrub and tree cover were also predicted for the test area under the two future climates.

Models and predictions
Distributions of vascular plant species with a minimum of ten presence observations in the data (resulted in 130 species) were modelled implementing five modelling algorithms (generalised linear mod-el=GLM, generalised additive models=GAM, generalised boosted method=GBM, random forest= RF and MAXENT) under BIOMOD2 platform [50]. First, the models were built with five abiotic predictors only (GDD, TCQ, TWI, RAD, CALC; hereafter abiotic models). The second set of models also included the shrub and tree covers (hereafter biotic models). Based on the five algorithms and implementing a majority vote ensemble technique [51], the distributions of vascular plant species were predicted to the test area under current and future climatic and biotic conditions, resulting in six projections (i.e. abiotic and biotic models under current and two future conditions); see figure S2.
The projections were compared between current and future conditions, and between abiotic and biotic models. For abiotic-biotic comparisons only the species not included in the shrub and tree cover estimates (116 species) were included to avoid circular reasoning. The comparisons were made to assess changes in mean vegetation height, functional group ratios, species richness (species gains and losses), species composition (species turnover as measured with Sørensen similarity coefficient between predicted species occurrence matrices) and extinction rate (proportion of species going locally extinct); figure S2. For a subset of analyses, vascular plant species were divided into separate groups: arctic (51 species) and boreal (65 species) to represent species' main geographical distribution area, and shrubs (14 species), graminoids (30 species) and forbs (63 species) to represent different plant functional types. Here, the shrub species to be modelled individually were species meeting the requirements specified above for shrub cover predictor with a minimum of 10 observations (see table S2). Mean heights of vascular plant species, based on [52], were used to assess changes in mean predicted vegetation height (i.e. averaging the heights of species predicted to occur).
All models were evaluated using four-fold cross validation and root mean squared error (RMSE), area under the curve (AUC), true skill statistic (TSS) and Cohen's kappa values [as done in 53]. Variable importance values were derived using the function provided in BIOMOD2 to compare the magnitude of influence of predictors in the models [50, 54 p 597]. Response curves produced by the models, representing the relationship between species presence and the predictor, were visually examined to check that models produced ecologically sound relationships between species and environmental variables.

Impact of warming on vegetation
The SDMs based on abiotic predictors perform well (mean AUC=0.835 across all 130 species; for other evaluation metrics, see table 1). Warming climate enlarges predicted distributions of a majority of studied species, and, thus, increases vascular plant species richness in the test area (table S3). Warming also alters species composition (as measured with turnover) and causes local extinctions: for example, in 2070, the mean extinction rate across the test area is almost 50% (table S3). Based on our predictions, mean species height (as averaged across species) increases subsequent to warming (figures 2(a)-(c)). Comparing ratios of predicted species richness of different functional groups (shrubs, graminoids, forbs) demonstrates warming-induced shrub-domination, especially at low altitudes (figures 2(d)-(l)).

Impact of warming and shrubification on vegetation
Shrub and tree cover models and predictions The models of shrub and tree covers perform moderately (RMSE=20.3 for shrub and RMSE=6.3 for tree cover), with a comparison of observed and predicted values revealing a bias to underprediction. Nevertheless, the predictions of current shrub and tree covers in the test area spatially match with the aerial photographs (National Land Survey of Finland; tiedostopalvelu.maanmittauslaitos.fi/tp/kartta), and all predictions correspond with earlier observations and estimates of woody plant abundance and expansion [3,8,10,55] . The warming expands distributions of both shrubs and trees as measured with the cover predictions ( figure 3, and table S4). For example, for 2050, the mean shrub cover is predicted to approximately double, and for 2070, the tree cover is predicted to reach a maximum cover of 37.2% at the lowest elevations, while higher elevations remain treeless.

Species distribution models and predictions based on biotic models
The biotic SDMs perform better than the abiotic SDMs as measured with significant improvements in all evaluation metrics; see table 1. Based on the mean variable importance values of biotic SDMs, the shrub cover is the second most influential predictor, while the tree cover has only minor influence in the models (table S5).
The area of predicted occurrence decreases for most studied species with additions of shrub and tree covers to the abiotic models: distributions of 70, 73 and 75 species (out of 116) decrease subsequent to incorporating shrubification when examining the predictions under current, 2050 and 2070 conditions, respectively. Consequently, including the shrub and tree covers to the SDMs results in an overall decrease in the predicted species richness (figure 4(a) and table S3). Shrubification also amplifies turnover and increases extinction risk of the 116 ambient species (figures 5(a), (b) and table S3).

Variation between species groups
Predicted vegetation changes vary between the species groups. Warming alone increases the species richness of boreal species but decreases the richness of arctic species ( figure 4(b) and table S3). For shrubs and graminoids, the effect of warming is generally positive; for forbs, it is negative (figure 4(c) and table S3). For all species groups, the biotic models predict lower mean species richness for the test area than abiotic models, yet the richness of arctic species and graminoids are more heavily decreased (figures 4(b), (c) and table S3).
The effects of warming and shrubification on the other two vegetation change parameters of different species groups follow the changes of all vascular plant species together, i.e. warming increases turnover and local extinction risk, with shrubification amplifying these impacts (see figure 5 and table S3). When comparing the arctic and boreal groups, warming causes bigger changes for arctic species, but shrubification more heavily alters the changes caused by warming for boreal species. From the other two species groups, forbs react more strongly to changes in temperature, while graminoids respond slightly more to increases in shrub cover.

Variation across landscape
The effects of warming and shrubification on vegetation show spatial variation, with the strongest influence at low and high elevations (figure 6). All spatial patterns hold in general for all species groups despite some fine-scale spatial variability (see figures S3-S8), yet the variation between predictions of the gains and losses of arctic and boreal species is noteworthy ( figure  S3). The spatial predictions indicate increases in vascular plant species richness at high elevations and decreases at low elevations subsequent to warming (figures 6(a) and (b)). Comparing the spatial predictions of plant species richness produced by abiotic and biotic models shows that shrubification can both reverse and amplify the changes generated by warming (figures 6(c)-(f)). At high elevations, the biotic models predict lower increases in species richness than the abiotic models, while at low elevations, shrubification amplifies the decrease in species richness. At midelevations, shrubification has locally heterogeneous impacts, both amplifying and decreasing the predicted warming-induced changes in species richness. Both the turnover and extinction risk caused by warming are strongly amplified by shrubification at the highest and lowest elevations of the test area (figures 6(g)-(r)).

Discussion
Our predictions demonstrate potential changes in tundra vegetation and call for acknowledging the mediating role of shrubification in forecasts of warming-induced environmental changes [5,53,56]. Primarily, a warming climate is predicted to affect species diversity and composition, yet subsequently the increase in the abundance of woody plants will further alter vegetation by modifying biotic interactions among species. The combined effects of warming and shrubification decrease the diversity of tundra vegetation, placing especially small arctic forbs at risk. Spatially, the predicted changes will be the strongest under extreme climatic conditions. All in all, our findings highlight the significant role of shrubification in driving vegetation change along with warming. In the next paragraphs, we present the predicted changes in more detail and discuss potential mechanisms behind the impacts. Yet, we remind, that due to the used methodological approach, the findings are contingent to three assumptions [57][58][59][60]: species are in equilibrium with their current abiotic and biotic environments, current species-environment relationships persist in future, and no temporal or spatial barriers are posed for universal dispersal. Our abiotic models predict changes in tundra vegetation structure, supporting warming-induced expansion of woody species [61,62]. In addition, warming increases total species richness, modifies species compositions and causes local extinctions. These predictions are all in line with earlier (experimental) studies demonstrating diversity loss [e.g. 62, 63]when recognising that here the increase in the total predicted vascular plant species richness is mostly due to the increase of shrub richness ( figure 4).
The effect of shrubification reflects, in general, competitive biotic interactions: including the shrub and tree covers in the SDMs causes a decrease in the mean species richness of ambient vegetation    compared to the predictions produced by abiotic SDMs. Presumably, shrubification intensifies competition for light, space, moisture and/or nutrients [14,63]. Decreasing species richness as a function of increasing productivity (assuming that shrubification increases total vegetation biomass) can also be postulated as a productivity-diversity-hypothesis [35,64]. The amplified turnover and extinction rate instigated by shrubification ( figure 5) could indicate, in addition to competition, indirect effects, such as rising summer microclimatic temperatures at ground level caused by increased shrub volume [13].
Our predictions also demonstrate that species groups can be expected to show a wide variety of responses to warming and shrubification depending on species' life-history characteristics. The varying outcomes between species groups might also reveal some mechanisms behind the impacts. For example, the increased richness of boreal plant species, and the decreased richness and high extinction rate of arctic plant species due to the warming may be caused by 'thermophilization', i.e. warmer climate favouring warm-adapted species [e.g. 61]. Warming also increases the richness of shrubs and graminoids at the cost of forbs [see also 65], which might reflect the different resource acquisition strategies of species [31,66]. Shrubification decreases the species richness of arctic plants and graminoids more heavily than the boreal plants and forbs, respectively. These outcomes might reflect varying competitive abilities between the species groups [e.g. 30,31]. For example, taller boreal plants are superior light-competitors compared to low-growing arctic plants. Contradictorily, the turnover and extinction rates are amplified more strongly by shrubification for boreal than arctic species, and for graminoids than forbs. This might, instead of an ecological rationale, be explained by spatial associations: the estimates of shrub cover and tree canopy are the highest at low elevations with higher ratios of boreal species and graminoids.
The predictions also demonstrate spatial variation in the effects of both warming and shrubification [see also 29,33], with effects being the strongest under extreme environmental conditions. At high elevations, the impacts of both warming and shrubification could be 'worsened' due to the low buffering ability of the low species richness [67,68] and/or species being specified to cold environments [61] and sensitive to competition [69]. At low elevations, the covers of shrubs and trees are predicted to be the highest with potential for a strong impact too. In addition, while in general the warming and shrubification cause changes of the same direction, the predictions reveal some spatial differences: at high elevations, warming increases vascular plant species richness that shrubification stabilises, and locally at midelevations shrubification increases species richness and diminishes turnover and extinctions. These outcomes could be explained by facilitative impacts: shrubification could, for example, locally sustain cool climatic conditions through shading or provide shelter and nutrients through trapping snow [17,70]. Nevertheless, the spatial variation in the predictions highlight how the original plant species pool and environmental conditions are crucial when estimating the effects of both warming and shrubification.

Concluding remarks
Based on our models and predictions, warming will severely change tundra vegetation. Shrub abundance will increase with further influences on ambient vegetation. In addition to decreasing species diversity, shrub expansion amplifies species turnover and local extinction rates. Arctic forbs and extreme environments particularly are at the risk. In general, the combined effects of warming and shrubification amplify one other, yet our results also demonstrate spatial variation in the effects of warming as, locally, shrubification can also diminish or even reverse the warming-induced changes. All in all, our findings call for acknowledging the mediating role of shrubification to enable more realistic forecasts of the future changes in tundra vegetation.