Resident vegetation modifies climate-driven elevational shift of a mountain sedge

Mountain plant species are changing their ranges in response to global warming. However, these shifts vary tremendously in rate, extent and direction. The reasons for this variation are yet poorly understood. A process potentially important for mountain plant re-distribution is a competition between colonizing species and the resident vegetation. Here, we focus on the impact of this process using the recent elevational shift of the sedge Carex humilis in the northern Italian Alps as a model system. We repeated and extended historical sampling (conducted in 1976) of the species in the study region. We used the historical distribution data and historical climatic maps to parameterize a species distribution model (SDM) and projected the potential distribution of the species under current conditions. We compared the historical and the current re-survey for the species in terms of the cover of important potential competitor species as well as in terms of the productivity of the resident vegetation indicated by the Normalized Difference Vegetation Index (NDVI). We found that Carex humilis has shifted its leading range margin upward rapidly (51.2 m per decade) but left many sites that have become climatically suitable since 1976 according to the SDM uncolonized. These suitable but uncolonized sites show significantly higher coverage of all dwarf shrub species and higher NDVI than the sites occupied by the sedge. These results suggest that resistance of the resident vegetation against colonization of migrating species can indeed play an important role in controlling the re-distribution of mountain plants under climate change.


Introduction
In temperate and boreal mountain ecosystems, organisms are shifting their ranges to higher elevations Pauli et al. 2012) and these shifts have been accelerating with accelerating rates of climate warming (Steinbauer et al. 2018). However, the extent and direction of elevational range shifts vary among individual species and populations (e.g. Lenoir et al. 2010;Chen et al. 2011). The reasons for this variation are not well understood so far. Traits related, for example, to plant architecture (Guittar et al. 2016), seed morphology (e.g. Vittoz et al. 2009;Matteodo et al. 2013) or variability in leaf traits (Henn et al. 2018) were found to correlate with species' responses to a warming climate, but differences in traits usually explain only a part, and often a small part of the differences in range shifts (Angert et al. 2011;Rumpf et al. 2018). The magnitude of climate change is related to the extent of range shifts but the scatter around this relationship is large Rumpf et al. 2019b). Interactions of climatic change with other facets of global environmental change, such as land-use change or atmospheric pollution, are certainly important in regional contexts (e.g. Crimmins et al. 2011;Angert et al. 2011), but often also explain only a low fraction of the observed variability in elevational shifts (Rumpf et al. 2018;Steinbauer et al. 2018). Pronounced microclimatic variation can certainly buffer species against climate warming Körner 2010, 2011), but if the underlying topographic complexity is related to variation in elevational range expansion rates over the recent past is poorly known.

3
Another factor little explored in this context is the resistance of the resident vegetation against plants that move up from below. Even if the competitive ability of species appears to decrease with elevation (Alexander et al. 2015), establishment odds of range expanders may nevertheless depend on species and trait composition of the communities growing above their (former) range margin (e.g. Dullinger et al. 2003Dullinger et al. , 2004. Indeed, the importance of resident vegetation for the likelihood that species new to a community may establish is well known from the field of biological invasions (e.g. Carboni et al. 2016;Conti et al. 2018). On the other hand, a recent experimental study suggested that alpine communities may generally be more open to seedling recruitment, and hence to the immigration of species, than communities of lower elevations, mainly because germination and establishment are easier in low-growing swards (Meineri et al. 2020). Nevertheless, competition has known effects on recruitment even in alpine communities (Alexander et al. 2018), and even a gradual decrease of this effect with elevation does not hence contradict the conjecture that resident mountain plant communities may substantially modify the pace of species' upward shifts under climate warming.
A way to more directly evaluate whether the resistance of resident vegetation may actually contribute to the variation in range shifts is to correlate the extent of observed shifts to properties of the invaded or uninvaded plant communities. Here, we developed such an approach by combining a re-survey of historical species distribution data with species distribution modelling and the use of remote sensing data. Our model system was the distribution of the sedge Carex humilis in a valley of the central European Alps. We, first, re-mapped the distribution of the species in this valley, based on data collected in 1976. Since then the average surface temperature in the study area has increased by 0.84 °C, equivalent to a temperature increase of 0.2 °C per decade, which is nearly double the global average rate between 1901(IPCC 2013. We documented the range shifts of C. humilis during these 42 years of rapid warming and evaluated the possible impact of climatic and biotic factors on these shifts. We, therefore, used the historical distribution of the species in combination with historical climate data to fit a species distribution model (SDM). We then applied the SDM to predict how the distribution of the sedge should have changed in response to climate warming over the past 42 years. We compared sites where predicted and observed the current presence of C. humilis as well as predicted and observed colonization between 1976 and 2018 match or diverge, respectively, in terms of the cover of four dwarf shrub species and the Normalized Difference Vegetation Index (NDVI). We particularly focused on dwarf shrubs because they represent an important growth form in the upper subalpine and lower alpine zone of the region, i.e. those elevational belts where a climate-driven range expansion of C. humilis is most likely. Moreover, their current spatial distribution in the region is much better documented than the distribution of other potential competitors among grasses, sedges or herbs. Alternatively, we used the NDVI, a standard index of vegetation cover, to relate the dynamics of C. humilis to properties of the entire recipient communities. We hypothesized that modelled and observed incidence and colonization by the sedge match the better the lower the dwarf shrub cover and the NDVI, i.e. that resistance of the resident vegetation, as indicated by these properties, has reduced the establishment odds of C. humilis when moving upward.

Study species
The sedge Carex humilis Leyss. is a perennial hemicryptophyte reaching 10-15 cm height, on average (Fischer et al. 2018). The species usually populates warm, dry, shallow and nutrient-poor sites like steep slopes of southern exposition and exposed crests with pronounced maxima of air and soil temperature. Nevertheless, C. humilis is not obligatorily xerothermic and has a broad ecological amplitude as well as a large geographical distribution from Central-and Southern Europe to Southern Russia and East-Asia. (Böttner et al. 1997;Fischer et al. 2018). In Europe, where C. humilis might be a relic species from Pleistocene steppes (Walter and Breckle 1994), the species mostly occurs in lowland and lower montane habitats, except on the southern side of the Alps, where it can also be found in the higher (sub)alpine regions. (Fischer et al. 2018). It clearly prefers grasslands and other open habitats, but may also occur in dry, gappy forests, especially pine forests (Fischer et al. 2018).

Study site
The study was conducted in the Matsch/Mazia-Valley in the Autonomous Province of Bolzano/Bozen, South Tyrol, (northern Italy), located in the central Alps (46°42′N 10°38′E). Since 2014 the catchment is part of the LTER network for ecological long-term research (LTER_EU_IT_097) and thus part of several monitoring programs and densely equipped with microclimatic stations (lter.eurac.edu).
The Matsch/Mazia-Valley is an inner-alpine dry valley with a continental climate characterized by distinct seasons 1 3 and wide day/night temperature amplitudes. The annual precipitation at 2000 m is around 700 mm (Niedrist et al. 2016), which makes the valley one of the driest areas in the Alpine region. The research area extends across the valley's orographically right slope, mainly exposed to the south and south-east, from the Saldur river at the bottom of the valley at about 1400 m a.s.l. to the uppermost crests at 2700 m a.s.l. The natural vegetation of the area mainly consists of montane and subalpine conifer forests dominated by Norway spruce (Picea abies), European larch (Larix decidua) or Stone pine (Pinus cembra). The undergrowth of these forests is dominated by various dwarf shrub species (Vaccinium vitis-idaea, Rhododendron ferrugineum and others). However, during a century-long history of livestock grazing parts of these forests have been replaced by pastures and meadows (Niedrist et al. 2016). While the use of lower-lying meadows has partly been intensified over the most recent decades, pastures and meadows near the tree-line as well as the natural grasslands above the tree-line are still used in traditional ways or have been abandoned.

Historical data
Within her master thesis in 1976, Ilse Rühmkorf collected data on the distribution of C. humilis in the study area (Rühmkorf 1977). These data consist of 304 occurrence records (presence only, cf. Figure 1) collected during systematic walks along parallel, equally spaced elevational transects and mapped onto a 1:25,000 terrain map, with the assistance of an altimeter. The aim of the study was to evaluate the upper elevation limit of the thermophilic vegetation (inter alias C. humilis) in the study area. As Mrs. Rühmkorf corroborated in personal communication, uppermost occurrence points of the sedge documented by her in 1976 were actually the uppermost points where she could find the species at that time.
In summer 2018, we re-visited all the occurrence points mapped by Mrs. Rühmkorf (called historical occurrences henceforth) and, in addition, did systematic walks along elevational isolines (Fig. 1). Contrary to the historical study we documented both absences and presences. We searched along isolines as far as the terrain allowed, with higher intensity in those broad habitat types where, according to the known ecology of the species, C. humilis might have found suitable conditions (e.g.: grasslands, sparse forests and forest clearings), and lower intensity in likely unsuitable habitats (closed forests). We also searched more intensively at and around the historical upper elevation limit to assess the species' current upper range limit (Fig. 1). When walking along isolines outside of forests, we recorded its presence or absence at points approximately 50-100 m apart (depending on terrain conditions). Presences were mapped at the place of the largest clone we found within a radius of 3 m at these points. To make sure that we did not overlook high-elevation  (1976) and the re-survey in 2018. Note that the historical survey only reported species occurrences while the re-survey documented presences and absences occurrences, we carried this mapping forward for another 200 m upslope the highest occurrence of C. humilis. In the closed forest, we haphazardly recorded small gaps that we encountered during isoline walks, i.e. we searched a radius of about 10 m around the center of such gaps for specimens. Additionally, we also searched for changes of lower range margins by continuing the historical transects of Mrs. Rühmkorf down to lower elevations. In total, we surveyed 1255 sites and observed the species at 687 sites. We recorded the geographical position of the center of each surveyed site by means of a hand-held GPS device.

Environmental data
Maps of abiotic habitat conditions of the study area were downloaded from the homepage of the Geoservice of the Autonomous Province of Bolzano/Bozen (https ://geose rvice s.buerg ernet z.bz.it/geoka talog /#!). These descriptors include topographic conditions (slope inclination and aspect derived from a digital elevational model (DEM, resolution: 2.5 m), solar radiation income (calculated for the years 2004-2013, based on a 25 m DEM), and terrain curvature (calculated using the curvature-function in ArcGIS® applied to the 25 m DEM). All variables were re-sampled to a common target resolution of 5 m.
We prepared two sets of climatic maps of the study area, one for current conditions, and one for the climate of the historical records. Calculations were based on two sets of raster data: the first one was a time series of monthly values of minimum, mean, and maximum temperature for the years 1979-2013, with a cell size of 100 m. The dataset was obtained from CHELSA (Karger et al. 2017) at a resolution of 0.5 arc min, and subsequently downscaled to a cell size 100 m as described in Dullinger et al. (2012) using ordinary kriging with elevation as covariable.
The second data set containing monthly values of the same variables for the time between 1943 and 2016 were derived from a global dataset with a spatial resolution of 0.5° provided by the Climate Research Unit of the University of East Anglia (CRU TS v. 3.23; Harris et al. 2014). Due to the coarse resolution of the earlier climatic data, the entire study area was included in one cell. Hence, temperature changes were the same between upper and lower elevations, meaning temperature difference due to elevation-dependent warming or micro-climatic variation (Scherrer and Körner 2010) could not be accounted for.
To generate climatic maps for the 35 years preceding the first and the second survey from these two datasets we, first, calculated changes in macroclimate at the coarse grid, and second, added these macroclimatic differences onto the fine grid. In more detail, we first defined the time span covered by our finer grid (1979-2013) as our reference period. For this time span, we calculated the means of all 12 (months) × 3 variables from the fine-scale grids and averaged these temporal means across space, i.e. across the entire study area (= coarse-scale resolution). Subsequently, we computed the difference (= delta) of each monthly variable between each year of the entire coarse-scale time series and the averages of the reference period. Third, we projected these deltas back onto the fine scale by adding (or subtracting) them to (from) the 1979-2013 mean of each fine-scale grid cell. In this way, we produced year-wise values of all monthly variables for the entire time series 1943-2016 at the fine resolution. Based on this time series, we calculated the bioclimatic variables "Mean Temperature of Warmest Quarter" (Abbreviation: "Bio10") and "Maximum Temperature of Warmest Month" for the periods 1943-1977 and 1982-2016 as biologically relevant climatic descriptors for the historical and current surveys of C. humilis, respectively. We used 1977 as the end of the historical period to have the same period of 35 years for both the survey and the re-survey. We did not include precipitation-related climatic descriptors as precipitation did not change significantly between original and re-survey years.

Vegetation data
We collected two different kinds of data as indicators of competitive resistance of the resident vegetation to colonization by C. humilis. The first dataset was rasterized distribution maps of the most frequent dwarf shrubs in the study area, Rhododendron ferrugineum, Vaccinium myrtillus & Vaccinium vitis-idaea and Calluna vulgaris.
These species can all form dense mats and may hence be considered major competitors for colonizing C. humilis. Moreover, their ecological profiles differ along a gradient of snow cover duration and ability to withstand dry conditions in summer (Ellenberg 2009). R. ferrugineum prefers sites that are covered by deep snow until mid-spring, at least, and often occurs in shady sites, whereas Calluna vulgaris prefers warm and dry slopes of southern exposition, often with a relatively short snow cover. The two Vaccinium species have intermediate profiles, with V. myrtillus more similar to R. ferrugineum, and V. vitis-idaea able to withstand harsher conditions at less snow-protected sites. The distribution maps of these species were obtained from the Geoservice of the Autonomous Province of South Tyrol (https ://geose rvice s.buerg ernet z.bz.it/geoka talog /#!) as shapefiles with a scale of 1:25.000. The maps have been created in 2010 based on a detailed vegetation map and data search, the participation of the local community (interviews) and analysis of orthophotographs. The data are percentage cover estimates in three classes for each of the four species (1-25%, 25-50%, > 50%). A numerical validation of map accuracy was not available from the data provider, but according to our personal experience in the Mazia Valley, the current distribution of dwarf shrubs is well described by this map. In general, the center of the mapped regional distribution of the shrubs is the subalpine belt and extends approximately from 1,300 to 2500, 1000 to 2400, 1000 to 2500 and 1600 to 2500 m asl. for Rhododendron ferrugineum, Vaccinium myrtillus, Vaccinium vitis-idaea and Calluna vulgaris, respectively. For all the following steps including Vaccinium spp. we used a newly generated shapefile that has been created from the average values of the two shapefiles of V. myrtillus and V. vitis-idae. We further converted the shapefiles to raster maps with a cell size of 5 m using the maximum area method to assign cover values to the individual cells. Subsequently, we re-classified the species' cover classes to either zero or the median value of the respective classes (0%-12.5%-37.5%-75%).
The second dataset was a map of the Normalized Difference Vegetation Index (NDVI) in the study area. This map has a spatial resolution of 10 m and was calculated based on the Sentinel-2 Multispectral Imager remote sensing images from 13 July 2018 (0% cloud cover over the study area). The NDVI is an optical indicator of total vegetation cover and/ or productivity that is based on the ratio of near-infrared to the red light which are reflected or absorbed by vegetation, respectively (e.g. Borowik et al. 2013;He et al. 2015). We resampled the map to the common target resolution of 5 m.

Species distribution model (SDM)
We used the 304 historical occurrence points of C. humilis in combination with the topographic descriptors and the climatic maps of 1943-1977 to parameterize an SDM of the species in the BIOMOD package in R (Thuiller et al. 2009). As in the historical study, only presence data were collected we additionally added 291 absences randomly chosen in areas where, according to the remembrance of Mrs. Rühmorf, C. humilis did not occur in 1976 (mostly concave parts of the terrain, or those exposed to the North, or areas above the uppermost occurrences in 1976). Prior to modelling, we screened all available predictor variables on correlations. We found that topographic variables (slope, aspect, curvature) and solar radiation were only loosely correlated among each other and with climatic descriptors while the three bioclimatic variables were highly correlated (Pearson r > 0.9). We hence decided to use, in addition to the topographic variables and solar radiation, only one of these climate variables as a predictor in the SDMs and selected the "Mean Temperature of Warmest Quarter" as the most straightforward indicator of thermal energy available during the vegetation period. The map of this climatic descriptor was resampled to the common target resolution of 5 m in ArcGIS using the maximum area method. With this approach, we ensured that all 5 m cells inside the resampled 100 m cells have the same value as the initial cell.
We fitted SDMs by means of five different algorithms: Generalized Linear Models (GLM), Generalized additive models (GAM), Generalized Boosted Regression Modeling (GBM), Random forests (RF) and Multivariate adaptive regression splines (MARS). In case of GLMs, we used polynomials of second order to model predictor-response relationships, in case of GAMs spline smoothers with four knots. We followed default options in BIOMOD for all other parameters of the modelling algorithms. We applied four evaluation runs for each algorithm using 70% of our data for testing and 30% for validation. Models were validated by means of the "True skill statistic" (TSS, Allouche et al. 2006), a measure of distinctive ability that varies between 0 and 1. To assess the importance of predictor variables in the fitted models, we used a built-in permutation function that scores the impact the respective variable has on the discrimination ability of the model on a scale ranging from 0 to 1 (Thuiller et al. 2009).
We used the models fitted with the historical distribution of the species and the historical climate to project the expected spatial distribution of the species under the current climate, i.e. to select the subset of our 1255 visited sites where we should have found the species in 2018 if it would have shifted its range in accordance with the warming climate. We used 'ensemble' predictions, i.e. we predicted occurrence probabilities as weighted averages of predictions from the five modelling algorithms, with weights defined as mean TSS-scores of the four validation runs per model. According to the default in BIOMOD, only models with a TSS > 0.6 were allowed to contribute to these ensemble predictions. We converted the resulting probabilistic predictions to 0/1, i.e. presence-absence predictions, using the probability threshold that maximized the TSS-score.

Statistical analyses
We assessed the shift of upper and lower range margins of the species. We, therefore, selected the 30 uppermost and 30 lowest occurrences of the species in both the survey and re-survey data and compared their elevations by means of two-sample t-tests. We calculated the "elevational-shift-perdecade" by dividing the difference of the present and the historical uppermost and lowermost elevational limit of the species by the number of the decades since the historical study (4.2 decades for the present study).
To evaluate whether the four dwarf shrub species as well as the NDVI may represent meaningful indicators of the resistance of the resident vegetation on C. humilis we used two different approaches. First, we focused on all sites predicted to be occupied by the SDM in 2018 (n = 659). For this subset we compared sites where we either found or did not find the species in 2018 in terms of the mean cover of Vaccinium, Rhododendron and Calluna as well as in terms of the NDVI, using two-sample t-tests. We expected that sites where we found C. humilis in 2018 should have the lower cover of these species, on average, if competition among them and the sedge is actually a process affecting the distribution of the latter. Dwarf shrub cover and NDVI were moderately correlated at the 1255 recorded sites (Kendall's τ = 0.42, 0.3, and 0.43 for Calluna, Rhododendron, and Vaccinium, respectively, p < 0.001 in every case).
Second, we more narrowly focused on predicted colonization events between 1976 and 2018, i.e. on the 490 sites for which the SDM projected absence in 1976, but the presence in 2018 (excluding all sites where the species had already been mapped in 1976). For these sites, we assessed whether the mapped presence of C. humilis in 2018 is related to the cover of the four indicator species and to the NDVI by means of a multiple logistic regression model (using the glm function in R with a binomial error distribution and a logistic link function). In addition, we did the complementary analysis, i.e. the same regression analysis, but focused on the sites for which no colonization is predicted by the model, i.e. where the model projected absence of the species in both 1976 and 2018 (again excluding all sites where the species had already been mapped in 1976).

Climatic trends
All 1255 mapped sites of the re-survey are located in one single cell of the climatic dataset we used to calculate the climatic conditions for the historical and recent study. The mean annual temperature of this cell increased by 0.84 °C, the rate of temperature increase was hence 0.2 °C per decade. At the same time, the maximum temperature of the warmest month increased by 1.2 °C, and the mean temperature of the warmest quarter by 1.13 °C, which means summer temperatures warmed more pronouncedly than the annual mean temperature.

Species distribution model
The SDM of C. humilis fitted with the historical distribution and climatic data had the excellent predictive ability (TSS = 0.846, Sensitivity: 89.769, Specificity: 94.845, cf. Table 1). In all single models, and hence also in the ensemble model, the mean temperature of the warmest quarter (Bio10) was by far the variable with the highest impact on the species' distribution within the study area, followed by solar radiation.
For the 973 sites that had not been suitable to the species in 1976 according to the SDM, i.e. all the surveyed sites that are candidates for a predicted colonization event, the model projected observed presence and absence of the species in 2018 correctly in 308 and 329 cases, respectively. Vice versa, the model predicted 182 false presences and 154 false absences (cf. Fig. 2).

Effects of resident vegetation
Of the total set of 1255 surveyed sites, the SDM predicted 659 to be suitable to the species in 2018. Among them, those 453 sites where we found C. humilis differed significantly from the 206 where the species was absent in terms of cover of the four dwarf shrub species and of the NDVI. Coverage of Rhododendron, Vaccinium and Calluna was on average higher by 12%, 14%, and 5%, respectively, on sites without C. humilis, and the NDVI was higher by ca. 10% (cf. Fig. 3, Table 2).

Multiple regression models
The multiple regression model focused on the 490 sites predicted to have become colonized between 1976 and 2018 demonstrated that the likelihood of such colonization depended on the NDVI (cf. LR1 in Table 3). More precisely, this likelihood slightly increased from very low to intermediate NDVI-values, and then strongly decreased with a further increase of the index, but it did not depend on the cover of the dwarf shrub species. The apparent contradiction to the results of the t-tests may arise from correlations between the cover of dwarf shrubs and the NDVI. Indeed, when evaluated in univariate models, the cover of each dwarf shrub species was negatively correlated with the occurrence of C. humilis (cf. Table 4 in Appendix). For those 483 sites which, according to the SDM, were unsuitable to the species both in 1976 and 2018, i.e. where colonization by the sedge was not predicted, the relationship between the likelihood of (nevertheless) finding C. humilis in 2018 and the NDVI was the same as in the case of the sites predicted to have become colonized (see LR2, Table 3). However, among the dwarf shrubs Rhododendron was negatively associated with the incidence of C. humilis while the presence of Vaccinium was positively correlated to the occurrence of the sedge in this multiple regression. Fig. 2 The 973 sites mapped in 2018 that had not been already colonized by Carex humilis in 1976, subdivided into two major groups (expected (dots) / not expected (triangles) to have become colonized according to a species distribution model). Each group was further divided into two subgroups (species observed / not observed in 2018). Green dots and black triangles represent sites where the species was observed, yellow dots and red triangles those where the species was not observed in 2018 Fig. 3 Average cover of dwarf shrub species at sites predicted to be suitable for Carex humilis by a species distribution model and found to be either occupied or unoccupied by the species in 2018 1 3 In univariate models, none of the dwarf shrubs was significantly correlated with the occurrence of C. humilis (cf. Table 4 in Appendix).

Shifts of range margins
The results of this re-visitation study confirm many other observations of upslope shifting elevational range limits of plant species in the European Alps and other mountain ranges, especially in the subalpine and alpine belts (e.g. Rumpf et al. 2018Rumpf et al. , 2019b. According to our data, the leading edge of the species moved upward at a rate of more than 50 m per decade. This is exceeding the velocity reported in the last global meta-analysis (12.2 m per decade; Chen et al. 2011) by a factor of four, and it is more than twice as fast as the average rate reported for plants of the European Alps by Rumpf et al. (2018), although within the variation reported in this study. However, it is comparable with the maximum migration velocity of 53.2 m per year reported by Walther et al. (2005) in his study conducted in the nearby south-eastern Swiss Alps. In addition, the observed 215 m upslope movement of the species is nearly twice of what one would expect based on the calculated temperature increase of 0.84 °C, at least when assuming a lapse rate of 0.65 °C per 100 m of elevation, which is the usual lapse rate of summer temperatures in the European Alps (Rubel et al. 2017). Several reasons could have contributed to this particularly fast shift of the upper range margin. First, species of subalpine and alpine areas are, outside of forests, largely decoupled from air temperatures and rather depend on the microclimate near the ground which can vary substantially over short geographical distances (Scherrer and Körner 2011). The effects of changing air temperatures on microclimates are poorly understood, but on dry, rocky and / or exposed sites, i.e. the preferred habitats of C. humilis, the increase in micro-temperature regimes might even exceed the increase of the air temperature measured 2 m above the ground due to higher radiation sums (Aalto et al. 2018). Moreover, these habitats often have shallow soils and are, in the study area, predominantly exposed towards southern directions. As a consequence, increasing top-soil temperatures will also result in drier conditions in summer, which may additionally favor drought-resistant species such as C. humilis and facilitate their colonization. Third, rising temperatures are associated with earlier snowmelt which benefits early flowering species such as C. humilis (Lenoir et al. 2008). Finally, livestock grazing has also extended to higher elevations with warming and this may have facilitated seed dispersal of the species to high lying sites and hence accelerated upward movement (Svenning et al. 2014).
In contrast to upper range margins, rear edges of C. humilis in the study area did not retreat to higher elevations. In fact, the species occurs at many warmer and lower sites in other parts of the European Alps (and outside of the Alps), i.e. it likely did not reach its warm limits within the study area in the 1970s and probably does not reach it under current conditions. Loss of the species from some sites near the valley bottom has hence probably not been triggered by climate warming but rather resulted from changing land use, and especially from intensified grazing, partly in combination with fertilization and irrigation.

Effects of resident vegetation on range shift
Despite the rapid upward movement of its leading edge, we could not find the species at all sites where it was predicted to be present in 2018 by our species distribution model, and neither at all those sites that became suitable between 1976 and 2018, i.e. where we would have expected a colonization event. In part, these mismatches are likely due to, first, inherent errors of the underlying SDM that may, for example, arise from a non-equilibrium distribution of the species during the historical survey (Guisan and Thuiller 2005;Early and Sax 2014), and, second, to dispersal limitation (Alexander et al. 2018;Rumpf et al. 2019a). However, the cover of potential competitor species as well as the total vegetation cover, as indicated by the NDVI, were significantly correlated to whether the species occurred at sites predicted to be suitable to it or not. We note that the distribution of dwarf shrubs and the NDVI will likely have changed since C. humilis has actually colonized most new sites over the past 42 years. However, past and current distribution of the shrubs and the NDVI certainly overlap to a considerable degree. Moreover, it is very unlikely that the correlation found developed post-hoc just by chance. We hence conclude that successful colonization by C. humilis is co-determined by site features that are related to these variables. Such a correlation could essentially result from two sources: the vegetation properties may either indicate some abiotic conditions that are important for the distribution of the sedge but not represented by the predictors in the SDM, such as e.g. nutrient or water availability, or length of the snow-free season; or they indicate the type and strength of biotic interactions that may restrict or facilitate colonization by C. humilis. We argue that interactions at least contribute to the mismatches observed because, first, the DEM underlying several predictor variables in the SDM has a relatively fine resolution and snowmelt times as well as water and nutrient re-distribution in alpine landscapes are normally closely related to fine-scale topography (e.g. Körner 2003). Second, the dwarf shrubs included as potential competitors have distinct ecological niches, but all of them were negatively related to the incidence of C. humilis. If the shrubs would indicate unrepresented abiotic conditions, correlations with C. humilis should vary among them, with positive correlations for species that have ecological requirements similar to the sedge, such as Calluna vulgaris, and negative correlations with species that have contrasting requirements such as Rhododendron ferrugineum. While it is true that differences in dwarf shrub cover of colonized and uncolonized sites are only moderate in magnitude, on average in our (coarsely measured) data, even such moderate differences can effectively increase the resistance of the resident vegetation, especially if micro-scale conditions restrict the incidence of both dwarf shrubs and C. humilis to only part of the 5 × 5 m cell Körner 2010, 2011). Related to that, the subset of sites where we found C. humilis and those without the species did not differ in elevation. It is hence also unlikely that differences in elevational niches between the dwarf shrubs and C. humilis were responsible for the pattern of mutual exclusion we found. Finally, in the multiple regression models the effect of the NDVI proved most important, independent of which species were dominating the vegetation, and hence also independent of the variation in ecological conditions that are associated with different dominant species in alpine environments (Ellenberg 2009). This effect was actually curvilinear, suggesting that a dense resident vegetation hindered the establishment of C. humilis, but that the species also had problems to colonize extremely sparsely vegetated sites. Both abiotic and biotic factors may contribute to the latter observation: sparsely vegetated sites often are rock faces or scree where a lack of soil or frequent substrate movement restrict plant growth. On the other hand, this can also be extremely exposed sites where neighbors may improve microclimatic conditions and hence play an important role as facilitators, even for stress-tolerant species such as C. humilis (Choler et al. 2001;Callaway et al. 2002;Cavieres et al. 2013).
The mechanisms by which the cover of dwarf shrubs, or of other dominant competitors indicated by higher NDVI values, affect the performance of C. humilis cannot be deduced from coarse data on incidence and abundance such as those we have used in this study. However, significant effects of the NDVI and the shrubs in the model restricted to sites that were predicted to have become colonized between 1976 and 2018 suggest an impact on vital rates involved in population establishments, such as germination or seedling and sapling survival. In case of dwarf shrubs, their effects on soil properties may matter in this respect. According to Landolt et al. (2010), C. humilis is a species that prefers base-rich soils with an intermediate to high pH and an intermediate humus layer. However, the litter of dwarf shrub species tends to decompose slowly under subalpine and alpine climatic conditions, and the accumulating litter contributes to humus accumulation and lowering of the pH (Ellenberg 2009). It may hence be that dwarf shrubs decrease establishment odds of C. humilis not, or not only, by competition for resources, such as light, but also by changing soil conditions to the disadvantage of the sedge. This probably affects sensitive life history stages such as germination and seedling establishment more strongly than the survival of adult tussocks.
Negative effects of dwarf shrubs on sedges and other herbaceous species were already shown in experimental studies. For example, the cover of Vaccinium myrtillus was found to change the composition and reduce the diversity of associated herbaceous plants in the Swiss Alps (Boscutti et al. 2018); elimination of dwarf shrubs also stimulated leaf production of Carex vaginata and generally increased recruitment of herbaceous species in the Scandes (Klanderud 2005(Klanderud , 2010. Furthermore, Carex bigelowii was associated negatively with dwarf shrubs in protected sites of alpine Norway (Olofsson 2004). As C. humilis is known as a rather weak competitor (Hroudova and Prach 1986) similar competitive restriction of our study species by dwarf shrubs is plausible. However, under abiotically marginal conditions, the neighborhood of dwarf shrubs was also found to increase the likelihood of finding species such as C. bigelowii (Olofsson 2004). This might be one reason why we found a positive effect of Vacccinum spp. in the multiple regression model when considering sites that were predicted to have not become colonized by the sedge. By definition, unsuitable sites are those where the species is at the margin or beyond its (realized) niche space and where dependence on facilitative neighbors should hence increase (Liancourt et al. 2005;Michalet et al. 2006). Indeed, Vaccinium includes the stress-tolerant V. vitis-idaea, a species that can also be found on exposed ridges and which could hence actually serve as a facilitator under such conditions. By contrast, R. ferrugineum is a species that often occurs in shady sites and requires a long-lasting snow cover. It is hence much less tolerant to harsh conditions than C. humilis itself and unlikely to protect the sedge against such conditions anywhere in the alpine landscape. Its negative association with the sedge on the sites unsuitable to C. humilis thus rather results from (stronger) sensitivity of the sedge against competition at sites where these two species co-occur. Alternatively, the cover of Rhododendron may distinguish those sites marginally suitable from those definitely beyond the niche of C. humilis.

Conclusions
While the 'average subalpine and alpine plant species' moved upward in elevation during the recent decades, the huge variation in range dynamics across species, as well as within species across different sites, remains only partly understood. Here, we show that interactions with the resident vegetation most likely play an important role in these dynamics. While this is in line with theoretical expectations, the resistance of resident vegetation has so far played a surprisingly marginal role in analyses of observed past or in modelling of expected future range shifts of mountain plants. This neglect is likely due to difficulties in spotting biotic interactions in re-survey studies  as well as in including biotic interactions into prevalent modelling tools, especially in a species-specific way (e.g. Dormann et al. 2018). However, competition and facilitation among plants are likely diffuse in most cases and coarse community-level indicators may hence often be sufficient to estimate resistance against new colonizers. This is particularly likely for the establishment of new species from seed that can effectively be constricted even by species that are competitively inferior in the adult stage (e.g. Dullinger et al. 2003;Alexander et al. 2018). Including non-species-specific proxies such as the NDVI or other optical indices may thus have considerable potential to improve both analyses and predictive models of mountain species' range dynamics under climate change.
provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.