Vegetation responses to snow cover removal and rainfall reduction in subalpine grasslands: Insights from a 4-year field experiment in the Spanish Pyrenees

ABSTRACT

Reduced precipitation as well as warming may result in less snow accumulation in seasonally snow-covered areas, leading to lower minimum soil temperatures and more frequent and severe soil frosts. Therefore, plant stress is increased not only by drought in warmer months, but also by increased exposure to frost in cooler periods. We conducted a 4-year field experiment to evaluate the effects of diminished snowpack accumulation (through snow removal) and rainfall reduction (through rain-out shelters) on aboveground plant productivity, diversity and species composition of two subalpine grassland plant communities from the central region of the Spanish Pyrenees. We found that the snow removal treatment decreased minimum soil temperature by 0.5°C. Plant diversity decreased by 16 percent, although this effect was only observed in one of the grasslands studied. Aboveground primary productivity seemed to be unaffected. In contrast, we found that the rainfall reduction treatment negatively affected aboveground productivity of leguminous forb species, yet no effect on plant diversity was observed. Both treatments were important drivers of changes in plant species composition. Overall, our results suggest that the resilience of subalpine grasslands to snow cover removal and rainfall reduction treatments may depend on the specific community composition and dominant plant groups.

KEYWORDS:

• Aboveground biomass
• plant diversity
• community composition
• leguminous forbs
• climate change

Introduction

Anthropogenic climate change constitutes one of the major threats to the functioning, management and conservation of the Earth’s ecological systems. In addition to the increase in atmospheric CO₂ levels and subsequent global warming (Blunden and Boyer 2022), climate change leads to an hydrological intensification, which involves the alteration of precipitation patterns worldwide (Dore 2005; Marvel and Bonfils 2013; Seneviratne et al. 2021). In this regard, both observations and model projections indicate a generalized increase in the frequency, magnitude and duration of precipitation extremes such as droughts (Dai 2013; Seneviratne et al. 2021). Understanding how ecosystems will respond to these extreme events is critical to properly anticipate the potential impacts of climate change. Accordingly, the effects of extreme droughts have been extensively examined in most terrestrial ecosystems, typically through manipulative rainfall experiments (e.g., by using rain-out shelters to impose artificial droughts; Knapp et al. 2015; Hoover, Wilcox, and Young 2018; Loik et al. 2019). However, these experimental studies often focus on the immediate effects of drought (Hoover, Wilcox, and Young 2018), while those that consider the effects of drought over multiyear study periods are much scarcer (Evans et al. 2011; Hoover, Wilcox, and Young 2018; Xu et al. 2021; Batbaatar et al. 2022). Therefore, it still remains unclear whether ecosystem responses to drought are immediate or progressive over time (Xu et al. 2021; Batbaatar et al. 2022).

There is considerable interest in understanding the effects of precipitation extremes in grassland ecosystems, as they cover a large proportion of the global land surface (Dixon et al. 2014) and provide essential ecosystem services (e.g., forage, carbon storage, etc.; Bengtsson et al. 2019). Specifically, how grassland primary productivity responds to changes in precipitation has been extensively assessed (e.g., see Wilcox et al. 2017; Zhang et al. 2021 for recent meta-analyses). Overall, grassland productivity is more sensitive to increased precipitation than to extreme drought (i.e., positive asymmetry; Wilcox et al. 2017; Zhang et al. 2021; Wang et al. 2022), particularly regarding aboveground net primary productivity. In addition, the sensitivity of grassland productivity can be related to mean annual precipitation. In mesic temperate grasslands, aboveground plant productivity often exhibits high stability under precipitation extremes, including droughts, whereas more xeric grasslands usually show more dramatic responses (Knapp et al. 2015; Stuart-Haëntjens et al. 2018; Muraina 2020; Liang et al. 2021; but see Byrne, Adler, and Lauenroth 2017). Furthermore, the response of grassland productivity to drought can be largely determined by the identity and traits of the dominant species (Evans et al. 2011; Byrne, Adler, and Lauenroth 2017). For example, forbs (i.e., herbaceous plants that are not graminoids), especially leguminous forbs, tend to experience greater reductions in aboveground biomass due to extreme drought, followed by poorer productivity recovery, than other more resource-acquisitive plant groups such as graminoids (Hoover, Knapp, and Smith 2014; Mackie et al. 2019; Tello-García et al. 2020; Xu et al. 2021; but see Zhang et al. 2021). Consequently, long-term droughts may become an important driver mediating shifts in plant community diversity and composition (Hoover, Knapp, and Smith 2014; Griffin-Nolan et al. 2019; Tello-García et al. 2020), to the point that extreme droughts may have more pronounced effects on plant community diversity and composition than on primary productivity (Batbaatar et al. 2022). Therefore, studies of multiyear drought events are needed to define generalized response patterns of plant diversity and to complement this knowledge.

In alpine and subalpine environments, where seasonal snow cover is common, reduced precipitation may also imply reduced snowpack accumulation (Rumpf et al. 2022). Snow cover insulates the soil from cold air, so a decrease in snow depth leads to lower minimum soil temperatures and more severe soil freezing (Brown and DeGaetano 2011; Wilson et al. 2020), ultimately affecting grassland productivity and composition. For example, grassland plant species may experience increased freeze damage in perennial vegetative structures due to increased exposure to frost (Gu et al. 2008; Lubbe and Henry 2020; Lubbe, Klimešová, and Henry 2021). However, the specific response to increased soil freezing may depend on plant traits related to bud depth and belowground organs of perennation. In this regard, hemicryptophytes (i.e., plants with overwintering buds on the soil surface) appear to be relatively vulnerable to increased freezing temperatures, especially tap-rooted and rhizomatous species, whereas geophytes (i.e., plants with overwintering buds/storage organs belowground) may be less exposed to frost stress (Lubbe and Henry 2020; Lubbe, Klimešová, and Henry 2021). In addition, reduced snow cover alters many ecosystem processes that can also affect plant productivity, diversity and composition. For example, more frequent soil freeze-thaw cycles due to reduced snow cover and earlier snowmelt induce shifts in soil fauna and soil microbial functioning and composition (Gavazov et al. 2017; Coq and Ibanez 2019; Garcia et al. 2020; Broadbent et al. 2021) and, subsequently, in nutrient cycling (e.g., reduced N-availability due to earlier nitrification, loss of soil N during snowmelt, increased plant-available P, reduced soil respiration, etc.; Buckeridge et al. 2010; Gao et al. 2021). Although earlier snowmelt and higher temperatures may benefit certain species by extending the growing season, this effect may still be counteracted by increased freezing due to reduced snow cover (Choler 2018). Therefore, in alpine and subalpine grassland ecosystems, reduced precipitation may increase plant stress not only due to drought in warmer months, but also due to frost in cooler periods (Henry et al. 2018). However, there is a lack of studies in the literature that assess vegetation responses to the prolonged and combined effects of reduced precipitation and snowpack accumulation in alpine and subalpine grassland ecosystems.

In this study, we conducted a 4-year field experiment to assess the effects of diminished winter snowpack accumulation (through snow removal) and rainfall reduction (through rain-out shelters) on aboveground plant productivity, diversity and species composition of two subalpine grassland plant communities in NE Spain, and to determine whether these effects are immediate or progressive over time. We were also particularly interested in examining the interaction between the effects of snow cover removal and rainfall reduction, in order to disentangle, which might be the main driver mediating future changes in subalpine grassland vegetation. We hypothesized that (i) both forms of drought negatively affect subalpine grassland productivity and diversity, especially when acting in combination. Specifically, we expected that (i.1) rainfall reduction would cause a greater decline in productivity within leguminous forbs than within more resource-acquisitive groups such as graminoids, and that (i.2) snow cover removal would cause a generalized decline in productivity and diversity, as most plant species in our study grasslands were hemicryptophytes. Accordingly, we also hypothesized that (ii) both forms of drought are important determinants of shifts in plant species composition in subalpine grasslands. Finally, we expected to find both immediate (i.e., effects observed from the first year) and progressive (i.e., cumulative effects over the study period) vegetation responses to snow cover removal and rainfall reduction.

Material and methods

Study area

This study was conducted in two subalpine grasslands from the central region of the Spanish Pyrenees (NE Spain). This mountain range, which rises to over 3000 m, forms a natural border between France and Spain. Overall, precipitation in the study area has decreased by about 2.5% per decade over the past 50 years, particularly in summer and winter (OPCC-CTP 2018). In addition, a decrease in snowpack accumulation has been observed from 1950 to the present (OPCC-CTP 2018). Future model projections indicate that the average snowpack thickness in the central region of the Pyrenees could be reduced by half by the year 2050, and that the duration of snow cover would be reduced by more than one month (OPCC-CTP 2018).

The very long history of agro-pastoral activity in the region has shaped the subalpine landscape into a mosaic of open areas (i.e., grasslands) alternating with areas of denser vegetation (i.e., forests; Montserrat and Fillat 1990; Ruiz and Lasanta 2018; Canals 2019). Specifically, the study grasslands were located in the Torla and Pineta valley bottoms (Figure S1A), both within the Ordesa and Monte Perdido National Park. They have had a similar land management regime since medieval times, being created, used and maintained mainly for extensive livestock production (Montserrat and Fillat 1990; Ruiz and Lasanta 2018; Canals 2019). Although the study grasslands are currently rarely grazed due to land abandonment, woody plant species are infrequent because of low forest regeneration. Calcareous soil substrate predominates in the study area, where cambisols and phaeozems develop. The topography is similar in both study sites, although the Pineta valley is narrower and more sheltered. Despite Torla’s higher elevation, Pineta was characterized by lower minimum temperatures (Table 1). Mean annual precipitation was also slightly lower in Pineta than in Torla (Table 1). Winters are snowy at both study sites (Figures S2 and S3).

Table 1. Location, elevation, climatic parameters (i.e., mean potential solar irradiation, mean of maximum and minimum temperatures, mean snow cover and mean annual rainfall) and percentage of rainfall blocked by rain-out shelters at each study site.

Study site	Location		Elevation (m.a.s.l.)	Mean solar irradiation (kWh m^−2 d^−1)	Mean Tmax* (°C)		Mean Tmin* (°C)		Mean snow cover** (cm)	Mean annual rainfall* (mm)	Rainfall reduction*** (%)
	Latitude	Longitude			Jan	Jul	Jan	Jul
Torla	42.651	−0.066	1363	4.84	5.96	23.41	−3.04	9.17	38.84	1405.4	19.55
Pineta	42.657	0.119	1239	4.73	6.29	23.57	−3.39	8.67	37.19	1356.6	32.82

*Mean value for the period 1971–2020; data obtained from the Digital Climatic Atlas of Aragón (https://www.aragon.es/-/atlas-climatico-de-aragon; Cuadrat, Saz, and Vicente-Serrano 2007).

**February average (peak month) for the period 2000–2015; data obtained from the nearest weather station.

***Equivalent to the driest year of the last 100 years.

The study grasslands supported different plant communities (Table 2). Specifically, the grassland vegetation in Torla consisted of nonleguminous (e.g., Fragaria vesca L., Galium verum L. and Plantago media L.) and leguminous (e.g., Lotus corniculatus L.) forbs, accompanied by some grasses (e.g., Brachypodium pinnatum (L.) P.Beauv.). In Pineta, the grassland plant community was composed mainly of graminoids (e.g., Agrostis capillaris L., Brachypodium pinnatum (L.) P.Beauv. and Carex sp.) and nonleguminous forbs (e.g., Epipactis palustris (L.) Crantz).

Table 2. Grassland biotic parameters (i.e., mean species richness, Shannon diversity index and relative abundance of graminoids, leguminous forbs, nonleguminous forbs and woody species) at the subplot level per study site.

Study site	Richness	Shannon index	Relative abundance*
			Graminoids (%)	Legumes (%)	Non-legumes (%)	Woody (%)
Torla	20 ± 0.66	2.56 ± 0.03	18.48 ± 2.92	19.84 ± 2.57	59.61 ± 3.51	2.07 ± 1.58
Pineta	11.67 ± 0.93	1.51 ± 0.13	68.75 ± 5.41	7.48 ± 2.24	23.59 ± 4.18	0.18 ± 0.09

*As computed from species plant cover.

Experimental design

We conducted a 4-year field experiment within the framework of the International Drought Experiment (IDE; Knapp et al. 2017). Specifically, six 3 × 3 m plots were established within an enclosure of homogeneous vegetation at each study site in the autumn of 2015. Rain-out shelters were then installed over three randomly selected plots at each site (Figures S1(b)) and remained in place throughout the experiment. Rain-out shelters consisted of metal structures, with an angled roof partially covered with transparent and v-shaped plastic strips (Figure S1C), which passively reduced rainfall events by a fixed amount. The amount of rainfall blocked by rain-out shelters was calculated independently for each study site (Table 1; Figures. S1C, S2 and S3) to impose a level of drought (in terms of mean annual precipitation) equivalent to the driest year of the last 100 years (Knapp et al. 2017). In addition, two 1 × 1 m subplots were established within each plot (n = 2 × 6 × 2 = 24). One of the subplots was randomly assigned for snow cover removal treatment, and the other subplot was left as a control (Figure S1B). To minimize potential disturbance to the soil and plant litter layer, a plastic mesh with 1 cm gaps was spread over the ground. Snow removal was then performed at least once every 15 days, especially following major snowfall events, from October/November (i.e., beginning of the snow season) until the end of February (Figure 1C). Shovels and/or brooms were used to remove snow down to the level of the plastic mesh. Late winter snow cover was maintained to minimize treatment effects on the growing season soil moisture and spring phenology (Rixen et al. 2022). Finally, a soil temperature sensor was buried at a depth of 2 cm below the soil surface in the center of each subplot. In addition, concurrent with the vegetation surveys, soil moisture measurements were taken in each subplot from the 0–15 cm depth layer (depth at which most of the root biomass is located; Jackson et al. 1996) using an HH2-type moisture meter (Delta-T Devices, Cambridge, England).

Figure 1. (a) February’s minimum soil temperature (mean ± SE) and (b) July’s soil water content (mean ± SE) per study site, treatment and year. Horizontal lines indicate mean values per treatment and year. Symbol shape indicates study site (circle and square for Torla and Pineta, respectively), symbol/line color indicates snow removal treatment (gray and black for no snow removal and snow removal, respectively) and symbol/line pattern indicates reduced precipitation treatment (solid and dashed for no reduced precipitation and reduced precipitation, respectively). C, control, gray solid line; S, snow removal, black solid line; R, reduced precipitation, gray dashed line; S + R, snow removal plus reduced precipitation, black dashed line.

Vegetation surveys

Data on vegetation cover and aboveground biomass was collected annually over the study period (i.e., 2016–2019) at the peak of the growing season (i.e., end of July). In particular, the percentage of aboveground cover was recorded for every plant species found at the subplot level. Visual estimation of plant species cover was facilitated by the use of a sampling quadrat divided into 100 squares. Plant diversity was then computed as the effective number of species from the Shannon-Wiener index (Jost 2006), using plant cover as a surrogate for species abundance. This can be interpreted as the number of species of equal abundance required to produce a given value in the traditional Shannon-Wiener index, which is easy to interpret (e.g., a plant community with an effective number of species of 4 is twice as diverse as another community with an effective number of species of 2). On the other hand, aboveground plant biomass was harvested from a 0.2 × 1 m strip at each subplot, thus avoiding resampling the same surface over the 4-year period. Biomass samples were taken to the laboratory, sorted by species and classified into three groups: graminoids (i.e., perennial grasses and sedges), leguminous and nonleguminous forbs. Finally, plant biomass samples were oven-dried at 60°C for 48 h prior to weighing.

Statistical analyses

The effects of snow removal and rainfall reduction treatments were evaluated by fitting linear mixed models (LMMs) for repeated measures with the lmer function of the lme4 package (Bates et al. 2015) in R software (R Core Team 2021). More specifically, minimum soil temperature of the coldest month (i.e., February), soil water content of the warmest month (i.e., July), aboveground dry biomass of graminoids, leguminous and non-leguminous forbs, and the effective number of species were specified as the dependent variables in their respective models. Treatment, study site and time, as well as their interactions, were set as fixed terms. Treatment was included as a four-level factor (i.e., C, control subplots with neither snow removal nor rainfall reduction; S, subplots with snow cover removal; R, subplots with rainfall reduction; and S + R, subplots with the combination of snow removal and rainfall reduction), and time was included as an ordinal variable. Subplot nested within plot were set as random intercept effects in the models to account for the split-plot design. Selection of the best random effect structure was based on Akaike’s Information Criterion (>2 AICc points difference; Zuur et al. 2009). Data on soil water content was arcsine-transformed, data on graminoids and leguminous forbs dry biomass was log-transformed and data on nonleguminous forbs dry biomass was square-root-transformed to meet model assumptions. Temporal autocorrelation of residuals was checked for all models based on acf plots and Durbin-Watson tests. As a result, an order 1 correlation structure was specified in the soil water content model.

The significance of differences in plant species composition among treatments, study sites and time was assessed by performing nonmetric multidimensional scaling (NMDS) ordination of Bray-Curtis dissimilarities and permutational multivariate analysis of variance (perMANOVA; 999 permutations) with the metaMDS and adonis functions, respectively, from the vegan package (Oksanen et al. 2020) in R software (R Core Team 2021). The study site was specified as a blocking factor to constrain permutations. Additionally, post-hoc pairwise comparisons were performed to test for significant differences in species composition among specific pairs of treatments at each study site by using the pairwise.adonis2 function.

Results

Change in minimum soil temperature and soil moisture in response to snow cover removal and rainfall reduction

We found that minimum soil temperature at 2 cm depth differed significantly among treatments (Table 3). A post-hoc Tukey’s HSD test indicated that minimum soil temperature was, in general, about 0.5°C lower in subplots where snow cover was removed (Figure 1A). This effect was similar in both study sites, as the interaction between treatment and study site was not significant (Table 3). Moreover, we found a significant interaction between treatment and time to explain differences in minimum soil temperature (Table 3; Figure S4A). Specifically, the difference in minimum soil temperature between subplots where snow cover was removed and those where snow cover was maintained increased from no difference at the beginning of the study period, to approximately 1°C difference at the end (Figure 1A; Figure S4A). Finally, we found that the interaction between the study site and time was also significant (Table 3). The minimum soil temperature decreased slightly in Pineta over the study period, while it rose about 1°C in Torla (Figure 1A). The absolute minimum soil temperature recorded in February was −4.26°C in Torla (2018) and −3.4°C in Pineta (2019).

Table 3. Results of statistical analyses for the effect of treatment (i.e., C, S, R and S + R), study site (i.e., Torla and Pineta) and time on minimum soil temperature of the coldest month (i.e., February), soil water content of the warmest month (i.e., July), aboveground productivity, plant diversity (i.e., effective number of species) and species composition. Significant effects (p < .05) are highlighted in bold. Residuals Df = 83.

	Treatment influence						Vegetation responses
	Minimum soil temperature (R^2 = 0.47)			Soil water content (R^2 = 0.6)			Graminoids biomass (R^2m = 0.38)			Leguminous forb biomass (R^2m = 0.81)			Non-leguminous forb biomass (R^2m = 0.91)			Plant diversity (R^2m = 0.68)			Species composition (R^2 = 0.52)
	Df	Fvalue	p (>F)	Df	Fvalue	p (>F)	Df	Fvalue	p (>F)	Df	Fvalue	p (>F)	Df	Fvalue	p (>F)	Df	Fvalue	p (>F)	Df	Fvalue	p (>F)
Treatment	3	12.84	<.01	3	0.11	.95	3	0.53	.66	3	1.96	.13	3	0.33	.80	3	.32	.81	3	2.15	<.01
Study site	1	0.04	.85	1	96.88	<.01	1	4.25	<.05	1	59.88	<.01	1	94.28	<.01	1	35.13	<.01	1	71.10	<.01
Time	1	5.67	<.05	1	24.68	<.01	1	0.74	.39	1	9.31	<.01	1	9.99	<.01	1	3.08	.08	1	3.66	<.01
Treatment:Site	3	0.58	.63	3	0.02	.99	3	1.27	.31	3	0.63	.60	3	1.22	.32	3	3.12	<.05	3	2.11	<.01
Treatment:Time	3	12.78	<.01	3	0.18	.91	3	0.67	.57	3	3.37	<.05	3	0.30	.82	3	0.20	.90	3	0.29	.99
Site:Time	1	13.27	<.01	1	14.11	<.01	1	0.76	.39	1	0.02	.89	1	6.67	<.05	1	0.99	.32	1	2.55	<.01

We did not find any significant effect of snow cover removal and rainfall reduction treatments on soil water content (Table 3). We only found a significant effect of the interaction between the study site and time on soil water content in the 0–15 cm depth layer (Table 3), where soil water content in July diminished over the study period in Torla, while it remained fairly constant in Pineta (Figure 1B).

Change in aboveground productivity and diversity in response to snow cover removal and rainfall reduction

Nonsignificant effects of treatment, time or their interaction, were observed on the aboveground productivity of graminoids species (Table 3; Figure 2A). We found that it differed only between study sites (Table 3). In particular, the dry biomass of graminoids was significantly higher in Pineta than it was in Torla (Figure 2A).

Figure 2. (a) Graminoids aboveground dry biomass (mean ± SE), (b) leguminous forbs aboveground dry biomass (mean ± SE), (c) non-leguminous forbs aboveground dry biomass (mean ± SE) and (d) plant diversity (effective number of species, mean ± SE) per study site, treatment and year. Horizontal lines indicate mean values per treatment and year. Symbol shape indicates study site (circle and square for Torla and Pineta, respectively), symbol/line color indicates snow removal treatment (gray and black for no snow removal and snow removal, respectively) and symbol/line pattern indicates reduced precipitation treatment (solid and dashed for no reduced precipitation and reduced precipitation, respectively). C, control, gray solid line; S, snow removal, black solid line; R, reduced precipitation, gray dashed line; S + R, snow removal plus reduced precipitation, black dashed line.

There was a significant interaction effect between treatment and time on aboveground productivity of leguminous forb species (Table 3). Although the aboveground dry biomass of leguminous forbs increased over the study period in subplots under ambient rainfall, it did not change in subplots with rainfall reduction (Figure 2B; Figure S4B). On the other hand, we found that aboveground dry biomass of leguminous forb species also differed significantly between study sites (Table 3). In particular, the dry biomass of leguminous forbs was significantly higher in Torla than in Pineta (Figure 2B).

Regarding aboveground productivity of nonleguminous forb species, we did not find significant effects involving snow cover removal nor rainfall reduction treatments (Table 3; Figure 2C). Furthermore, we found that aboveground dry biomass of nonleguminous forbs was significantly higher in Torla than in Pineta (Table 3; Figure 2C). Moreover, the interaction term between study site and year was significant in explaining the differences in nonleguminous forb productivity (Table 3). In particular, the aboveground biomass of nonleguminous forbs species increased slightly in Torla during the study period, while it remained fairly constant in Pineta (Figure 2C).

Finally, we found a significant interaction between treatment and study site explaining differences in grassland plant diversity (Table 3). A post-hoc Tukey’s HSD test revealed that, in Torla, subplots with snow cover removal (i.e., S) had significantly lower plant diversity than control subplots (i.e., C), while the remaining treatments had intermediate values (Figure 2D; Figure S4C). Conversely, nonsignificant differences in plant diversity were found among specific pairs of treatments in the Pineta study site (Figure 2D; Figure S4C). On the other hand, we found that grassland plant diversity was significantly higher in Torla than in Pineta (Table 3; Figure 2D). The interaction effects between treatment and time, and between study site and time, were nonsignificant (Table 3).

Change in species composition in response to snow cover removal and rainfall reduction

We found that plant species composition varied significantly among treatments (Table 3). Interestingly, there was a significant interaction effect between treatment and study site (Table 3), indicating that differences in plant species composition between a specific pair of treatments depended on the study site. More specifically, post-hoc pairwise comparisons revealed that, in Torla, subplots under ambient rainfall (i.e., C and S) on the one hand, and subplots with rainfall reduction (i.e., R and S + R) on the other, shared species composition. Accordingly, the 95 percent confidence ellipses overlapped following such subplot grouping (Figure 3). In the Pineta study site, significant differences in species composition were only found between C and S + R subplots (Figure 3). The interaction term between treatment and time was nonsignificant (Table 3), indicating that differences in plant species composition among treatments responded similarly over the whole study period (Figure 3). Finally, plant species composition also changed significantly between study sites (Table 3). Accordingly, the NMDS plot showed a clear segregation among subplots located in Torla and Pineta, respectively (Figure 3).

Figure 3. Nonmetric multidimensional scaling (NMDS) ordination plot of Bray-Curtis dissimilarities among subplots’ species composition. Symbol shape indicates study site (circle and square for Torla and Pineta, respectively), symbol/line color indicates snow removal treatment (gray and black for no snow removal and snow removal, respectively) and symbol/line pattern indicates reduced precipitation treatment (solid and dashed for no reduced precipitation and reduced precipitation, respectively). C, control, gray solid line; S, snow removal, black solid line; R, reduced precipitation, gray dashed line; S + R, snow removal plus reduced precipitation, black dashed line. Higher drawing intensity indicates more recent data and vice versa. Ellipses indicate 95 percent confidence intervals.

Discussion

In this work, we investigated the vegetation responses of two subalpine grassland communities from the central region of the Spanish Pyrenees to diminished snowpack accumulation (through snow cover removal) and reduced rainfall (through rain-out shelters) over a 4-year field experiment. First, in line with previous studies, we found that snow cover removal treatment resulted in lower minimum soil temperatures (Brown and DeGaetano 2011; Gavazov et al. 2017; Henry et al. 2018). However, the observed decline in our study was about 0.5°C, while an average drop of 2°C in minimum soil temperature has been reported from globally distributed temperate grasslands following the same experimental protocols (i.e., IDE; Knapp et al. 2017; Henry et al. 2018). More interestingly, unlike the large-scale study carried out by Henry et al. (2018), our multiyear study period allowed for the evaluation of whether treatment effects were immediate or progressive over time. Specifically, we found that the effect of snow cover removal on minimum soil temperature increased over the study period. The effect of snow cover in insulating soil from low ambient temperatures is positively related to snowpack thickness (Zhang 2005; Wilson et al. 2020). In our study, the mean snowpack depth in 2016 was lower than the mean quantity of snow accumulated in subsequent years (Figs. S2 and S3). This could have diminished the effect of the snow cover removal treatment during the first experimental year, helping to explain our finding and highlighting the importance of multiyear studies to properly evaluate the effects of reduced snowpack accumulation in these environments.

We did not find any generalized effect of snow cover removal treatment on the aboveground grassland productivity. This result contrasts with the increasing and progressive effect of snow cover removal on minimum soil temperature. At a global scale, a negative relationship has been established between lower minimum soil temperatures from reduced snowpack accumulation and the effect of snow cover removal on aboveground biomass (Henry et al. 2018). Therefore, the lack of a treatment effect might simply be explained by the relatively small reduction in minimum soil temperature observed in snow removal subplots. In addition, plants from alpine and subalpine environments have different physiological mechanisms to escape, avoid and tolerate freezing. For example, some species can produce and accumulate non-structural carbohydrates (e.g., oligosaccharides of the fructans family) that enhance freezing tolerance (Valluru and Van den Ende 2008). Accordingly, low magnitude soil frosts would only slightly affect the perennation organs of the most vulnerable or rarest plant species, which in turn would have little impact on the aboveground productivity of grasslands (Stover and Henry 2019). Interestingly, plant mechanisms against freezing and drought can be very similar (Volaire et al. 2023), which would help explain why we found limited evidence of a combined effect of snow removal and rainfall reduction treatments in our study.

Alternatively, the reduction in minimum soil temperature due to snow cover removal may have been enough to increase freezing damage in some plant species. In addition, snow removal might have also increased the number of soil freeze/thaw cycles (Figure S5), which would increase plant stress and susceptibility to soil freezing. This could be the case for some plant species in the Torla study grassland, whose cover is reduced (e.g., Achillea millefolium L.) or even disappears from the community (Table S1). On the other hand, there were also species that responded positively (e.g., Plantago media L. or Centaurea nigra L.; Table S1). Temperate mesic grasslands are usually non-water-limited and their vegetation dynamics are strongly driven by competition (Grime 1973; Michalet et al. 2006). Consequently, plant species better adapted to freezing could expand under more favorable conditions (e.g., less competition for space or light) during the growing season, thus potentially counteracting the detrimental effects of increased frost exposure on aboveground productivity.

In fact, in accordance with previous studies reporting that decreasing snow cover affects diversity in biotic communities (Chen et al. 2008; Gavazov et al. 2017; Niittynen, Heikkinen, and Luoto 2020), we found that snow cover removal reduced grassland plant diversity. This finding provides evidence of changes in plant species richness and/or evenness, supporting the idea introduced above and pointing out that, in subalpine environments, reduced snowpack accumulation could be a more relevant driver of plant diversity loss than of a general decline in aboveground grassland productivity, especially when larger differences in snowpack accumulation and/or minimum soil temperature reduction might be expected under future climate change (Rixen et al. 2022).

Curiously, we did not find a significant effect of the snow cover removal treatment on Pineta grassland diversity. Indeed, there seemed to be no drastic changes in species dominance or richness among treatments at this site (Table S1). Such a lack of effect could be related to the community specific composition. In this sense, the grassland community in the Pineta study site was dominated by Carex sp. (Table S1). It has been reported elsewhere that Carex species (i.e., Carex curvula All.) can occur on cold sites without snow cover, can resist very low soil temperatures and have a high freezing resistance (von Büren and Hiltbrunner 2022). In addition, we suspect that the water table at Pineta might be shallower than at the Torla study site, which might have ultimately driven vegetation dynamics rather than the treatments. In support of this possibility, we found that soil temperature and water content changed over the study period in the Torla study grassland, but not in Pineta. Similarly, a shallower water table might also help to reduce an increased number of freeze/thaw cycles resulting from snow cover removal (Figure S5).

In line with previous studies, we found that extreme rainfall reduction imposed for four consecutive years negatively affected the aboveground productivity of leguminous forb species (Carlsson et al. 2017; Tello-García et al. 2020). Moreover, this effect seemed to accumulate progressively over the study period, as differences between treatments ultimately increased with time. Leguminous forbs constitute a key group of plant species in many grasslands around the world. For instance, they play a prominent role in the nitrogen cycle due to their ability to fix atmospheric nitrogen and enhance soil nutrient availability (Craine et al. 2002; Rodríguez, Canals, and Sebastià 2022). Legumes are also important in determining forage quality because of their low C:N ratio and the amount of protein that they provide (Craine et al. 2002; Amiri and Shariff 2012; Elgersma, Søegaard, and Jensen 2014) compared to other species groups.

Therefore, under future climate change scenarios in which extreme droughts become more intense and frequent, a progressive reduction could be anticipated in some of the ecosystem services provided by these grasslands. Nevertheless, we found that reduced rainfall treatment had no significant effect on grassland plant diversity. This finding makes sense, since leguminous forb species were the least abundant species group in our study sites. In contrast, we did not find any effect of rainfall reduction treatment on aboveground productivity of graminoid and non-leguminous forb species. These species, particularly graminoids, can be characterized by rapid growth and efficient use of resources (e.g., in water uptake; Armas and Pugnaire 2011; Pueyo et al. 2016). Therefore, even though the rainfall reduction treatment was equivalent to the most severe drought in the last 100 years, it is likely that these plant species were not stressed by water scarcity.

We did not find lower soil moisture in the 0–15 cm depth layer under reduced rainfall treatment. This may lead to thinking that the observed effect of reduced rainfall treatment on aboveground productivity of leguminous forb species might be due to undesirable side effects of rain-out shelters. However, it has been documented that infrastructure effects (e.g., blocking solar radiation or altering the microclimate under the shelter) are negligible or minimal (Yahdjian and Sala 2002; Loik et al. 2019). Moreover, these minimal effects may eventually appear when the shelters are expected to intercept a relatively high amount of rainfall (e.g., >60 percent; Yahdjian and Sala 2002; Loik et al. 2019), while in our case the rainfall interception was less than 30 percent (Figure S1C; Table 1).

Therefore, we discarded this possibility. On the other hand, the lack of effect of rain-out shelters on soil moisture could simply be explained by the fact that rainfall reduction plots behave as “islands” surrounded by a normal rainfall landscape. It is possible that subsurface water movement may have slightly reduced the effectiveness of the rainfall reduction treatment (Kröel-Dulay et al. 2022). Nonetheless, it should also be noted that our measurement of soil moisture may not be fully informative because it was measured in the driest month (Figs. S2 and S3), whereas it is possible that the largest differences in soil moisture between drought treatments occur immediately after rainfall events (Kreyling et al. 2008).

Finally, we found that both snow cover removal and rainfall reduction treatments altered grassland species composition and points to their potential role in jeopardizing the conservation of subalpine grassland ecosystems. Curiously, the major driver mediating changes in species composition differed between study sites. In Torla, significant differences in species turnover were determined by the rainfall reduction treatment, while in Pineta these differences were only found with combined snow cover removal and rainfall reduction treatments. This finding could be explained based on the composition of the initial plant communities at each site.

In Pineta, grasses and other graminoids were the most abundant group of species, which are known to better tolerate changes in precipitation (Hoover, Knapp, and Smith 2014; Mackie et al. 2019; Tello-García et al. 2020; Xu et al. 2021), while forbs were more abundant in Torla. Interestingly, effects on plant species composition were not associated with substantial changes in grassland plant diversity, which has also been reported in other temperate grasslands (Batbaatar et al. 2022). Additionally, we found that the response of grasslands regarding plant species composition was immediate. That is, differences among treatments were observed in the first year, while they remained constant over the study period. Further studies are needed to determine the extent to which these compositional changes may be irreversible. For example, it is necessary unraveling what are the thresholds of rainfall and snow cover reduction, or how long a given reduction in rainfall or snow cover is enough to make it difficult to restore the current conditions of diversity and ecosystem services of alpine and subalpine grasslands.

In conclusion, our results highlighted that, in the future context of climate change, reduced precipitation could be more relevant than diminished snow cover in determining progressive changes in aboveground plant productivity, especially in those grasslands dominated by leguminous forb species. Conversely, lower snowpack accumulation could be more relevant than rainfall reduction in determining changes in grassland plant diversity. Both, rainfall and snow cover reduction can be important drivers of rapid changes in plant species composition, although their relative importance may be linked to the specific traits of the dominant plant species groups (e.g., forbs vs graminoids). Relatively small vegetation responses in our subalpine grasslands were associated with relatively small effects of treatments on abiotic soil conditions (i.e., minimum soil temperature and soil moisture), suggesting that more dramatic grassland vegetation responses might be expected under greater climate change impacts.

Finally, we found little evidence of the combined effects of snow cover removal and rainfall reduction treatments in our subalpine grassland vegetation. Further research is needed to elucidate the exact causes of this apparent lack of interactive effects and to distinguish the effects of snow removal and rainfall reduction over time due to different weather conditions between years from those due to the extension of the experimental study period itself. In addition, future studies should further investigate the effects of reduced snow cover and precipitation on soil freeze/thaw cycles and its relationship to the vegetation response in alpine and subalpine grassland ecosystems.

Authors’ contributions (CRediT)

Antonio I. Arroyo: Conceptualization, Formal analysis, Visualization, Writing – Original Draft; Miguel Castillo-Garcia: Investigation, Writing – Review and Editing; Yolanda Pueyo: Conceptualization, Writing – Review and Editing; Concepción L. Alados: Conceptualization, Investigation, Formal analysis, Writing – Review and Editing, Funding acquisition.

Data availability

Data supporting the results of the study will be made available in an appropriate public repository (e.g., Dryad [https://datadryad.org/stash] or Digital CSIC [https://digital.csic.es/]) upon article acceptance.

Acknowledgments

This work was partially funded by the Spanish National Research Council (Proyecto Intramural Especial, PIE, no. 202230E013). M. Castillo-Garcia was supported by the Spanish Ministry of Science, Innovation and Universities (grant FPU16/05508; FPU Program 2016). We would like to thank J. Escós, P. Sánchez-Navarrete and M. Gartzia for their assistance with the fieldwork and data collection. We are also especially grateful to the Department of Agriculture, Farming and Natural Environment of Gobierno de Aragón for granting permission to carry out the study in a protected area.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15230430.2023.2292339

Additional information

Funding

This work was supported by the Consejo Superior de Investigaciones Científicas [PIE, no. 202230E013]; Ministerio de Ciencia, Innovación y Universidades [FPU16/05508; FPU Program 2016].