Introduction

High mountain environments generally reach higher degrees of naturalness and are exposed to fewer direct human impacts than lowlands (Rodríguez-Rodríguez and Bomhard 2012), show high habitat variation over short distances (Scherrer and Körner 2010), but are often strongly fragmented (Antonelli et al. 2018). A number of studies in European mountain ranges found rapid directional shifts in plant species compositions of alpine communities and increasing species numbers in synchrony with rising temperatures (e.g. Gottfried et al. 2012; Pauli et al. 2012; Steinbauer et al. 2018), thus corroborating the expectation of a high sensitivity of mountain floras to climatic changes (Lenoir and Svenning 2015). These local increases in species richness are the result of an upward shift of species’ distribution ranges (McCain and Colwell 2011), where their lower range margins were reported to move at least at similar velocities as the upper ones in the Eastern Alps (Rumpf et al. 2018). Of particular concern is the situation in small mountain ranges with a limited vertical extent, where suitable cold-enough habitats are expected to rapidly become scarce in the usually cone-shaped mountain topographies. This applies to most of the Mediterranean mountain ranges (Jiménez-Alfaro et al. 2014). Owing to the orographic isolation of their alpine refugia, which were even ice-free during Pleistocene glacial periods (Flantua et al. 2020; Gomez-Ortiz et al. 2013), Mediterranean mountains host an exceptionally high number of cold-adapted endemic plants (Fernández Calzado et al. 2012; Kazakis et al. 2007; Noroozi et al. 2018; Stanisci et al. 2005), and, thus, they are exposed to a high risk of climate change-induced biodiversity loss on the species level. Although observational studies in high-mountain vegetation in the context of climate change are still underrepresented in semi-arid areas (Giménez‐Benavides et al. 2018), upward shifts of species including drought-tolerant shrub species were found in Mediterranean mountains by re-visiting historical survey sites (Evangelista et al. 2016).

A prominent example of a small-scaled isolated Mediterranean high-mountain range is the Sierra Nevada of southern Spain, exhibiting the only true alpine region between the North African Atlas mountains, the Sistema Central, and the Pyrenees, all being several hundreds of kilometers away. The Sierra Nevada marks the southernmost limit of the influence of the Quaternary glaciations in Europe, when it was covered with glaciers only in areas above 2500 m a.s.l., while large areas remained free of permanent ice even at higher elevations (Gomez-Ortiz et al. 2013). This has promoted speciation, resulting in the highest level of vascular plant endemism on the Iberian Peninsula (Buira et al. 2020), and probably in Europe (Lorite et al. 2020). For the Sierra Nevada, a mean temperature increase of 4.8 °C at the end of the twenty-first century was simulated, indicating a vertical shift of the suitable habitats for Sierra Nevada key species at a rate of 12 m/year (Benito et al. 2011). If this holds true, extinction events may occur within a few decades, because many of the endemic species are concentrated in the uppermost bioclimatic belts. A pan-European study (including Sierra Nevada) detected widespread thermophilisation of alpine vegetation, i.e. a plant community transformation towards more thermophilous species assemblages, after a period of only seven years (2001–2008; Gottfried et al. 2012). Contrary to temperate and boreal high mountains, a decline in the overall species number per summit was found for all assessed Mediterranean sites (Pauli et al. 2012), suggesting increasing water stress as explanation. For the simple reason that plants are not only constrained by energy input but also by water supply, temperature changes alone will not be sufficient to understand biodiversity changes (Crimmins et al. 2011), especially not in water-limited environments (Cowles et al. 2018). For example, McCain and Colwell (2011) predicted a tenfold increase in extinction risk for vertebrates in mountain areas, when considering combined effects of warming and decrease in precipitation, which also may be relevant for plants as sessile organisms. For the Mediterranean high mountains, a general warming trend coupled with a reduction in the average annual rainfall is projected, resulting in longer periods with drought stress (Nogués Bravo et al. 2008; Sillmann et al. 2013). As a consequence, a higher vulnerability of alpine Mediterranean vegetation is expected through declines in snow cover duration (Beniston 2003), and shifts from snowfall to rainfall (Pérez-Palazón et al. 2018). To assess relationships between plant diversity and the two most vital aspects of climatic changes, we used permanent plot data sets of plant assemblages on Sierra Nevada summits of the international GLORIA monitoring network (Global Observation Research in Alpine Environments; www.gloria.ac.at) with three surveys spanning a total period of 18 years (2001–2019). Combined with continuous measurements of soil temperature at the plot locations and large-scaled temperature and precipitation data, we address the following questions:

  • Do the previously found patterns (Fernández Calzado and Molero 2013; Pauli et al. 2012) of change in species richness, species gains and losses, and species’ abundances constitute ongoing trends?

  • Have the observed changes been driven by particular growth forms or vertical species distribution types?—i.e., can processes of ‘thermophilisation’, ‘shrubification’ and/or loss of endemic species be observed?

  • Can the observed changes be attributed to either rising temperatures, to water-related factors (precipitation, snow cover duration), or both?

Methods

Site description

The Sierra Nevada is a European and Mediterranean biodiversity hotspot (Blanca et al. 1998; Cañadas et al. 2014) located in southern Spain, extending from 36° 50′ 24″ to 37° 15′ 0″ N in latitude and 3° 44′ 24″ to 2° 35′ 24″ W in longitude. Within its limited surface area of 2100 km2, it has a complex orography with a wide altitudinal range from 200 to 3482 m a.s.l. (Fig. 1a). The climate is Mediterranean, characterized by cold and wet winters and hot and dry summers (pronounced summer drought in July–August; Gómez 2002). Average temperatures are below 0 °C during winter with a snow cover that can persist up to 8 months in the highest areas (occasionally up to 10 months in small patches of snowbeds). The average annual rainfall is highly irregular, with values ranging between 400 and well over 1000 mm in areas above 2500 m (Polo et al. 2019), depending mostly on elevation and the topographic position. Precipitation is mainly in the form of snow above 2500 m a.s.l. The core area is composed of siliceous rocks, mainly mica-schists, surrounded by limestone and dolomite at lower elevations (Jabaloy et al. 2008).

Fig. 1
figure 1

Location of the GLORIA summits in Sierra Nevada, Spain, and the climatic conditions during the recent decades. a Geographic location of eight GLORIA summits of the two study regions Sierra Nevada—West (SNE) and Sierra Nevada—Northeast (SNN). Map source: modified after Esri et al. (2019). b Days with snow cover, i.e., days with a maximum of 0.5 °C and a minimum of − 5 °C and a maximum daily difference of 2 °C (after Liberati et al. 2019; Teubner et al. 2015), derived from soil temperature measurements in each cardinal direction on each summit (dots and triangles), as well as the loess local weighted regression (lines). c Mean annual air temperature and d annual precipitation of the Sierra Nevada, derived from ERA5 data. The grey line represents the 20-year Gaussian low-pass filter smooth trend. On the x-axis, vegetation survey years in the two study regions are shown (dots and triangles)

Site design

The two GLORIA study regions of the Sierra Nevada (SNE: Sierra Nevada—West, SNN: Sierra Nevada—Northeast) include eight summits, arranged along an elevation gradient between 2668 and 3327 m (Fig. 1a, Table 1). On each summit, eight summit area sections (SAS), and four 3 m × 3 m quadrat clusters (in each cardinal direction) were established following the standardized protocol of the GLORIA Multi-Summit approach (Pauli et al. 2015). The presence of vascular plant species was recorded for each SAS and the four corner quadrats of each quadrat cluster (i.e., per summit 16 quadrats of 1 × 1 m). The percentage cover of each vascular plant species was estimated visually only in the quadrats. The summits of SNE were surveyed in 2001, 2008 and 2015, those of SNN in 2004, 2011 and 2019, except summit MIR, where the baseline survey was in 2006.

Table 1 Number of vascular plant species per summit and survey in the Sierra Nevada, Spain

Climate data

In the centre of each quadrat cluster, a data logger (Onset TidbiT and later GeoPrecision M-Log5W) was buried at 10 cm soil depth, measuring soil temperature at hourly intervals since the setup of the permanent plots. Measurements 24 h before and after each logger change were discarded. All temperature values outside the 0.05–99.95 percentile, temperature differences between two subsequent measurements exceeding a threshold of 4 °C and plotted time series were checked for plausibility. Obvious inconsistencies, most likely deriving from logger failure, were excluded. Gaps were imputed based on the measurements of the remaining loggers on the same summit using the function amelia (R-package Amelia II; Honaker et al. 2011) which applies an EM (expectation–maximisation)—algorithm on multiple bootstrapped samples of the incomplete original dataset. Imputation of missing values was repeated 30 times. Soil temperature data were used to calculate days with snow cover, i.e., days with a maximum of 0.5 °C and a minimum of − 5 °C and a maximum daily difference of 2 °C (after Teubner et al. 2015).

For longer time series of climate-related indices, we used ERA5 data which provide hourly estimates of air temperature and precipitation on a 30 km grid at different elevation levels (from 1979 to 2019 for the coordinates 37° 2′ 44.988″ N, 3° 21′ 7.9992″ W; https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5). We used temperature at 700 hPa which is equivalent to an elevation of 3000 m a.s.l. ERA5 combines available observation data from satellites and in situ stations, which are assimilated and processed using ECMWF's Integrated Forecast System (IFS) Cycle 41r2. Despite the low spatial resolution of ERA5, the time series patterns are consistent with those of local soil temperatures (Fig. S1).

Analyses were conducted with prior periods of one, five and seven years (i.e., means of mean annual air temperature, annual precipitation and mean number of days with snow cover per year over a given period prior to the vegetation record). For the concept of prior periods, see Gottfried et al. (2011) and Gottfried et al. (2012). Yearly means were calculated from June to May of the following year to incorporate climate factors until the onset of a growing season and, thus, the period immediately preceding a vegetation survey, as well as the entire winter season (e.g., for mean precipitation of a period of five years prior to survey year 2001 (Prec01p(5)) = sums from June 1996 to May 2001 divided by five). The model with the best performance was selected using the Akaike Information Criteria (AIC; Burnham and Anderson 2002).

Vegetation data

All thirteen annual species were removed from the dataset (Table S1, Table 1) because of high inter-annual fluctuation, which can influence the results in a random way. Thirty records identified only to the genus level (3 in SNE and 27 in SNN) were also removed from the data, however, these cases only represented poorly developed individuals with low cover (median = 0.1%) which may well belong to a determined species.

All statistical analyses were performed using R (v. 3.6.1; R Core Team 2019). Generalized linear mixed models were derived using template model builder (R-package glmmTMB; Brooks et al. 2017). Table S2 gives an overview of the statistical models employed, including the response variable, fixed effects, random intercept terms and distribution family.

For the response variables, the appropriate distribution was determined by building a histogram (package lattice; Deepayan 2008) and performing likelihood ratio tests (anova function from stats package; R Core Team 2019) between models using different distribution families. To avoid multicollinearity between predictors, the Variance Inflation Factor (VIF) of the function corvif (Zuur et al. 2009) was tested (VIF > 2: concerned predictor was removed). Model assumptions were checked visually by plotting fitted values and residuals of model data (plot function from base; R Core Team 2019) and for small sample sizes additionally with simulated residuals (function simulateResiduals from package DHARMa; Hartig 2019). Correlations between raw data and predicted model data were tested with the function simulate (package glmmTMB; Brooks et al. 2017) and visualised with ggplot. If required, overdispersion was tested with the functions check_overdispersion (package performance; Lüdecke et al. 2020) and testDispersion (package DHARMa; Hartig 2019). Significance of the predictor variables’ effects of all models was extracted with the function summary (R Core Team 2019) and for pairwise comparisons the function emmeans was used (package emmeans; Lenth 2019). All graphs were generated with ggplot (package tidyverse; Wickham 2017).

Species richness was calculated as the number of species per SAS and survey year. Colonisation events were defined as the number of species per SAS present at the time of the resurvey, which were absent in the respective SAS at the previous survey. Disappearance events were defined as the number of species per SAS absent at the time of the resurvey, which were present at the previous survey. Colonisation and disappearance events per species were treated as Bernoulli trials from a local species pool, i.e., all species ever found on a given summit (1: successful colonisation of a species into a not-occupied SAS or disappearance from an occupied SAS; 0: species not colonised SAS or not disappeared from a SAS).

Species cover sum per quadrat was calculated as the cumulative cover of all species present in a quadrat. Cover per species and quadrat was used for species-specific analyses.

The thermic indicator was calculated after Gottfried et al. (2012) for each quadrat and survey as an averaged composite score of the species-specific altitudinal distribution values (altitudinal ranks (ARs): a proxy for the thermal preferences of species and their distribution along the elevation gradients after Gottfried et al. (2012) and expert judgment of authors for new species), weighted by the cover of the respective species.

Species characteristics were used according to expert judgment of authors: growth form (gr: graminoid; he: forb-hemicryptophyte; cu: cushion; sh: shrub (including dwarf-shrubs and matorral species); su: succulent), and according to Lorite et al. (2020): endemic restricted to Sierra Nevada (yes/no), spinescence (yes/no) and hygrophilous (yes/no). For the change of the proportional cover sums per species group (shrubs, spiny, endemic), the proportional cover sums per quadrat c within an interval of zero and one were compressed to avoid zeros and ones after Smithson and Verkuilen (2006) by c′ = \(\frac{c\times \left(n-1\right)+0.5}{n},\) where n is the number of observations.

To investigate the effects of climatic drivers, the average vegetation change per year (i.e., divided by the number of years between the respective surveys) was related to the change of climatic factors in the prior periods (e.g., for change of precipitation from 2001 to 2008 (ΔPrec08p(5)) = Prec08p(5) − Prec01p(5)).

Results

Climate data

Snow cover duration on study summits changed considerable, most pronounced on the two higher summits of the SNE with a strong decrease between 2001 and 2006 followed by an increase until around 2010 and again decrease thereafter (Fig. 1b). The highest summit of SNN fitted into this pattern, and showed an increase again after 2013. All lower summits showed less pronounced changes in snow cover duration. Annual precipitation varied between ~ 350 and 1100 mm over the past forty years, with a peak in 1996 followed by years with lower values, and further peaks in 2010 and 2018 (Fig. 1d). Mean annual air temperature ranged from ~ 1.2 to ~ 3.25 °C (Fig. 1c). The years prior to the baseline survey showed stronger inter-annual fluctuation than in the following years, but a maximum in 2018.

Vegetation changes with time and space

The total number of species in Sierra Nevada—West (SNE) dropped from 76 in 2001 to 71 species in 2008 and increased to 75 in 2015. In Sierra Nevada—Northeast (SNN), the total number of species increased from 67 species in the first survey (2004, 2006) to 73 in 2011 and to 77 in 2019. For numbers of species per summit and survey, see Table 1; for species’ presence/absence over all SASs per study region, see Table S1. In total, 13 woody plant species occurred on the summits (shrubs in Table S1). Neophytes did not occur on nor colonised the summits.

In the study region SNE, mean species richness per SAS decreased significantly on three of the four summits and overall from 2001 to 2008, and increased significantly between 2008 and 2015 on all four summits. In SNN, the study region and its individual summits (except MIR) showed significant increases in mean species richness over the entire observation period, but the increase over the periods 2004–2011 and 2011–2019 was only significant on DIE (Fig. 2a, b, Tables S3a, S4a). Over all eight summits, the species richness per SAS increased slightly but significantly by 0.01 species per year (GLMM, p = 0.006, Table S5a). Irrespective of time, species richness decreased significantly by 0.2 species per 100 m of altitudinal gain (GLMM, p < 0.0001, Table S5a) over the eight summits. They showed, however, no significant differences between the upper and the lower sections of a summit (i.e. the SASs 0–5 m versus 5–10 m below the highest summit point, Table S5a), nor between the aspects (N, E, S, W; Tables S5a, S6a).

Fig. 2
figure 2

Vascular plant species richness on GLORIA summits in the Sierra Nevada, Spain. Mean richness (± standard error SE) per a study region (SNE: Sierra Nevada—West and SNN: Sierra Nevada—Northeast), each including four summits consisting of eight summit area sections, respectively, and b per summit. Summit elevation (m) is colour-coded. Different lowercase letters in (a) denote significant differences between the survey years within a study region based on generalised mixed-effects models (glmmTMB; Brooks et al. 2017). In SNN the first survey took place in 2004 and 2006 and was therefore set to the year in-between (2005) in (a). For further details see Tables S3a, S4a and S5a

The mean number of colonisation events per SAS increased significantly on each of the four summits in SNE from the first to the second period. In SNN, colonisations remained constant, except on summit MIR, which showed a significant decrease from the first to the second period (Fig. 3a, b, Tables S3b, S7a). Over all eight summits, colonisations per SAS decreased significantly with increasing elevation with − 0.14 species per 100 m of altitudinal gain (GLMM, p = 0.042), but did not change over the entire timespan (2001–2019), aspects or position of SAS (Tables S5b, S6b).

Fig. 3
figure 3

Colonisation and disappearance events of vascular plant species on the GLORIA summits in the Sierra Nevada, Spain. Mean values (± standard error SE) of (a, b) colonisation and (c, d) disappearance events per year for each (a, c) study region (SNE: Sierra Nevada—West and SNN: Sierra Nevada—Northeast) and (b, d) summit, calculated over eight summit area sections per summit of the four summits per study region. The number of colonisation and disappearance events is calculated for each period between subsequent surveys, and divided by the number of years between surveys to account for different lengths of periods. Summit elevation (m) is colour-coded. Different lowercase letters in (a) and (c) denote significant differences between the periods within a study region, based on generalised mixed-effects models (glmmTMB; Brooks et al. 2017). For further details see Tables S3b, c, S5b, c, and S7

The mean number of disappearance events per SAS decreased significantly from the first to the second period on three summits in SNE (except on the highest summit MAC with the lowest number of disappearances in both periods) and in the whole study region. In SNN, there was no change in disappearances from the first to the second period (Fig. 3c, d, Tables S3c, S7b). Over all eight summits, disappearances per SAS decreased significantly with increasing elevation with -0.15 species per 100 m of altitudinal gain (GLMM, p = 0.042), but did not change over the entire timespan (2001–2019), aspects or position of SAS (Tables S5c, S6c).

The mean species cover sum per quadrat in SNE was 11.89 dm2 in 2001, 10.01 dm2 in 2008, and increased significantly to 14.73 dm2 in 2015. In SNN, the mean species cover sum per quadrat increased significantly from 14.79 dm2 at the first survey (2004, 2006) to 19.78 dm2 in 2011, and to 20.25 dm2 in 2019 (Fig. 4a, b, Tables S3d, S4b). Over all eight summits, the mean species cover sum per quadrat increased slightly but significantly with 0.02 dm2 per year (GLMM, p < 0.0001, Table S8a). Along the elevation gradient, the mean species cover sum per quadrat decreased significantly with 0.4 dm2 per 100 m of altitudinal gain (GLMM, p < 0.001, Table S8a), without significant differences among aspects (Table S6d).

Fig. 4
figure 4

Vascular plant cover and thermic indicator on GLORIA summits in the Sierra Nevada, Spain. Mean (± standard error SE) of (a, b) cover sums (dm2) and (c, d) thermic indicator (Gottfried et al. 2012) for each (a, c) study region (SNE: Sierra Nevada—West and SNN: Sierra Nevada—Northeast) and (b, d) summit, calculated over 64.1 m × 1 m quadrats per study region (16 per summit). Summit elevation (m) is colour-coded. Different lowercase letters in (a) and (c) denote significant differences between the survey years within a study region, based on generalised mixed-effects models (glmmTMB; Brooks et al. 2017). In SNN the first survey took place in 2004 and 2006 and was therefore set to the year in-between (2005) in (a) and (c). For details see Tables S3d, e and S4b, c

The thermic indicator per quadrat in SNE increased significantly from 2001 to 2008. In 2015, the thermic indicator was intermediate between the values of the two preceding surveys. In SNN, the thermic indicator per quadrat remained constant over all surveys (Fig. 4c, d, Tables S3e, S4c). Along the elevation gradient over all eight summits, the thermic indicator decreased significantly with − 0.18 per 100 m of altitudinal gain (GLMM, p < 0.001, Table S8b). The thermic indicator did not change among the aspects (Table S6e) and not over time (2001–2019; Table S8b).

Species groups

Forb-hemicryptophytes represented the highest number of species in each survey, followed by graminoids and shrubs (Table S1, Fig. S2a). The number of colonisations of shrub species was not different to other growth forms (Table S9a), whereas the number of disappearances of shrub species was significantly lower than that of hemicryptophytes (Table S9b). Shrubs had the highest total cover in all surveys, followed by graminoids and hemicryptophytes (Table S1, Fig. S2b). The average cover of shrub species was significantly higher than the cover of other growth forms (Table S9c) and increased slightly but significantly with 0.02 dm2 per quadrat and year (GLMM, p = 0.002). The proportion of shrub species of the total cover sum increased marginally significantly with 1.6% per year (GLMM, p = 0.064).

Twelve species out of the species pool of 103 species were spiny. Spiny species colonised (Table S9a) and disappeared (Table S9b) significantly less often than non-spiny ones, but, over time, their cover significantly exceeded that of non-spiny species (Table S9c). The significant increase of the cover of spiny species was 0.03 dm2 per quadrat and year (GLMM, p < 0.001). However, the proportion of spiny species of the total cover sum did not change significantly (GLMM, + 1% per year, p = 0.203).

In total, 40 Sierra Nevada endemics and 63 non-endemic species occurred on the summits. In SNE, the mean cover of endemic species per quadrat decreased from 5.1 dm2 in 2001 to 4.3 dm2 in 2008 and increased to 5.6 dm2 in 2015. In SNN, the mean cover of endemic species per quadrat increased from 3.0 to 4.4 dm2 in 2011, and decreased slightly to 4.0 dm2 in 2019 (Table S1, Fig. S3). Over all eight summits, there were significantly more disappearance events of endemic species than of non-endemic species (Table S9b), but no difference in the number of colonisations (Table S9a), nor in species cover (Table S9c). The cover of endemics did not change over the entire time span (GLMM, 0.005, p = 0.24). However, the proportion of endemic species of the total cover sum decreased significantly by 2.4% per year (GLMM, p = 0.002).

Nine out of 103 species were hygrophilous. There were no differences in colonisations, disappearances and cover between hygrophilous/non-hygrophilous species (Table S9).

Relationship with climatic factors

In all cases, models using a 5-year or a 7-year period prior to the surveys had the lowest AIC (Table S10).

Rising mean annual air temperatures between the periods prior to the surveys were highly positively related to changes in species richness (GLMM, estimate = 0.88, p < 0.0001), but there was no significant relationship with colonisation, disappearance or changes in cover sums (Table S10).

Increasing mean annual precipitation sums between the periods prior to the surveys were highly significantly and positively related to changes in species richness (GLMM, estimate = 0.005, p < 0.0001), colonisations (GLMM, estimate = 0.003, p < 0.0001), and changes in cover sums (GLMM, estimate = 0.005; p < 0.0001), and negatively with disappearances, i.e. increasing precipitation led to a decrease in species losses (GLMM, estimate = − 0.002, p < 0.0001; Table S10).

An increasing number of days with snow cover per year of the periods prior to the surveys were positively related to changes in species richness (GLMM, estimate = 0.005, p = 0.001), and colonisations of species (GLMM, estimate = 0.003, p < 0.0001). Changes in disappearances and in cover sums of species did not show significant relationship with increasing snow cover duration (Table S10).

Discussion

Contrary to previous findings (Fernández Calzado and Molero 2013; Pauli et al. 2012), we found an increase in vascular plant species richness and vegetation cover on alpine summits of the Sierra Nevada in Spain. Generally, in European mountain systems, rising species numbers were found far more often than decreasing ones. Repeated surveys of historical vascular plant species inventories on alpine to subnival summits of temperate, boreal and arctic Europe even showed an acceleration of increase in species richness during the last decades in synchrony with rapid warming (Steinbauer et al. 2018). Increases in species numbers were mainly attributed to upward range shifts of species, which were found to be greatest where the highest levels of warming were observed in a global study across different organism groups (Chen et al. 2011). Upward species shifts as response to climate warming were also reported for Mediterranean alpine grassland communities in Sierra de Guadarrama, Spain, resulting in a general increase in species richness (historical surveys from 1962, 1977 and 1983 resurveyed in 2012; Jiménez-Alfaro et al. 2014). However, on Mediterranean summits of Europe (including Sierra Nevada), species numbers decreased from 2001 to 2008. Combined effects of warming and a reduction in precipitation in Mediterranean areas were hypothesized to be the main drivers of species losses (Fernández Calzado and Molero 2013; Pauli et al. 2012). Our study shows not only an increase in species richness after 2008, but a significant relationship between increasing species richness and rising temperatures, increasing precipitation, and more days with snow cover. Previous studies suggested that temperature alone is insufficient to explain diversity changes in a water-limited environment (Crimmins et al. 2011). For instance, consistent and negative relationships between drought stress and species richness (Giménez-Benavides et al. 2007; Walck et al. 2011), and positive effects of snow cover duration (Niittynen et al. 2018; Wipf and Rixen 2010) have been reported. After a general downward trend of precipitation during the second half of the twentieth century (Sinoga et al. 2011), precipitation tended to rise in the Sierra Nevada during the last decade according to ERA5 data used in this study. However, this has to be treated with caution, because (1) the rather coarse resolution of our precipitation data must be taken into account, (2) we could not evaluate the actual increase in precipitation with local data because the available high-elevation station data showed gaps, measuring errors or did not cover the recent years, and (3) a considerable inter-annual variability of precipitation patterns has been reported (Norrant and Douguedroit 2006; Polo et al. 2019). Published time series do not extend to 2019, however, data in Polo et al. (2019) at least indicate a precipitation increase between the two survey periods of SNE (2001–2008 and 2008–2015). As an additional measure of water availability, we used the number of days with snow cover, derived from temperature loggers, positioned at our plots. Both, coarse precipitation data and the duration of snow cover, clearly indicate a significant positive relationship between species numbers and water supply in the Sierra Nevada.

Colonisation and disappearance events

Increasing number of colonisation events, combined with a decreasing number of disappearances were found more often on the lower summits. The lower alpine belt of Sierra Nevada is populated by more species (Fernández Calzado et al. 2012) and thus the species pool for upwards migration is larger, similar to the situation in other European mountains (Bruun et al. 2006; Vittoz et al. 2010). Colonisation showed a highly significant positive relation to the variables for water supply: precipitation and duration of snow cover. While longer growing seasons were suggested as one of the main factors for higher community productivity and associated species distribution changes under climate change, and at the same time implicate an increased risk of frost damage (Bueno de Mesquita et al. 2020; Wipf and Rixen 2010), a higher snowpack generates an increased water supply during the vegetation period, which is especially relevant in water-limited environments. Seed dormancy, germination and establishment are mainly driven by temperature and water supply (Bernareggi et al. 2016; Walck et al. 2011), so that in water-limited environments climate change may have an immediate effect on seedling success, in addition to a higher probability of a reduced reproductive success (Giménez‐Benavides et al. 2018). The high seedling mortality due to summer drought is considered a primary constraint for recruitment in Mediterranean mountains (Cavieres et al. 2005; Giménez‐Benavides et al. 2018). In turn, Crimmins et al. (2011) documented a significant downward shift in the optimum elevation of mountain plants in Mediterranean California, following increasing water availability that outpaces evaporative demands. In line with this, species colonisation on our summits in Sierra Nevada could represent horizontal shifts from microrefugia (Scherrer and Körner 2010) in response to ameliorating climatic water-deficit. Besides, colonisations could also result from upper range shifts of lower-elevation species (Jiménez-Alfaro et al. 2014; Rumpf et al. 2019).

Species disappearances were highly significantly negatively related to an increasing amount of precipitation, i.e., the number of disappearances decreased with increasing precipitation. However, disappearance events were not related to changes in temperature or days with snow cover. Disappearance events in high-mountain areas can indicate retracting lower range margins, as expected through warming-driven range shifts (Giménez‐Benavides et al. 2018; Ugarte et al. 2019) or competitive displacements (Lenoir and Svenning 2015). Rumpf et al. (2018) reported from the Alps that the trailing lower margins of species ranges are contracting at least as rapidly as the leading edges are expanding towards higher elevations. In the central Alps, cold-adapted species, being restricted to the high-elevation zones, have experienced strong losses in cover in the lower part of their distribution range (Lamprecht et al. 2018; Steinbauer et al. 2020). In the Sierra Nevada, endemic species disappeared significantly more often than more widespread species. The majority of the Sierra Nevada endemics almost exclusively occurs in the upper vegetation belts (Fernández Calzado et al. 2012) and thus most of them are low-temperature specialist species. The relative contributions of different processes leading to disappearances were probably different in the climatically contrasting mountains of the Alps and the Sierra Nevada, even though apparently driven by climatic changes in both. Water supply, however, is confirmed to be a crucial component for the survival or disappearance of species in periods of drought.

Changes in cover and community-weighted thermic indicator

A thermophilisation of alpine vegetation, i.e. an increase of the ratio of warmth-demanding to cold-adapted species, was detected on the SNE summits only during the first seven-year period (2001–2008; Gottfried et al. 2012), accompanied by an insignificant decrease in the mean total species cover. Similar to the situation in the central Alps, where accelerating thermophilisation occurred with progressing losses of vegetation cover over a 20-year period (Lamprecht et al. 2018), thermophilisation in Sierra Nevada resulted from slightly stronger losses of cold-adapted than gains of warmth-demanding species. In the high central Alps, however, this was far more pronounced with strong and consistent losses of all subnival–nival species (Steinbauer et al. 2020). Periods without signals of thermophilisation (2004/2006–2011 in SNN; 2008–2015 in SNE), in contrast, showed increases of total species cover. Species cover did not show a relation to temperature, which generally was rising in the Sierra Nevada (Pérez-Luque et al. 2016), but a significant positive relationship with increasing annual precipitation. This would indicate a high sensitivity of plant growth to changes of water availability, however, an increase of days with snow cover was not related with an increase in species cover. Impacts of different snow cover durations are only little studied in Mediterranean mountain environments, but positive effects of meltwater on plant growth can hardly be negated (Giménez‐Benavides et al. 2018; Ugarte et al. 2019). In temperate and boreal mountains, an earlier snowmelt extends the humid growing season, whereas in Mediterranean mountains, it would expand the dry summer period and thus the period of potential drought stress. Besides, positive effects of an extended growing season are hypothesised to be counteracted by the detrimental effects of an increasing frequency and intensity of frost (Choler 2018; Klein et al. 2018).

Changes in endemic species and shrubification

The cover of endemic species did not change over time, however, the proportion of endemics of the total cover sum decreased significantly by 2.4% per year. This is in line with declines in endemic dry alpine grassland specialists in alpine vegetation assemblages of central Spain, where also an ongoing colonisation of shrubs into alpine areas was found (Jiménez-Alfaro et al. 2014). Similarly, we found not only a continued increase in shrub cover, but also an increase of their proportion of the total cover sum. This increase of either locally common but endemic shrub species dominating the vegetation zone between the treeline and high-mountain dry grasslands (i.e., Cytisus galianoi, Genista versicolor, Thymus serpylloides subsp. serpylloides), or generally widespread shrub species (i.e., Hormathophylla spinosa, Juniperus sabina) indicates a consistent expansion of their upper distribution ranges to higher vegetation zones (cf. Fernández Calzado and Molero 2013), leading to a ‘shrubification’ of high-elevation habitats of the Sierra Nevada. Shrub advances are known from the alpine belt of the central Apennines and the southwestern Alps (Stanisci et al. 2014). Unfortunately, studies on changes in the abundance of high-mountain species are far rarer than on species numbers, and even more underrepresented in Mediterranean ecosystems (Giménez‐Benavides et al. 2018), which is particularly unfortunate, since Mediterranean high-mountain regions host many locally endemic and at the same time cold-adapted species (Blanca et al. 1998; Pauli et al. 2012). The vulnerability of these species through climate warming increases with the growing risk of summer drought, caused by the rise in average summer air temperature and a reduction in annual rainfall (Nogués Bravo et al. 2008).

Land use changes

Change in land use, including reductions in livestock density or abandonment, which may have significant effects on the performance of high-mountain plants, is unlikely as an explanation for the current observed vegetation changes. While livestock grazing in Sierra Nevada, based on traditionally established rights and thus tolerated even in the core zone of the national park, tends to be more relevant in lower areas and in depressions with snowbed vegetation and shows slight declining tendencies (National park staff personal communication), the population of wild ungulates (in particular Iberian ibex, Capra pyrenaica hispanica) increased tenfold between 1960 and 2012, but individual numbers are considered to be more or less stable over the last twenty years (Granados et al. 2020). This may imply an overall pressure on the vegetation, yet, less so for spiny species, such as abundantly occurring Hormathophylla spinosa, Arenaria pungens and Festuca indigesta, having xeromorphic leafs with spiny stems or tips. Indeed, over time the cover sum of spiny species increased more than the cover sum of non-spiny species and they disappeared significant less than non-spiny species. We therefore cannot exclude any interference through ibexes, but positive trends also of the non-spiny species during the supposedly less dry periods, as well as stable numbers of ungulates over the last twenty years, suggests, that the influence was not a dominant factor. By the way, there were no obvious impacts of ibex grazing discernible in the plots. Tourism was also rapidly growing, with the number of visitors tripled since the establishing of the national park in 1999 (MITECO 2020). However, trampling is expected to reduce cover of most species and only may promote some trampling-resistant graminoids (Barros and Pickering 2015). A new phenomenon of touristic donkey or horse riding in the national park is noticed with concern, as this activity not only damages alpine vegetation, but also supports eutrophication of high-mountain areas as well as the dispersal of low-land plant seeds (Barros et al. 2020); for example, Urtica dioica grows on the summit of Mulhacén, the highest peak of the Sierra Nevada (pers. obs. Lamprecht and Pauli 2016). However, the GLORIA summits are less attractive secondary peaks situated off the main riding trails.

Limitations

Limitations of the datasets and results discussed above are mainly associated with the poor availability of local climate data, in particular indicators for water supply. This study therefore had to be based on rather large-scaled climate data, thus, we must be cautious with interpretations of the underlying mechanisms behind the observed plant diversity changes. New well-managed weather stations in the higher areas are required for a better understanding of the impacts of both temperature and precipitation, as suggested by Pérez-Palazón et al. (2018). Further, investigations in stress physiology of alpine plant species, using experimental studies (e.g., warming manipulation and irrigation; Pugnaire et al. 2020) or biotic proxies, such as dendrochronology (e.g., Juniperus sabina; García-Cervigón Morales et al. 2012), would help to assess and interprete vegetation changes in relation to key climate components.

Conclusion

Plant species composition and diversity on Sierra Nevada summits are strongly influenced by climate change, where both warming and water supply are crucial components for understanding and predicting biodiversity losses and gains. Between 2001 and 2008, the retreat of low-temperature specialist species, which are to a great proportion Sierra Nevada endemic, was rapidly leading to a thermophilisation of the summit vegetation. Vice versa, in periods with more precipitation, vegetation seems to ‘recover’, since it increased in cover, but to the disadvantage of the endemics. For the twenty-first century, not only a pronounced warming is projected (Pérez-Palazón et al. 2018), but also a severe decrease in mean precipitation and an increase in interannual variability for most Mediterranean regions (Nogués Bravo et al. 2008). Giorgi (2006) considered the Mediterranean region to be the largest ‘climate change hot-spot’ in the world. If these climate projections come true, our results suggest that detrimental effects on this water-limited high-mountain flora of the Sierra Nevada are to be expected. Linked to that, Mediterranean ecosystems are among the ‘hyper-hot candidates’ for conservation support (Myers et al. 2000). In this context, continued vegetation monitoring and in situ climate measurements are of high priority for an effective evaluation of the status of ecosystems and their biodiversity, of changing processes and their relevance for conservation. Rapid responses of Mediterranean alpine species indicate a tight synchronisation with climate changes, notably with water availability, which reinfores concerns of endemic biodiveristy losses in the view of regional climate projections.