Snowbed communities and soil C and N dynamics during a four-year investigation in the NW-Italian Alps

ABSTRACT

The snowbed habitats represent a relevant component of the alpine tundra biome, developing in areas characterized by a long-lasting snow cover. Such areas are particularly sensitive to climate changes, because small variations in air temperature, rain, and snowfall may considerably affect the pedoclimate and plant phenology, which control the soil C and N cycling. Therefore, it is fundamental to identify the most sensitive abiotic and biotic variables affecting soil nutrient cycling. This work was performed at seven permanent snowbed sites belonging to Salicetum herbaceae vegetation community in the northwestern Italian Alps, at elevations between 2,686 and 2,840 m.a.s.l. During a four-year study, we investigated climate, pedoclimate, floristic composition, phenology, and soil C and N dynamics. We found that lower soil water content and earlier melt-out day decreased soil N-NH₄⁺, N-NO₃⁻, dissolved organic carbon (DOC), dissolved organic nitrogen (DON), total dissolved nitrogen (TDN), microbial nitrogen (Nmicr), microbial carbon (Cmicr), and C:Nmicr ratio. The progression of the phenological stages of Salix herbacea reduced soil N-NH₄⁺ and increased DOC. Our results showed that the snow melt-out day, soil temperature, soil water content, and plant phenological stages were the most important factors affecting soil biogeochemical cycles, and they should be taken into account when assessing the effects of climate change in alpine tundra ecosystems, in the framework of long-term ecological research.

KEYWORDS:

• LTER
• standard monitoring variables
• Salix herbacea
• topsoil
• phenology
• carbon and nitrogen
• alpine tundra

Introduction

The presence of a substantial and long-lasting snow cover holds a prominent relevance for the ecology of a large part of the Earth’s biomes, especially at high latitude and in mountain regions; in particular, 63 percent of the land masses in the Northern Hemisphere are at least temporarily covered by snow (Hammond, Saavadra, and Kampf 2018). In fact, the seasonal snowpack depth, duration, and melt-out control the onset and duration of the growing season in cold regions (e.g., Cooper, Dullinger, and Semenchuk 2011; Freppaz et al. 2018), mitigate negative soil temperature, and affect microbial activity (e.g., Mikan, Schimel, and Doyle 2002; Freppaz et al. 2018). In addition, it affects soil nutrient cycling (e.g., Hiltbrunner, Schwikowski, and Körner 2005; Filippa et al. 2010), soil gas fluxes (e.g., Filippa et al. 2009; P. D. Brooks et al. 2011), pedogenesis (e.g., Holtmeier and Broll 1992), and other soil properties (e.g., Hiller et al. 2005; Pintaldi et al. 2016; Freppaz et al. 2018), along with plant composition (Carbognani, Tomaselli, and Petraglia 2014).

Snowbed plant communities represent a relevant component of the alpine tundra biome (Björk and Molau 2007), developing in small areas where winter snow accumulates with great thickness, mainly due to topographical factors and wind redistribution (e.g., Gjaerevoll 1956; Billings and Bliss 1959; Johnson and Billings 1962). Snowbeds have a dense vegetation cover, and one of the most relevant species is the dwarf shrub Salix herbacea L. (Quaglia et al. 2020), which, despite the small size of the aerial parts, can develop an impressive and extremely ramified root system, showing also pronounced longevity (Wijk 1986). These habitats, located above the treeline, are characterized by specific environmental properties, such as a limited growing season, extended periods with air temperatures below freezing, and prolonged snow-covered soils (Edwards, Scalenghe, and Freppaz 2007; Knowles, Blanken, and Williams 2015). Thus, snow represents a key ecological factor in mountains due to its capability to shape the spatial patterns of soil properties and alpine vegetation (Braun-Blanquet 1932; Gjaerevoll 1956; Bliss 1963; Scott and Billings 1964; Freppaz et al. 2018). However, the predicted increase in temperature (IPCC 2021) makes seasonally snow-covered regions, such as the Alps, particularly sensitive to climate change (Cebon et al. 1998; Beniston 2006), because small changes in air temperature and/or precipitation may result in large modifications of the pedoclimate and plant phenology, with potential effects on soil C and N cycling. In particular, global warming might introduce a number of different effects on the seasonal snow cover; for instance, increasing air temperatures can result in a higher percentage of rain falling instead of snow, causing a thinner and shorter-lasting snow cover, a phenomenon already visible in many mountain ranges (e.g., Serquet, Marty, and Rebetez 2013; Notarnicola 2020). As remarked by Hock et al. (2019), higher spring and summer temperatures lead to earlier melt-out dates and longer growing seasons for vegetation. Snowmelt timing is known to influence soil water availability during the summer months (Wilson et al. 2020). A warmer climate will likely deplete soil moisture in the summer by increasing evapotranspiration but also indirectly by advancing the timing of seasonal snowmelt (Blankinship et al. 2014). Snowpack changes can thus have cascading effects through ecosystems due to their influence on soil temperatures and water availability (Wilson et al. 2020).

Variations in the thickness and duration of snow cover result in remarkable differences not only in pedoclimatic conditions and nutrient cycling but also in the composition (D. A. Walker et al. 1993; Fisk, Schmidt, and Seastedt 1998) and phenology of plant communities (Khorsand Rosa et al. 2015).

Snowbed communities (i.e., Salicetum herbaceae Braun-Blanquet and Jenny 1926) are especially vulnerable to the forecasted decrease in snowfall and change in snow cover duration, coupled with increased temperature (Carbognani, Petraglia, and Tomaselli 2012; Bender, Lehning, and Fiddes 2020). In fact, the moment of snowmelt and the related beginning of the snow-free period sets the start and the duration of the whole growing season, conditioning plant survival, distribution, growth, and phenology (D. A. Walker et al. 1993; Stanton, Rejmanek, and Galen 1994; Huelber et al. 2006; Quaglia et al. 2020). For instance, Quaglia et al. (2020) noted that the interannual variability of the snowmelt date and the length of the snow-free season markedly affected the temporal distribution of the phenological phases, and Petraglia et al. (2014) showed, through a manipulation experiment, that an earlier snowmelt time determines an earlier occurrence of the flowering date of several snowbed species. A longer growing season could make grassland species more competitive than snowbed ones, thus changing the community to grasslands dominated species (Matteodo et al. 2016). In turn, these biological responses can alter biogeochemical cycles and ecosystem services, causing an increased decomposition of carbon stored in soils (Schuur et al. 2008), with a loss of biodiversity (Pauli, Gottfried, and Dullinger 2012) and changes in trophic dynamics (Post and Pedersen 2008; Wookey et al. 2009).

However, because the uptake of elements depends both on the phenology of plants and the availability of the resource, it is not easy to discriminate between the contribution of the factors affecting the availability of soil nutrients (e.g., pedoclimatic conditions) and the plants acquisition of nutrients. For instance, nitrogen (N) supplies, including inorganic (NH₄⁺ and NO₃⁻) and organic (dissolved organic nitrogen, DON) forms, show a pronounced seasonality, which is run by variations of soil temperature and moisture levels, whereas its uptake is related to different phenological phases of plants (Tian et al. 2017). Because the early stages such as foliation and flowering are the most demanding for the majority of alpine plants, during this phenological phase the nutrients uptake is the highest, whereas it reduces near the end of the plant cycle (Adair and Burke 2010). The balance of nutrient uptake, assimilation, and losses through turnover, leaching, exudation, and other disturbances (e.g., grazing, parasites, etc) affect nutrient concentrations (such as nitrogen) in plant biomass (Chapin and Shaver 1989; Aerts and Chapin 2000; Eckstein and Karlsson 2001). Güsewell (2004) noted that nitrogen concentration in herbaceous plant biomass was positively correlated with the Ellenberg et al. (1991) nutrient (N) indicator value, which can be considered as an estimator of nutrient requirements (Ryser and Urbas 2000). This indicator value, among others, belongs to the wide family of ecological indicator values (EIVs) widely used by plant ecologists because they can characterize the ecological condition of a site using species composition without the need of direct measurements (Diekmann 2003).

For the flora of Central Europe and for the Alps, Ellenberg et al. (1991) and Landolt et al. (2010), respectively, attributed to each plant species a numerical value referred to an ordinal range (Ellenberg: from one to nine; Landolt from one to five) to express the ecological optimum of a species for a given environmental factor (e.g., for soil moisture, nutrients, reaction, etc). For instance, high average N values suggest that vegetation communities occur in nutrient-rich conditions, whereas low average N values suggest nutrient-poor conditions (Güsewell, Peter, and Birrer 2012; Pittarello et al. 2016). The accuracy of such EIVs lies in the knowledge and extensive field experience of botanists and ecologists (Ellenberg et al. 1991; Pignatti et al. 2005; Landolt et al. 2010). Moreover, EIVs are a valid tool for habitat estimation and could be used to indirectly describe the corresponding pedoclimatic properties (Schaffers and Sýkora 2000), which could be affected by changes in climate conditions and influence the soil C and N cycling.

Despite the reliability of such a tool, few studies utilized EIVs, and few of them were performed in the snowbed habitat (e.g., Robbins and Matthews 2014; Matteodo et al. 2016). As supported by several authors (e.g., Björk and Molau 2007; Schöb et al. 2008; Carbognani, Petraglia, and Tomaselli 2012; Ernakovich et al. 2014), we hypothesize that the predicted change in the amount and timing of snowpack due to climate change will influence snowbed ecosystems, altering pedoclimatic condition and phenology, which in turn will affect the C and N soil biogeochemical cycles. Furthermore, given that snowbeds are more vulnerable to climate change effects than other subalpine and alpine plant communities (Matteodo et al. 2016), they are often considered as a target habitat where the effects of global warming on alpine vegetation can be thoroughly explored (Björk and Molau 2007; Schöb et al. 2008). Therefore, for four years we studied the floristic composition, phenology, and soil C and N forms of seven homogeneous permanent plots (all belonging to Salicetum herbaceae communities), identifying standard variables most sensitive to environmental changes, supporting the long-term ecological research in high-elevation ecosystems. In particular, our work aimed at (1) evaluating the interannual variability (2016–2019) of climatic and pedoclimatic conditions, soil C and N forms, and plant phenology dynamics; (2) assessing the effect of climatic, pedoclimatic, and vegetation variables (i.e., plant phenology and EIVs) on soil C and N dynamics; and (3) evaluating the considered abiotic/biotic drivers as possible monitoring variables of the response of C and N biogeochemical cycles to climate change in the framework of the long-term ecological research.

Materials and methods

Study area

The study was carried out at the Long-Term Ecological Research (LTER) Site Angelo Mosso Scientific Institute (https://deims.org/17210eba-d832-4759-89fa-9ff127cbdf6e), located in the Western Italian Alps at the foot of the Monte Rosa Massif (4,634 m.a.s.l.), on the border between Piemonte and Valle d’Aosta Regions (northwest Italy; Figure 1). The research was performed at seven permanent study sites (named site 1, 2, 6, 7, 8, 9, 10, according to the LTER site numerations), located between 2,686 m.a.s.l. (site 10) and 2,854 m.a.s.l. (site 6).

Figure 1. Location of the research study area in Italy and in the northwestern Italian Alps (www.pcn.minambiente.it), aerial overview of the study area (orthoimage year 2006; coordinate system WGS84/UTM 32 N), and location of the seven permanent study sites at the LTER Site Angelo Mosso Scientific Institute.

Since 2007, meteorological parameters of the study area were continuously recorded by an Automatic Weather Station located at 2,901 m.a.s.l. and belonging to the Italian Army (Comando Truppe Alpine–Servizio Meteomont). From 2007 to 2019, the area was characterized by a mean annual air temperature of −2.2°C and a cumulative annual snowfall of 805 cm, and the mean annual liquid precipitation was ca. 370 mm. Generally, the snowpack accumulated by late October–early November, and snowmelt started typically in late April to early May. The mean snow cover duration was 263 days, ranging between a minimum of 236 to a maximum of 312 days. The mean soil temperature during the snow-covered season was 0°C, and it was 7°C during the snow free season. Generally, the mean duration of soil freezing lasted for 52 days, ranging between 0 and 162, with a mean temperature of −1°C.

During the time span 2016–2019 (hydrologic years; i.e., 1 October 2015–30 September 2016), the mean annual air temperature recorded by the Automatic Weather Station was −1.8°C, with mean daily values ranging from a minimum of −21.1°C (16 January 2017) to a maximum of +15.8°C (12 June 2018; Figure 2). The mean air temperature during 2017 and 2019 was warmer than the historical time series (2007–2019 time frame; −2.2°C), with values of −1.4°C and −1.5°C, respectively (Table 1). The cumulative liquid precipitation during the snow-free season ranged between 202 mm in 2017 and 380 mm in 2016, with a maximum daily liquid precipitation of 50.0 mm (16 September 2016). In 2017 and 2019 the cumulative liquid precipitation was below the mean value of the historical time series (283 mm; 2007–2019 time frame). The snow cover duration ranged between 254 days (2018) and 295 days (2016; Table 1); 2018 was the only year with a snow cover duration below the mean value of the historical time series (263 days). The earliest snowmelt day occurred in year 2017 (day of year [DOY] 184) and the latest occurred in 2016 (DOY 202). The mean snow density before snowmelt onset was equal to 341 kg m⁻³, and the mean snow water equivalent was around 1,113 mm; the N-NH₄⁺ and N-NO₃⁻ stocks in the snowpack ranged between 0.34–0.97 and 0.53–0.88 kg ha⁻¹ while dissolved organic carbon (DOC) and DON stock ranged between 7.4–10.9 and 0.85–2.34 kg ha⁻¹.

Table 1. Weather conditions recorded by the Automatic Weather Station during the experimental period (2016–2019) and historical time frame (2007–2019).

	2016	2017	2018	2019	2016–2019	2007–2019
Mean air temperature (°C)	−2.1	−1.4	−2.2	−1.5	−1.8	−2.2
Cumulative rain (mm)	380	202	319	253	289	283
SCD (days)	295	265	254	275	272	263
Cumulative snowfall (cm)	756	608	824	666	714	805
MOD (DOY)	202	184	200	193	195	190
MTFS (°C)	9	8	7	9	8	7
DSF (days)	96	0	153	21	68	51

MOD: melt-out day of snow; MTFS: mean soil temperature during snow-free season; DSF: duration of soil freezing. Further details about acronyms are shown in Table 2.

Figure 2. (A) Air temperature (mean daily values), snow depth, and rainfall recorded at the Automatic Weather Station (AWS) from 1 October 2015 to 30 September 2019 and (B) soil temperature (10 cm depth, mean daily values) at the seven study sites.

The bedrock was primarily micaschists, with some inclusions of amphibolites and calcschists. Soils, described and classified according to the Food and Agriculture Organization (2006) and IUSS Working Group WRB (2015), belonged mainly to Regosol (sites 1, 8), Leptosol (site 2), Cambisol (sites 6, 7), and Umbrisol (sites 9, 10) great groups (IUSS Working Group WRB 2015; Supplementary Material, Table S1). Soil total organic C (TOC) and N (TN) ranged from 0.1 to 18.4 percent and from 0.02 to 1.23 percent, respectively; soil pH ranged from 4.2 to 6.1. Detailed soil classification and physical and chemical parameters are summarized in Supplementary Material, Table S1.

The vegetation of the sites was included in the “Siliceous alpine and boreal grasslands” (habitat 6150, according to the EU Habitat Directive), belonging to the association Salicetum herbaceae (Quaglia et al. 2020). The vast majority of the species could be described, based on vegetative optimum, either as snowbed specialists, debris, or grassland species, according to Aeschimann, Lauber, and André Michel (2004).

Site-specific climatic-pedoclimatic variables

For the measurement of hourly soil temperature, each study site was equipped with thermistors combined with dataloggers (HOBO Pro v2 U23-00x, accuracy: ±0.1°C) placed within soil at a 10 cm depth, from fall 2016 until fall 2019. Moreover, several site-specific variables were calculated (Table 2): (1) proportion of cumulative rain days from the melt-out day to each sampling day (PRD; extracted from the Automatic Weather Station data); (1) snow cover duration (SCD), based on each hydrological year (as the following variables), using the daily soil temperature data: when the daily variation of soil temperatures remained below 1°C, the day was considered as a “snow-covered day” (Danby and Hik 2007); (3) melt-out day of snow (MOD), calculated as the first day of each hydrological year in which the daily amplitude of soil temperature went above 1°C; (4) duration of soil freezing (DSF); that is the number of days in which soil temperatures remained below 0°C; (5) number of soil freeze–thaw cycles (FTCs), estimated when the daily mean soil temperature reached values below and above 0°C (Phillips and Newlands 2011); (6) mean temperature of soil during the freezing period (MTF); (7) mean temperature of soil during the snow-free season (MTSF); (8) intensity of soil freezing (ISF); and (9) gravimetric water content of soil (GWC) calculated at every sampling time (Black 1965).

Table 2. Selected climatic, pedoclimatic, and vegetation variables used in generalized linearmixed models (GLMM).

Variables	Term	Definition	Unit
Climatic and pedoclimatic variables
Proportion of raining days	PRD	Proportion of cumulative raining days between melt-out day of snow and sampling day	%
Snow cover duration*	SCD	Sum of “snow-covered days” in each hydrological year	Days
Melt-out day of snow* or snowmelt day	MOD	Date of complete snowmelt (indicated as DOY)	DOY
Duration of soil freezing	DSF	Cumulative number of days, from 1 October to the melt-out day, when mean daily soil temperature <0°C	Days
Soil freeze–thaw cycles	FTCs	Number of soil freeze–thaw cycles in each hydrological year	Number
Mean soil temperature during soil freezing	MTF	Mean daily soil temperature when the soil is frozen, from 1 October to the melt-out day	°C
Mean soil temperature during the snow-free season	MTSF	Mean daily soil temperature between MOD and sampling	°C
Intensity of soil freezing**	ISF	Minimum soil temperature when soil is frozen	°C
Gravimetric water content	GWC	Soil gravimetric water content measured at each sampling time	%
Vegetation variables
Salix herbacea phenophase	SP	Median observed phenophase for each survey date	Number
Landolt moisture indicator value	Fm	Soil moisture based on vegetation preferences weighted with species percentage cover	Number
Landolt nutrient indicator value	Nm	Soil nutrients based on vegetation preferences weighted with species percentage cover	Number
Landolt reaction indicator value	Rm	Soil reaction based on vegetation preferences weighted with species percentage cover	Number

*MOD and SCD are considered pedoclimatic indices because they were calculated from soil temperature data.

**ISF was excluded from the GLMM r > |0.70|.

Soil sampling and analysis

At each site, inside a homogenous snowbed pertaining to the Salicetum herbaceae association, one 9 m² plot was established. Three topsoil samples (A horizon, 0–10 cm depth) were collected within the plot each month during the snow-free season from 2016 until 2019. The three topsoil samples were analyzed separately, as analytical replicates, and then we obtained a single average value for a single plot. More specifically, three soil samplings were performed each year in July, August, and September. Each soil sample was sieved at 2 mm within 24 hours since collection and oven-dried at 105°C for 8 hours to obtain the GWC. An aliquot of 20 g of fresh soil was extracted with 100 mL K₂SO₄ 0.5 M, and 10 g was fumigated using chloroform for 18 hours before extraction with 50 mL K₂SO₄ 0.5 M. The concentration of DOC both in not fumigated and fumigated soil extracts (extractable DOC) was determined with a TOC analyzer (Elementar, Vario TOC, Hanau, Germany) after filtration with 0.45-μm nylon membrane filters. The microbial carbon (Cmicr) was estimated as the difference in extractable DOC between fumigated and nonfumigated samples, corrected using a recovery factor of 0.45 (Brookes et al. 1985). Ammonium (extractable N-NH₄⁺) concentrations in soil extracts were measured spectrophotometrically (U-2000, Hitachi, Tokyo, Japan) using a modified Berthelot method based on the reaction with salicylate in the presence of alkaline sodium dichloroisocyanurate (Crooke and Simpson 1971). Nitrate (extractable N-NO₃⁻) concentrations in soil extracts were measured spectrophotometrically (U-2000, Hitachi, Tokyo, Japan) using the Greiss reaction (Mulvaney 1996) modified according to Cucu et al. (2014). Total dissolved nitrogen (extractable TDN) in the extracts was measured as reported for DOC. Extractable DON was determined as the difference between extractable TDN and inorganic nitrogen (extractable N-NH₄⁺ + N-NO₃⁻) in the extracts. Microbial nitrogen (Nmicr) was estimated from the difference in extractable TDN between fumigated and nonfumigated samples corrected using a recovery factor of 0.54 (Brookes et al. 1985).

Vegetation variables and phenology surveys

From summer 2016 to 2019, surveys were carried out to describe the site vegetation composition and identify the most abundant species. Within each site, a 16 m² (4 × 4 m) square plot (marked by metal pegs) for vegetation and phenology surveys was placed close to the soil sampling plot. Along the two diagonals of each plot, the vertical point-quadrat method (Daget and Poissonet 1971) was used to detect botanical composition. At 5-cm intervals, plant species touching a steel needle were identified and recorded, for a total of 200 measurements per plot. Because occasional species are often missed, a complete list of plant species found in the plot was recorded. The frequency of occurrence of each species, which is an estimate of species canopy cover, was calculated and converted to species percentage cover (%SC; Pittarello et al. 2016); that is, by dividing the frequency of occurrence by 2. A %SC = 0.3 was attributed to all occasional plant species (Tasser and Tappeiner 2005). The demand of each plant species for moisture level (F), soil pH level (R), and soil nitrogen content (N) was estimated according to the Landolt EIVs (Landolt et al. 2010; Supplementary Material, Table S2). Such values range from 1 (low demand) to 5 (high demand). Then, each plant species was associated with the corresponding F, R, and N Landolt value and the mean indexes (Fm, Rm, Nm) were calculated for each plot, separately for each year and site, by averaging species values weighted on their %SC. Plant phenology was recorded for S. herbacea, which not only was the most abundant species across all plots (%SC ranged from 37  to 73 percent, being lower than 50 percent in one single plot), but it had an extensive radical mass (Wijk 1986) when compared with the other species, and it was thus believed to have a greater interaction with soil chemicals. The main plot was divided into sixteen adjacent 1 × 1 m subplots, used as pseudo-replicates, where the phenophase of the most advanced individual was recorded seven to nine times per year throughout the whole growing season (Quaglia et al. 2020). The phenophase was assigned with a visual observation of the individual using an adaptation of the Biologische Bundesanstalt, Bundessortenamt and CHemical industry scale (Hack et al. 1992) fitted for the snowbed plant communities as a reference (Quaglia et al. 2020). This scale is characterized by a series of progressive integer numbers ranging from zero to fifty-nine, where the tens represent the major phenophases in chronological order (shoot emergence, leavf development, flower buds, flowering, fruiting, seed dispersal) and the units represent the advancement of such phenophases. The parameter S. herbacea phenophase (SP) was calculated for every soil sampling date as the linear interpolation of the median of the sixteen subplots of each phenological survey date.

Statistical methods

Statistical analyses were performed using R software v3.6.0 (R Core Team 2019). To assess the relative importance of climatic-pedoclimatic and vegetation variables (i.e., explanatory variables, Table 2), nine generalized linear mixed models (GLMMs) were fitted, one for each response variable related to soil C and N forms: N-NH₄⁺, N-NO₃⁻, DOC, TDN, Cmicr, Nmicr, DON, CNmicr. Before modeling data with GLMMs, a correlation analysis (Pearson coefficient) was performed on explanatory variables to exclude highly collinear variables (r > |0.70|; ISF excluded). Explanatory variables were standardized (Z-scores) to allow the assessment of their size effect by scrutinizing model parameters (β coefficients). To account for the repeated measure structure, the soil sampling period, year, and site were considered as crossed random factors. Because all response variables had continuous and positive values (Y > 0), models were built both with Gaussian (identity- and log-link functions) and Gamma (log-link function) distributions and then the best fitting one was chosen based on the lowest Akaike information criterion (Zuur et al. 2009). For each GLMM, explanatory variables with a minimal effect on the overall model fit were identified with the “drop1” function from “stats” R package and removed from the model afterward. Through this procedure, each factor was iteratively removed at a time while holding the remaining factors constant and, using an analysis of deviance test, a report of the extent to which each factor affected the overall model fit was provided. Model assumptions were checked using “DHARMa” R package v0.2.7. (Hartig 2020). Modeling was carried out using the “glmmTMB” function from the “glmmTMB” R package v1.0.1. (M. E. Brooks et al. 2017).

Results

Interannual variability of site-specific climatic and pedoclimatic conditions

During the snow-covered season, the topsoil temperatures recorded in the study sites were usually close to 0°C (Figure 2), with exceptions in early winters 2015–2016 and 2017–2018. It reached a minimum of −7.5°C in 2018 at site 6. ISF was always classified as mild freezing, except for two events of mild/hard freezing (cf. Tierney et al. 2001) recorded at sites 1, 6, 8 in 2017 (−6.6°C, −7.5°C, −6.8°C) and at site 10 in 2016 (−7.3°C). The lowest MTF (−2.4°C) occurred at site 10 in 2016. DSF ranged between 0 days (site 8, 2019) and 247 days (sites 6, 9, 2018; site 10, 2017). The maximum number of FTCs (two) was recorded at sites 1 and 10 in 2018 and at site 7 in 2019. SCD ranged between 230 (site 8, 2018) and 283 days (site 1, 2016; Table 3). The SCD was weakly positively correlated with the MOD (r = 0.38, p < .001). Earliest MOD occurred at site 6 on DOY 163 (2017) and the latest MOD occurred at sites 9 and 10 on DOY 201 (2018; Table 3).

Table 3. Descriptive statistics for all variables used in the generalized linearmixed models (GLMMs).

	n	Mean	SE	Median Min	Max	Range
Climatic-pedoclimatic variables
PRD	77	0.3±0.02	0.3	0.0	0.7	0.7
SCD	77	261.2±1.67	263.0	230.0	283.0	53.0
MOD	77	184.4±1.24	188.0	163.0	201.0	38.0
DSF	77	167.4±7.26	166.0	24.0	247.0	223.0
FTCs	77	1.1±0.04	1.0	1.0	2.0	1.0
MTF	77	−0.5±0.06	−0.3	−2.4	0.0	2.4
MTSF	77	8.6±0.17	8.8	2.1	11.4	9.3
GWC	77	24.9±0.73	23.1	14.3	43.8	29.5
Soil C and N forms
N-NH4^+	77	3.9±0.27	3.4	0.6	11.8	11.2
N-NO3^−	77	0.5±0.04	0.4	0.1	2.0	1.9
DOC	77	85.1±3.88	76.3	29.4	165.1	135.7
TDN	77	20.2±1.24	18.1	4.4	54.0	49.5
DON	77	15.9±1.25	13.3	1.9	50.5	48.7
Cmicr	77	484.1±26.60	445.8	143.0	1,166.2	1,023.2
Nmicr	77	36.7±2.37	31.4	7.4	86.8	79.3
C:Nmicr	77	14.5±0.57	14.1	5.0	29.4	24.4
Vegetation variables
SP	77	46.1±1.59	51.0	9.0	59.0	50.0
Fm	77	3.5±0.01	3.4	3.3	3.6	0.3
Nm	77	2.3±0.02	2.4	2.0	2.5	0.5
Rm	77	2.2±0.01	2.2	2.1	2.3	0.2

The unit of climatic-pedoclimatic and vegetation variables are reported in Table 2.

Soil C and N forms are expressed as milligrams per kilogram, except for C:Nmicr.

During the snow-free season, the topsoil temperatures reached the maximum of 16.7°C at site 8 in 2019 (Figure 2), and the lowest and highest MTSF were 2.1°C and 11.4°C, respectively, at site 10 (2019).

Rainy days (PRD variable) were not equally distributed in the different snow-free seasons. In particular, in 2016 rainy days were mostly concentrated in September, close to the third sampling, accounting for 40 percent of the total liquid precipitation of the snow-free period. In 2017 rainy days were equally distributed between August and September, during which 60 percent of the total rainfall occurred. Otherwise, rainy days were slightly concentrated in August in 2018, whereas in 2019 liquid precipitation was quite well distributed across the whole snow-free season.

The GWC was between 27 and 28 percent in 2016 and 2017, respectively, whereas in 2018 and 2019 it was lower, between 20 and 23 percent, respectively. The lowest GWC (14 percent) was recorded at site 8 (September 2019) and the highest (44 percent) at site 9 (September 2017). Overall, during the four years of monitoring, site 7 was the driest, with a mean GWC of around 19 percent, whereas sites 2 and 6 were the wettest, with values around 30 percent.

Interannual variability of soil C and N dynamics

On an interannual basis (Figure 3), N-NH₄⁺ reached minimum (0.6 mg kg⁻¹) and maximum (11.8 mg kg⁻¹) values in 2017 (site 8) and 2019 (site 9), respectively (Figure 3); the N-NO₃⁻ content reached a minimum value (0.1 mg kg⁻¹) in 2016 (sites 1 and 7) and 2017 (sites 7 and 8), and the maximum value (2.0 mg kg⁻¹) was recorded in 2019 (site 10). DOC content ranged between 29.4 and 165.1 mg kg⁻¹, both in 2016 (sites 8 and 2, respectively). TDN reached minimum (4.4 mg kg⁻¹) and maximum (54.0 mg kg⁻¹) values in 2017 (site 8) and 2016 (site 6), respectively. DON reached minimum (1.9 mg kg⁻¹) and maximum (50.5 mg kg⁻¹) values in 2019 (site 7) and 2016 (site 6). Cmicr ranged between 143.0 and 1,166.2 mg kg⁻¹, both in 2017 (sites 8 and 9, respectively), and Nmicr reached minimum (7.4 mg kg⁻¹) and maximum (86.8 mg kg⁻¹) values in 2019 (site 7) and 2016 (site 2), respectively. C:Nmicr reached a minimum value (5.0) at site 1 in 2018 and a maximum value (29.4) at site 7 in 2019.

Figure 3. Interannual variability of soil C and N forms in the seven snowbed communities. Sampling periods refer to different moments of sampling, which occurred every year within July, August, and September.

Vegetation and phenology dynamics and interannual variability

Although from 2016 to 2019 sixty plant species were identified within the seven permanent sites, S. herbacea was the most abundant (56.6 percent ± 0.02, mean %SC ± standard error). Overall, the other species had a moderate abundance, such as Alchemilla pentaphyllea L. (11.3 percent ± 0.35), Poa alpina L. (11.2 percent ± 0.02), Luzula alpinopilosa (Chaix) Breistr. (6.5 percent ± 1.68), and Gnaphalium supinum L. (5.0 percent ± 0.02).

The growing cycle of S. herbacea showed some variations throughout the monitoring years, with a perceivable shift in the timing of the events. In general, the first phenophases (from emergence to flowering) proceeded in a fast sequence right after the snowmelt, whereas the late phenophases (fruiting to seed dispersal) occurred at a slower rate (Figure 4). The mean onset of flowering (phenophase 35) was on DOY 212 in 2018, with twenty-one days of delay on the same phenophase in 2017 (Figure 4). The interannual difference in the period of appearance of the pre-flowering phenophases (i.e., less than thirty-five) is greater than the difference in the post-flowering phenophases (more than thirty-five; Figure 4). Beginning of seed dispersal started the earliest in 2017 (DOY 213) and occurred the latest in 2018 (DOY 257) with a maximum forty-four-day gap. Mean seed dispersal started on DOY 217 in 2017 and DOY 247 in 2018, whereas it occurred on DOY 241 in 2016 and DOY 234 in 2019.

Figure 4. Phenology of the species Salix herbacea throughout the four-year monitoring period. Phenophases are indicated on the y axis as a progressive number following the adaptation of the Biologische Bundesanstalt, Bundessortenamt and Chemical industry scale (Hack et al. 1992); day of year is reported as a variable on the x axis. For every year, the lines represent the median phenophase of the seven sites on each survey date.

Influence of climatic-pedoclimatic and vegetation variables on soil C and N forms

GLMMs revealed that several climatic and pedoclimatic variables significantly influenced soil extractable C and N forms in snowbed areas (Figure 5 and Supplementary Material, Table S3), although with a different weight. In particular, considering the ranking of relative importance provided by β coefficients for each soil C and N forms (Figure 5), the GWC resulted overall the most important variable, affecting primarily, with a positive relation, extractable DOC, TDN, Cmicr, and C:Nmicr. Moreover, it affected, as a second variable, N-NH₄⁺ and Nmicr and, as a third variable, N-NO₃⁻. The MOD played a fundamental role as well, resulting in one of the most important variables, influencing primarily, with a positive relation, extractable N-NO₃⁻ and Nmicr. MOD also affected positively, as a second variable, extractable DOC, TDN, and Cmicr and, as third variable, N-NH₄⁺. The DSF influenced positively, as a first and third variable, respectively, extractable DON and TDN. The climatic variable PRD negatively influenced N-NO₃⁻ (as a second variable) and Nmicr (as a fourth variable). Although with a less important ranking, FTCs positively affected, as a second variable, C:Nmicr, whereas MTSF affected, as third variable with a positive relation, both Cmicr and Nmicr. SCD played a marginal role, affecting as a third variable only extractable DOC, with a negative relation.

Figure 5. Standardized effects of climatic-pedoclimatic and vegetation variables on soil C and N forms derived from generalized linear mixed models. For each model, the effects (standardized coefficient) of each independent variable are represented by a point symbol and their confidence intervals by a horizontal line. The independent variables are listed on the y axis: PRD, SCD, MOD, DSF, FTCs, MTF, MTSF, ISF, GWC, SP, Fm, Nm, and Rm. When the horizontal line intersects the vertical black line (representing a null effect), it means that the independent variable has no significant effect. Conversely, when the confidence interval line does not intersect the vertical black line, the independent variable has a significant effect. Significance: *p < 0.05, **p < 0.01, ***p < 0.001; the absence of asterisks means that there is no significant effect.

Concerning the vegetation variables, SP and Nm were the most important variables, whereas the effect of Rm was less pronounced. In particular, SP primarily influenced N-NH₄⁺ with a negative relation, also affecting DOC (positive relation) as a fourth variable. Nm negatively influenced DON and Nmicr, respectively, as second and fifth variables and Rm affected C:Nmicr as a second variable with a negative relation.

Because EIVs are attributed on a scale from 1 to 5 with a 0.5 interval precision, the required variation range for each index significance is 0.5 (Table 3). Nm variation is slightly higher, whereas both Rm and Fm have lower variation rates (Rm min = 2.06, Rm max = 2.30, Fm min = 3.29, Fm max = 3.59), and despite the barely significant relation of Rm with C:Nmicr, neither one provides significant interactions with any other measured soil parameter (Figure 5).

Discussion

Climatic-pedoclimatic drivers on soil C and N forms

In the time span 2016–2019, during both the snow-covered and snow-free seasons, the climatic and pedoclimatic variables had a fundamental effect on soil C and N forms in our research area, strongly influencing the soil biogeochemical cycles of the snowbed communities. In particular, soil GWC during the snow-free season, which ranged between 14.3 and 43.8 percent, had a first-order effect on soil C and N forms: the greater the GWC, the higher the content of N-NH₄⁺, N-NO₃⁻, DOC, TDN, Nmicr, Cmicr, and C:Nmicr. Our results confirmed the importance of soil water content as one of the main variables favoring microbial activity and C and N transformations (e.g., Haynes 1986; Stark and Firestone 1995; Lipson, Schmidt, and Monson 1999; Magnani et al. 2017; Pintaldi et al. 2019), in both snow-covered (Larsen, Jonasson, and Michelsen 2002) and growing seasons (Lipson, Schmidt, and Monson 1999). Indeed, within the climate change perspective, even small changes in water content availability could strongly influence the soil microbial biomass and N availability (Rui et al. 2011; Wang et al. 2014). The results are in agreement with those reported by Magnani et al. (2017) in the same study area but under different vegetation cover, which showed a positive correlation between soil water content and the microbial biomass, N-NH₄⁺, and DOC concentrations, confirming that water content is of great importance to the C and N dynamics in these high-elevation ecosystems. This is in accordance with Thornton et al. (2021), who identified the soil hydraulic properties among the essential mountain climate variables, such as snowmelt.

In our work, the snowmelt timing (MOD) played an important role as well, significantly influencing, with a positive relation, several C and N forms. In particular, our results indicated that an earlier MOD caused a decrease in soil N-NH₄⁺, N-NO₃⁻, DOC, TDN, Nmicr, and Cmicr, therefore representing one of the most important abiotic variables in our snowbed areas. The strong effect of MOD on soil C and N could be explained by the interaction of several variables such as the snowpack nutrient input (Filippa et al. 2010; Freppaz et al. 2018), earlier plant development and nutrient uptake (Cooper, Dullinger, and Semenchuk 2011), and the related competition dynamics with microbial communities (Legay et al. 2013). As demonstrated by several authors (P. D. Brooks and Williams 1999; Hiltbrunner, Schwikowski, and Körner 2005; Petrone, Buffam, and Laudon 2007; Filippa et al. 2010; Freppaz et al. 2018; Westergaard-Nielsen et al. 2020), snowpack is considered a remarkable source of both dissolved inorganic and organic N; therefore, the timing of snowmelt defines the start and duration of the growing season, and water and nutrients released from the snowpack influence soil moisture, nutrient status (Kuhn 2001), and microbial activity until later in the summer (Freppaz et al. 2018). An earlier snowmelt day is a forecasted effect of global climate change (Beniston et al. 2003), which may affect the temporal N dynamics, limiting the important ecosystem functions performed by winter soil microbial communities, such as the breakdown of complex litter residues and associated buildup of high soil DON and available NH₄⁺ (Broadbent, Snell, and Michas et al. 2021). Thus, in a warming climate scenario, we can assume an overall future decrease in soil available nutrients and microbial biomass in more snow-free high-elevation snowbed habitats.

Based on our results, the DSF significantly affected DON and TDN with a positive relation. The results agree with those reported by Pintaldi et al. (2019), although in a subalpine forest catchment, where DSF significantly influenced, with a positive relation, both organic and inorganic N forms. As shown by several authors (Kalbitz et al. 2000; Grogan et al. 2004; Groffman et al. 2011; Viglietti et al. 2014), DSF can cause an increase in soil organic labile forms, which is mainly attributable to the physical disruption of the litter layer (and/or rootlets) and the consequent leaching from the disrupted materials. The cryogenic disruption of aggregates (Kværnø and Øygarden 2006) and microbial mortality (Groffman et al. 2001) can have the same effects. Furthermore, Vestgarden and Austnes (2009) found that prolonged soil freezing increased soil C and N forms as well, probably partly because of starvation and exhaustion and increased lysis of microorganisms. Less snow and a later start of protective snow cover associated with climate change can be associated with deeper frost penetration in soils at the beginning of winter, thus increasing soil freezing effects on soil nutrient cycling and microbial communities. In fact, the effect of colder soils in a warmer world is well known (Groffman et al. 2001). Higher labile C and N contents at the beginning of the snow-free season, associated with early leaching, will also lead to a nutrient loss in these high-elevation soils.

The number of FTCs was a weak variable, affecting only C:Nmicr with a positive relation. This result can be explained by the fact that, during the four years of monitoring, a very limited number of FTCs were recorded. In addition, snowbed areas are, by definition, characterized by long-lasting snow cover (even eight to ten months; e.g., Hiller et al. 2005), which decouples soil temperature from the extreme cold air temperatures (Körner 2003).

The climatic variable rainfall days (PRD) significantly influenced both N-NO₃⁻ and Nmicr, with a negative relation. The negative relation between PRD and N-NO₃⁻ was probably linked with leaching phenomena governed by precipitation (e.g., Egli, Mirabella, and Sartori 2008), as well as with a sort of dilution effect of infiltrating water caused by increasing rain (e.g., Balestrini et al. 2019; Pintaldi et al. 2019), even if this variable is not related with the GWC of each plot, which seems to favor the nitrification processes. Regarding the relation between PRD and Nmicr, a similar trend was reported in the same area but for different soil types by Freppaz et al. (2019). Although this variable is not necessarily related with the GWC of each plot, it is reasonable to suppose that greater PRD may also indicate a greater soil water content (at least in the short term), which may reduce the microbial N assimilation, even if the site-specific GWC was positively related with the Nmic.

MTSF significantly influenced only the microbial biomass, positively affecting both Cmicr and Nmicr. These results are partially in agreement with Magnani et al. (2017), who found a positive correlation between MTSF and Cmicr in the same area. The results are consistent with the high levels of sensitivity of soil microorganisms to variations in temperatures, especially in high-elevation ecosystems (e.g., Bing et al. 2016). Wang et al. (2014) and Rui et al. (2011) also reported an increased microbial biomass after a soil warming experiment in alpine meadows. Pintaldi et al. (2019) also found that soil temperature significantly affected Cmicr and Nmicr in a subalpine forest.

Surprisingly, SCD significantly influenced only extractable DOC. Thus, our results differed from those reported by previous works in the same study area (e.g., Magnani et al. 2017; Freppaz et al. 2019), which highlighted the strong influence of SCD on several C and N forms, although under different vegetation covers. As shown by several authors (e.g., Fisk, Schmidt, and Seastedt 1998; D. A. Walker 2000; Edwards, Scalenghe, and Freppaz 2007), in both alpine and arctic tundra, SCD affects the ecosystem functions and structure; however, the intensity of its effects is still unclear, because it varies among different ecosystems. Concerning DOC, our results agree with those found by Freppaz et al. (2019), who indicated a negative correlation between SCD and DOC. The soil temperature usually remained close to 0°C, which promoted the subnival degradation processes, allowing the progressive depletion of organic substrates caused by microbial respiration (Lipson, Schmidt, and Monson 2000). Based on the conceptual model proposed by P. D. Brooks and Williams (1999), our research area could fall in a transition zone between Zone II and Zone III, where small changes in SCD could have remarkable effects on the soil N and C dynamics (Magnani et al. 2017). However, our results indicate a less pronounced effect of SCD, suggesting that other pedoclimatic variables had a first-order effect on C and N forms (e.g., GWC and MOD).

Phenology dynamics and vegetation variables: Plants–soil interaction

Our observations confirmed soil water content and melt-out day as main factors influencing phenology in alpine ecosystems (Ernakovich et al. 2014). More specifically, alpine plants typically depend on snowmelt water (Fisk, Schmidt, and Seastedt 1998), which may influence growth and growing season (Berdanier and Klein 2011). Thus, reduced water rather than nutrient availability is more likely to alter the timing of plant life cycles and productivity (Ernakovich et al. 2014). Furthermore, our results are in agreement with those reported by other research (Prevéy et al. 2019; Quaglia et al. 2020), which indicated snowmelt day as one of the main factors influencing phenophases of snowbed communities. Moreover, the potential change in snowmelt timing and snow thickness as a consequence of climate change would alter season length and nutrient availability (Semenchuk, Elberling, and Cooper 2013; Semenchuk et al. 2015) and produce several long-term plant species-specific effects, affecting future plant composition through possible cascading competition.

Concerning the biotic drivers of soil C and N cycling, the phenology of S. herbacea significantly affected N-NH₄⁺, with a negative relation, whereas it influenced DOC with a positive relation. This relationship could be related to the plant uptake and development: the plant uptake increases with increasing phenophases until senescence, reducing N-NH₄⁺ in soils, whereas the DOC content increases as a result of the release of root exudates into soils (Hütsch, Augustin, and Merbach 2002). In fact, up to 20 percent of the photosynthetically fixed C is released by roots into the soil during the vegetation period as C-rich rhizodepositions, which play several functions in both plant nutrition and soil ecology (Hütsch, Augustin, and Merbach 2002). Regarding N-NH₄⁺, the i norganic nitrogen input from melting snowpack represents a large part of the annual atmospheric nitrogen input in snowbed plant communities (e.g., Bowman 1992; Björk and Molau 2007). Furthermore, N-NH₄⁺ is considered the dominant inorganic N form during the early stage of snowmelt (Mullen, Schmidt, and Jaeger 1998) and, given that nitrates are more leachable, ammonium represents the main inorganic form in soils of alpine and tundra cenoses (Makarov et al. 2010). Ammonium is an important nitrogen source in cold wet climates with acidic soils, and where ammonium is the primary source of nitrogen, also associated with organic N, this may regulate growth (Björk et al. 2007; Nord and Lynch 2009). Our results are in agreement with those reported by Makarov et al. (2010), who observed a gradual reduction of N-NH₄⁺ during the growing season, caused by plant and microorganism uptake. This means that an advanced phenology of plants determined by climate change could lead to an early deployment of available N forms, influencing microbial activity and, thus, decomposition rate and soil fertility.

When considering Nm, its negative relation with DON and Nmicr highlights the vegetation influence on soil N dynamics. Greater values of the parameter Nm correspond in vegetation surveys to a considerable occurrence of plant individuals with high N demands, especially the grass P. alpina, which is one of the most frequent species in these snowbed communities and, to a lesser extent, the forb Taraxacum alpinum (Ten.) DC., which similarly is a eutrophic species but is less common. Therefore, our results verify that the content of microbial and organic N in the soil is reduced with an increase in eutrophic species in the vegetation composition. Our findings align with many literature sources in confirming the valuable contribution of EIVs for the description of indirect properties (Delgado and Ederra 2013; Szymura et al. 2014; Perotti et al. 2018). Even though we were not able to confirm the accuracy of more indexes, we argue that the structure of plant composition and the derived nitrogen EIV are relevant for inferring the nutrient content in snowbed soils and its seasonal fluctuations.

Conclusion

In this study, through four-year monitoring of soil C and N forms and phenology of seven homogenous snowbed areas located at the LTER Site Angelo Mosso Scientific Institute, we investigated the effect of climatic, pedoclimatic, and vegetation variables on soil C and N forms. Our results indicated that in these high-elevation microhabitats, among the climatic and pedoclimatic variables, soil water content, snow melt-out day, and duration of soil freezing during the snow-covered season are the main abiotic factors driving soil C and N dynamics. More specifically, a greater soil water content caused an increase in N-NH₄⁺, N-NO₃⁻, DOC, TDN, Nmicr, Cmicr, and C:Nmicr, whereas an earlier melt-out day decreased soil N-NH₄⁺, N-NO₃⁻, DOC, TDN, Nmicr, and Cmicr. In addition, the duration of soil freezing in the previous winter played an important role, causing an increase in DON and TDN. Among the vegetation variables, the phenology of S. herbacea represented an important biotic driver, suggesting that the plant–soil interaction played a crucial role, significantly influencing the C and N forms in these high-elevation soils. The identified abiotic (snow melt-out day, soil temperature, soil water content) and biotic (plant phenological stages) parameters could be suggested as standard observation variables for the evaluation of the response of the soil biogeochemical cycles in snowbed communities to changes in climatic conditions. These habitats are considered more vulnerable to climate warming than other subalpine and alpine ecosystems and should be carefully monitored in a long-term perspective.

Acknowledgments

Thanks to the Comando Truppe Alpine-Servizio Meteomont for the data from the Automatic Weather Station at Col d’Olen and to Monterosa 2000 SpA. for the logistic support.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary material

Supplemental material for this article can be accessed on the publisher’s website.

Additional information

Funding

This work was supported by the Monterosa 2000 SpA.