Successional Herbaceous Species Affect Soil Processes in a High-Elevation Alpine Proglacial Chronosequence

The study investigated plant-soil interactions along a proglacial chronosequence in the Italian Alps, with a specific focus on pioneer and grassland species structure and biogeochemical processes, with the aim to evaluate the biotic patterns in ecosystem development. We recorded vascular plant frequencies and the mean diameter of one pioneer and one grassland target species in 18 permanent plots distributed along six different stages encompassing a 170-years chronosequence in the Lauson Glacier forefield (NW Italy). We evaluated the main soil properties and measured the C:N:P stoichiometry in the biomass of pioneer and grassland target species and in the underlying soil. For comparative purposes, we analyzed also bare soils sampled near the sampled plant individuals. Pioneer species number and cover significantly increased 10 and 40 years after deglaciation respectively, while alpine grassland species cover and number peaked only after 65 and 140 years, respectively. Along the chronosequence, soils beneath vascular plants were enriched in nutrients, especially under individuals of alpine grassland species, with total organic C contents ranging between 1.3 and 8.9 g·kg−1 compared to 0.2 and 3.3 g·kg−1 in bare soils. Nitrogen content in bare soils was nearly undetectable, while it increased in the plant-affected soils, leading to a more balanced C:N:P stoichiometry in the oldest stages. The colonization of alpine grassland species started immediately, although species number and cover increased only when the soil acquired sufficient nutrient supply and functionality. Although the ecosystem remained C and N limited, the soil could provide adequate conditions for more competitive species establishment, as confirmed by the increasing number and cover of alpine grassland species. Thus, soil nutrient dynamics were strongly influenced by plants, with a major influence triggered by late-successional grassland species.


INTRODUCTION
Throughout the last 170 years, the alpine region has been affected by an enhanced glacier retreat, with the exception of limited periods of glacier advance (Dyurgerov and Meier, 2000;Pörtner et al., 2019). This phenomenon accelerated at historically unprecedented rates during the last few decades, with alpine glaciers that have already lost nearly 50% of their total surface area since 1850 (Paul et al., 2004;Zemp et al., 2006). Consequently, new ice-free surfaces can nowadays provide the opportunity to investigate primary successions, soil development and emerging patterns of ecosystem structure and functions (Matthews, 1992;Egli et al., 2010).
The glacier retreat modifies the hydrological response and the sediment transport within glacier-covered basins, while the newice-free areas are destabilized by paraglacial processes. Mineral substrate left free from ice lacks in biological materials and the net primary production is limited by extreme abiotic constraints, such as low temperatures, intense solar irradiation and strong temperature variations (Stöcklin and Bäumler, 1996;Vitousek et al., 1997;Jones and del Moral, 2009). Receding glaciers open up surfaces for the succession of biota, which protects the released detritus. In turn, the stabilized till supports storage and interception of water and elements that allow plant colonization and soil development (Matthews, 1992;Matthews, 1999;Walker and del Moral, 2003). The combination of these abiotic and biotic processes and their interactions change with time since deglaciation, creating unique chronosequences (Egli et al., 2001;Cannone et al., 2008;Eichel, 2019), characterized by a short space-for-time substitution, from young, recently deglaciated debris to the oldest moraines (Pickett, 1989). Time since deglaciation is indeed assumed to be the key driver for ecosystem development, even though local processes at a small scale are not always negligible (Burga et al., 2010;Dümig et al., 2011).
One of the typical vegetation dynamics in these environments leads to species replacement along successional stages, from pioneer to early, mid-and late-successional species, where the latter are typical alpine grassland species, above treeline (Eichel, 2019). In order to identify successional trajectories and biodiversity dynamics, the correct attribution of the ecological needs to the species is essential, based on phytosociological findings or species traits. Grouping species in functional pools and evaluating their patterns along a chronosequence can provide a more comprehensive interpretation of vegetation dynamics and an easier tool for comparative community ecology (Caccianiga and Andreis, 2004). Additionally, the study of species traits, such as demographic structure of target species belonging to different successional stages (i.e., pioneer or late-successional species) or Competitor, Stress-tolerator, Ruderal (CSR) strategies (Grime, 2006), can allow deep insights into the reproductive mechanisms or the community interactions such as facilitation (Caccianiga et al., 2006;Těšitel et al., 2014).
The progressive plant colonization of high-altitude environments leads to soil formation, including mineral weathering Bernasconi et al., 2011), changes in microbial communities (Nemergut et al., 2007;Jumpponen et al., 2012), carbon (C) accumulation and nutrient cycling, especially nitrogen (N) and phosphorus (P) (Egli et al., 2006;Bernasconi et al., 2011). The combination of pedogenic processes with the C:N:P stoichiometry in the soilmicroorganisms-plant system can furnish key information for understanding the ecological relationships between plant communities and nutritional requirements. In particular, pioneer communities can play an essential role in the kick-off of soil biogeochemical processes, while late-successional ones exert a stronger influence on nutrient dynamics (Matthews 1992). It should be further noticed that, in high altitude primary successions, the effect of single herbaceous species may be easily identified, because of the low and discontinuous vegetation cover. Bonanomi et al. (2016), for instance, described the effect of single plants of Silene acaulis on soil C and N contents along an altitudinal gradient, proving its beneficial effects compared to soil without vegetation.
Soil nutrient dynamics and ecosystem development may in turn drive the C:N:P stoichiometry in photosynthetic tissues (Sardans et al., 2012;Zhang et al., 2018). If C and N generally remain the limiting factors in alpine ecosystems (Körner, 2003), in proglacial environments, the short spatial-temporal span can lead to faster changes in C:N:P patterns in plant species, although nutrient limitation and co-limitation in these environments are still poorly understood. However, a comprehensive approach which encompasses functional species pool cover and successional single species dynamics with soil biogeochemistry still remains largely unknown in high-altitude proglacial chronosequences.
Based on these considerations, we hypothesized that the succession and the related actions of plant species occurring in a short spatial-temporal span since deglaciation play a pivotal role for soil formation and development, leading to favourable ecological conditions for progressively hosting more competitive species. Thus, this study aimed at evaluating pioneer and late-successional grassland species patterns and their role in soil development and nutrient biogeochemistry along a high-altitude proglacial chronosequence in the western Italian Alps (Lauson glacier). To deeply understand the soil-plant interplayed relationship, we further investigated the role of two herbaceous single species in nutrient dynamics, i.e., the pioneer Saxifraga oppositifolia L. subsp. glandulifera Vacc. (hereafter Saxifraga) and the late-successional S. acaulis (L.) Jacq subsp. bryoides (Jord.) Nyman (hereafter Silene).

Study Area
The study was carried out in the proglacial foreland of the Lauson glacier, a small glacier with a 0.23 km 2 surface and 700 m length (Smiraglia and Diolaiuti, 2015; Figure 1), located within Gran Paradiso National Park (SAC/SPA IT1201000) in the upper Cogne Valley (Aosta Valley Region, North-Western Italy). The considered proglacial foreland was between the glacier terminus in 2016 and the frontal moraines left around 1820 during the Little Ice Age, located at 3,030 and 2,750 m a.s.l., respectively. The linear distance between the two moraines is 1,300 m, with a surface of about 0.65 km 2 . The climate of the Cogne Valley is endalpic, with mean annual precipitation of about 700 mm and mean annual temperature of +4.1°C (mean values of Valnontey and Lillaz weather stations, at 1700 m a.s.l.). The proglacial foreland is usually covered by snow between October and mid-July, and it has mean annual temperatures between −1 and −3°C (Mercalli and Berro, 2003). The soils are Eutric Skeletic Regosols (IUSS Working Group WRB, 2015), characterized by homogeneous texture and parent rock composition (mainly paragneiss with small quantities of amphibolites, Le Bayon and Ballevre, 2006), with disturbances limited to cryoturbation, solifluction and a weak water erosion. The vegetation is mainly composed of sparse individuals of alpine pioneer and grassland species, such as S. oppositifolia, Artemisia genipi Weber ex Stechm., S. acaulis and Poa alpina L. Out of the proglacial area, the climax vegetation is alpine grassland, dominated by Carex curvula All. The potential treeline in the area is around 2,350 m a.s.l. (Pecher et al., 2011). Climax soils under the climax alpine prairie are Skeletic Umbrisols or Skeletic Dystric Cambisols (IUSS Working Group WRB, 2015;D'Amico et al., 2020a).
We retraced the retreat of the proglacial terminus since 1820 until 2016 by analyzing historical field surveys carried out by the Istituto Geografico Militare in 1820, 1882, and 1931 and by performing photointerpretation of Regional Technical Maps (1975), ortophotos (1999, and SPOT satellite images (2009) (SctGeoVdA, 2020). A six stage chronosequence was then identified, encompassing a temporal range of 170 years and a spatial distance of 1,000 m (Figure 1), with the lowest limit at 2,850 m a.s.l. (Table 1), ca. 500 m above potential treeline (according to Pecher et al., 2011). Half of the stages were located within 200 m from the 2016 glacial terminus in order to deeply examine the early phases of the primary succession, i.e., the ones occurred in the last 40 years.

Vegetation Surveys
Three permanent quadrat plots (5 × 5 m) were placed at each chronosequence stage, avoiding areas visibly disturbed by water erosion or deposition processes. Plant species composition was then described following the vertical point-quadrat method (Daget and Poissonet, 1971) at every node of a 25 × 25 cm  grid (accounting for a total of 441 sampling points). At each point, all vascular plants were identified at species or subspecies level, while cryptogams were pooled in an aggregated group. A complete list of all occasional species occurring within the plot but not found at the nodes was recorded as well, following the same approach as the phytosociological one (Braun-Blanquet, 1932). The species cover was calculated as its frequency of occurrence across the plot divided by the total number of sampling points. A percentage cover of 0.1% was attributed to all occasional species not found at the nodes but within the plot (Tasser and Tappeiner, 2005). Taxonomic nomenclature followed the new checklist of the Italian vascular flora (Bartolucci et al., 2018). Plant diversity of vascular plants was assessed for each plot in terms of species richness and Shannon-Wiener index (Magurran, 1988). Moreover, we associated a phytosociological optimum (at class level) to each vascular plant species according to Aeschimann et al. (2004) in order to identify different functional pools of species, which are characterized by similar ecological needs. We identified two functional species pools, corresponding to different vegetation successional stages: 1) pioneer species (belonging to Thlaspietea rotundifolii class) and 2) alpine grassland species (belonging to Caricetea curvulae, Carici rupestris-kobresietea, Elyno-seslerietea, Molinio-arrhenatheretea, Nardetea strictae, and Salicetea herbaceae classes). Following this approach, adopted in many different alpine contexts (e.g., Pittarello et al., 2016;Moris et al., 2017;Perotti et al., 2018), we computed species number and cover for each functional species pool in the different plots. Additionally, we studied the population structure of two target species, Saxifraga and Silene, chosen according to the criteria of abundance and ubiquity along the chronosequence. Saxifraga was selected as the pioneer species, with a phytosociological optimum corresponding to Thlaspietea rotundifolii (Aeschimann et al., 2004). Silene was chosen as an alpine grassland species, with a phytosociological optimum corresponding to Caricetea curvulae (Aeschimann et al., 2004). For both target species, we determined the diameter (mean value of the 2 min and max perpendicular diameters, in cm) of all individuals found in each vegetation plot, as a proxy of individual age (Benedict, 1989). We then assessed the population structure by calculating the frequencies of occurrence (number of individuals in 100 m 2 ) of each 1 cm diameter class.

Soil and Plant Tissue Analyses
A soil profile close to each vegetation plot was opened and described following the FAO guidelines (FAO, 2006). Soil samples were then collected from each genetic horizon, airdried and sieved at 2 mm. Soil pH was measured potentiometrically in deionized water (soil:water 1:2.5). Total carbon (TC) and nitrogen (TN) were determined by dry combustion with an elemental analyzer (NA2100, CE Instruments, Rodano, Italia). Carbonate content was assessed by volumetric analysis. The content of total organic carbon (TOC) was determined by difference between TC and carbonate-C.
To evaluate the influence of vegetation on soil properties, in each plot three replicates of Saxifraga and Silene individuals with the largest diameters were sampled, as well as the underlying soil at 0-10 cm depth. For comparative purposes, three soil replicates without vegetation, close to the eradicated plants, were also collected at each stage.
The eradicated individuals of Saxifraga and Silene were placed in sealed polyethylene bags, immediately stored at 4°C in a portable refrigerator, and transported to the laboratory, where green stems and leaves were separated. Aboveground green tissues (i.e., stems and leaves) were dried, grounded to 0.5 mm and analyzed to assess total C, N, and P contents following the above-mentioned procedures.
Soil samples were separated in two aliquots: one was air-dried, sieved, and analyzed to determine TOC, TN, TP, and P av as described above. The other aliquot was processed within 24 h from sample collection, in order to measure microbial biomass C (C micr ) and N (N micr ). In particular, a 30 g aliquot was extracted with 100 ml K 2 SO 4 0.5 M, while a 10 g aliquot was fumigated with chloroform for 18 h before extraction with 50 ml K 2 SO 4 0.5 M. The DOC concentration in non-fumigated soil extracts was determined with a TOC analyser (Elementar, Vario TOC, Hanau, Germany) after filtration with 0.45 μm membrane filters). The same procedure was applied to fumigated samples; the difference in DOC between fumigated and non-fumigated samples, corrected by a recovery factor of 0.45, was C micr (Brookes et al., 1985).
In order to measure N micr , ammonium (extractable N-NH4 + ) concentrations in soil extracts were determined spectrophotometrically (U-2000, Hitachi, Tokyo, Japan) by the method described by Crooke and Simpson (1971). Nitrate (extractable N-NO 3 − ) was determined following Cucu et al. (2014). Total dissolved nitrogen (extractable TDN) in the extracts was determined as reported for DOC. Dissolved organic nitrogen (extractable DON) was determined as the difference between extractable TDN and inorganic nitrogen (extractable N-NH 4 + + N-NO 3 − ). Microbial nitrogen (N micr ) was then calculated as the difference in extractable TDN between fumigated and non-fumigated samples, corrected by a recovery factor of 0.54 (Brookes et al., 1985). All analyses were carried out in triplicate.

Statistical Analyses
The differences among the six stages of the chronosequence in terms of the considered vegetation variables (i.e., cryptogam and vascular plant cover, the two plant diversity indices, and species Frontiers in Environmental Science | www.frontiersin.org January 2021 | Volume 8 | Article 615499 number and cover for the two functional species pools) were tested performing a one-way ANOVA. To investigate the variations within the population structures of Saxifraga and Silene along the chronosequence we performed a PERMANOVA for each species on the frequencies of diameter classes among stages. The similarity index was calculated using Euclidean distance method and 9,999 permutations were set. The soil property changes triggered by Saxifraga and Silene were analyzed with two separate one-way ANOVAs. The same analysis was carried out on chemical properties (i.e., TOC, TN, TP, C:N, and N:P) of Saxifraga and Silene tissues, among stages and between the two plant species.
For univariate tests, assumptions of normality and homoscedasticity were checked with Shapiro-Wilk's and Levene's tests, respectively and, in case of assumption violation, logarithmic and square-root transformations were applied to variables prior to perform the ANOVAs. Whenever normal distribution or homoscedasticity did not occur, even after transformations, we performed the non-parametric Kruskal-Wallis test. Tukey's and Dunn's post-hoc tests were adopted for one-way ANOVA and Kruskal-Wallis, respectively, in case of significant differences.
In order to describe the interactions between soil properties and vegetation cover along the chronosequence, two matrices were arranged: 1) a soil property matrix, with TOC, TN, TP, P av , C micr , and N micr of soil samples without vegetation, and 2) a vegetation matrix, including the covers of cryptogams, total vascular plants, pioneer species, and alpine grassland species. A Mantel test was used to calculate the correlation between the soil and vegetation matrices. A preliminary detrended crosscorrelation analysis (DCCA) was performed to assess the lengths of gradients (Ter Braak and Smilauer, 1998). The DCCA revealed the presence of short (linear) gradients (<4 standard deviations) so we performed a redundancy analysis (RDA) between the two matrices (Ter Braak and Smilauer, 1998).

Vegetation
A total of 65 vascular plants, belonging to 19 botanical families, were recorded within the 18 plots. The most abundant species were Saxifraga bryoides L., Cerastium uniflorum Clairv., Silene acaulis subsp. bryoides, Poa alpina L., and Saxifraga oppositifolia subsp. glandulifera. No species belonging to the Fabaceae family were found in the study area. Five years since deglaciation, seven species colonized the siliceous parent material of the Lauson foreland and the pioneer Saxifraga oppositifolia was the most widespread in terms of both frequency and cover. Vegetation cover increased significantly along the chronosequence, for both cryptogams and vascular plants (Figures 2A,B). In particular, vascular plant cover increased from values close to 0 at 5-10 years since deglaciation to 24% at the 170 years stage (Figure 2b). Species richness and Shannon-Wiener index increased along the chronosequence as well ( Figures 2C,D). The number of species belonging to both functional species pools increased along the chronosequence as well ( Figures 3A,B). In particular, pioneer species number significantly increased already after 10 years and, even if not statistically significant, tended to decline at older stages, while alpine grassland species number showed a slower but progressive increase along the whole chronosequence. Species cover significantly differed among stages for both pioneer and alpine grassland species, reaching the highest values since 40 and 65 years after deglaciation, respectively ( Figures 3C,D).
Overall, we studied 1931 plant individuals for the population structure analysis, 775 were Saxifraga (observed in 14 plots out of 18) and 1,176 Silene (observed in 12 plots out of 18).
In the 5 years deglaciated plots, many Saxifraga individuals (ca. 70 in 100 m 2 ) were counted, with the highest frequencies in the 1 cm diameter class ( Figure 4A). The highest number of individuals and cover was found in plots deglaciated for 40 years, with plants occurring in all diameter classes up to 25 cm and abundance of large individuals (mean diameter 13 cm and several individuals with diameter over 30 cm). This species almost disappeared in plots deglaciated for longer times. Conversely, no individuals of Silene were found within the first deglaciation stage, while it sparsely colonized young moraines only after 10 years ( Figure 4B). The species became more frequent after 40 years (mainly with small individuals) and spread on the 65 years-old areas, with diameter classes up to 19 cm constantly present. In the last two stages, the frequencies of regenerating individuals (1 and 2 cm classes) doubled and large individuals (45 cm class) were detected.
According to PERMANOVA results, the distribution frequency of individual diameters differed among the chronosequence stages for both target species (Figure 4).

Soil Development
According to the WRB classification (IUSS Working Group WRB, 2015), the soils of the first and third stages of the chronosequence were classified as Eutric Skeletic Regosols, while the others as Eutric Skeletic Regosols (Turbic). Soil chemical properties are shown in Table 2. We observed two to five genetic horizons in the soil pits of each stage, including A, AC, CA, and C horizons, with soil depths (A + AC or CA horizons) ranging from 15 to over 40 cm. Soil pH values ranged from 8.9 (C horizon deglaciated for 10 years) to around 6.0 in the upper horizon of the 170 years old stage. pH values decreased below 7.0 only in the surface horizons of the 65 years-old stage. Despite the high pH values, no carbonates were detected. Total organic carbon (TOC) was always greater in surface than in deep soil horizons and increased along the chronosequence, with a maximum content of 26.2 g·kg −1 in the A1 horizon deglaciated for 170 years. Organic horizons were never detected. Total nitrogen (TN) was below the limit of quantification in soils up to 40 years since deglaciation and in all C horizons of the chronosequence (except the 140 years-old one), while the greatest contents were detected in A horizons, ranging between 0.2 and 1.5 g·kg −1 . As a consequence, the TOC: TN ratio varied between 8 (at 40 years) and 17 (at 170 years).
Younger soils (i.e., up to 40 years) showed higher Fe o concentrations than older ones and both Fe o and Fe d were higher in surface horizons than in deeper ones. The Fe o :Fe d ratio, which indicates the ratio between poorly crystalline Fe (hydr)oxides and pedogenetic Fe (hydr)oxides, decreased along the chronosequence with the lowest value in the soils deglaciated for 170 years, in A2 and AC horizons (0.4 g·kg −1 ). Total phosphorus (TP) ranged between 593 and 776 mg·kg −1 without recognizable age or depth trends, whereas the available form (P av ) was 0.5 mg·kg −1 in the youngest soil stages and reached 4 mg·kg −1 in the A horizons of the 170 years-old stage. Cation-exchange capacity (CEC) ranged between 0.89 cmol (+) ·kg −1 in the C horizon of the 65 years-old stage and 10.1 cmol (+) ·kg −1 in the oldest topsoil (Table 3). Base saturation (SATB) was greater in younger stages, ranging between 67% and 100% in the first 40 years-old stages. In the oldest stages, it was generally lower (except for A horizon from the 65 years-old stage with 83%) with a minimum value of 23.7% in the AC horizon of the 140 years ice-free stage. Calcium ion was the main element saturating the exchange complex, both in terms of concentration and ratio, especially in deeper soil horizons ( Table 3).

Soil and Plant Interaction
Chemical composition in the topsoil layers sampled under Saxifraga (SaxS), Silene (SilS) and in the corresponding nonvegetated soils (BareS) significantly differed within the same age and among stages ( Figure 5). The comparisons among all the three target soils were possible only at the third and fourth stages, where we found individuals of both Saxifraga and Silene with adequate development. Total OC showed a lag time in the first stages and then increased, especially with Silene establishment, reaching 6-8 g·C·kg −1 after 140 years, whereas the C content in the non-vegetated soil reached no more than 4 g·C·kg −1 ( Figure 5A). Total N followed the same trend with more pronounced differences between the vegetated and non-vegetated soil samples ( Figure 5B). At the 170 years-old stage SilS showed the highest N content. The stoichiometric TOC:TN ratio did not differ among target soils nor among stages, except for a higher value of SilS compared to BareS at the fifth stage ( Figure 5C). Microbial C and N in vegetated and non-vegetated soils was comparable in the youngest stages and then significantly increased along the chronosequence in SaxS and SilS ( Figures  5D,E). Total phosphorus content did not show significant trends in the three series of soils, while an increasing nutrient availability along the chronosequence, due to the presence of the two species, was observed ( Figures 5F,G). Lower P av concentrations were found in BareS if compared to SaxS at 10 years since deglaciation and to SilS at 140 and 170 years since deglaciation, respectively. In plant tissues, a decreasing trend occurred in P and N with time for Saxifraga and Silene, respectively, leading to a Frontiers in Environmental Science | www.frontiersin.org January 2021 | Volume 8 | Article 615499 corresponding significant N:P and C:N increase ( Figure 6). No differences were found between the two target species for any of the considered chemical features at the third and fourth stages, i.e., where both species were found. A significant correlation was detected between the soil and vegetation matrices by Mantel test (r 0.45, p < 0.001), highlighting that chronosequence stages with similar soil properties had similar vegetation covers. Significant correlations between soil and vegetation variables were observed (Figure 7), explaining 74.2% of the distribution fitting with the first axis and 2.8% with the second axis.
Cryptogam cover appeared weakly associated with soil properties, as highlighted by the arrangement of the arrow relatively to the soil ones. Pioneer species cover showed an intermediate trend between cryptogam and vascular plant covers, this latter being strongly related to TN, C micr , and N micr . Alpine grassland species cover was, instead, strongly associated with TOC and P av . The only variable arranged in the left part of the bi-plot was TP, showing a contrasting distribution compared to most of the other variables. The chronosequence stages separated well along the two axes of the bi-plot. More specifically, the first two stages were at the

DISCUSSION
The univariate and multivariate analyses highlighted a clear pattern of progressive soil development and vegetation successional dynamics throughout the six stages of the Lauson  proglacial foreland, at increasing ages from deglaciation. More specifically, the colonization patterns along the chronosequence of the foreland by vegetation, considering the large difference in elevation between the lower limit of the study area and the potential treeline (500 m), was extremely rapid in the first forty years following deglaciation, in terms of both species richness and cover. An average of four and twelve different species were able to colonize the Lauson proglacial till after 5 and 10 years since deglaciation, respectively. This early process was more rapid than the one observed in the Austrian Alps by Tscherko et al. (2005) and Raffl et al. (2006), in study areas located around the potential treeline (study areas at 2,300-2,450 m a.s.l., with a potential treeline at 2,400 m a.s.l. according to Pecher et al., 2011). Conversely, the studied area was comparable to the results obtained by Burga et al. (2010) in Switzerland, in a proglacial area located at 1900-2,100 m a.s.l., i.e., completely below the potential treeline, which is located at 2,450 m a.s.l. (Pecher et al., 2011). Similar results were also obtained by Cannone et al. (2008) in the Italian Alps, in a study area located at 2,700-2,850 m a.s.l., 300-450 m above the potential tree line, located at 2,400 m a.s.l. (Pecher et al., 2011), and by Matthews and Vater (2015) in southern Norway (study area at 1,420 m a.s.l., 400 m above the local tree line). In all these environments, substrate stability, together with the ability of seeds to reach colonization sites, created a suitable environment for the establishment of pioneer species. The lack of seeds, stated by several studies as a fundamental limiting factor (Chapin et al., 1994;Jones and del Moral, 2009), probably was not the major constraint for plant colonization in the young Lauson foreland, likely because of the short distance between bare sites, freshly released by the glacier, and plant communities growing on nearby slopes, which acted as seed sources. We suggest that, under the studied conditions, other factors, such as drought, frost and cryoturbation may have played a major role limiting plant establishment, according to Erschbamer and Caccianiga (2016). These factors are naturally interrelated with proglacial micro-topography and with the occurrence of protected microsites, important also for nutrient accumulation. In the Lauson chronosequence, after 65 years since deglaciation, vegetation cover increased due to the additive contribution of the alpine grassland species that partially replaced pioneer ones, contributing to an overall rise in vascular plant cover. Despite the FIGURE 6 | Chemical properties of plant tissues for the target species along the six stages of the chronosequence. SaxT, Saxifraga oppositifolia L. subsp. glandulifera Vacc. L. tissues; SilT, Silene acaulis (L.) Jacq subsp. bryoides (Jord.) Nyman tissues; C, carbon (A); N, nitrogen (B); C:N, carbon to nitrogen ratio (C); P, phosphorus (D); N:P, nitrogen to phosphorus ratio (E). Different letters indicate significant differences (p < 0.05) among stages within each target species (uppercase letters) and between target species within each stage (lowercase letters) according to Tukey's HSD test. ns, p ≥ 0.05. increase in species richness, the vegetation cover was very scarce even after 170 years since deglaciation. At the last stage of the chronosequence, cryptogams and vascular plants indeed occupied only 34% of the surface, while Schumann et al. (2016) recorded mean vegetation cover values of 70% in 123 years-old stages from 16 proglacial forelands in the Alps (located between 2000 and 2,600 m a.s.l.). This difference may be related to the complete absence of trees and shrubs within and around Lauson foreland, which is far above the treeline.
The dynamics observed for pioneer and alpine grassland species were confirmed by the population structure of the two target species, Saxifraga and Silene. In particular, Saxifraga showed a typical pioneer behavior, spreading over the study area since the very first stages of the chronosequence and drastically reducing its cover after 65 years, while the more nutrient-demanding Silene firmly established only after (40-) 65 years of ice-free development, indicating a complementary successional pattern. Nevertheless, both species showed an uneven-aged demographic structure, with higher frequencies in the smallest diameter classes and lower frequencies in the largest ones. Saxifraga and Silene seem indeed to perform ruderal strategies (Grime, 2006), which implies a great amount of seed production, resulting in the abundance of seedlings and showing fast-colonization ability, rather than a stress tolerant behavior, as it can be expected in extreme environments such as glacier forelands. The results obtained by our demographic approach, indicating a huge investment in plant reproduction, align with the findings of many authors in proglacial environments, where ruderal strategies allow for plant colonization in the early successional stages (Tscherko et al., 2005;Caccianiga et al., 2006;Gobbi et al., 2010;Erschbamer and Caccianiga, 2016).
Frequency variations in diameter classes were statistically significant for both species throughout the entire chronosequence even if a complete transformation of the demographic structure, from uneven-to even-aged, was not reached in the Lauson glacier foreland.
A clear pattern of a primary succession was observed in both functional species pool dynamics and target species population structure, with an increase in alpine grassland species number and cover over time, and a simultaneous decline in pioneer species. Overall, considering the limited vegetation cover and the partial plant community substitution among successional stages, the Lauson chronosequence time span (170 years) was not sufficient for the establishment of a stable alpine grassland climax community (highlighted also by the Silene demographic structure), supporting the 500 years-minimum age required for Carex curvula All. climax type stated by Andreis et al. (2001).
Although all soils found in the Lauson glacier foreland referred to the WRB group of Regosols, we observed clear soil-forming processes such as a weak mineral weathering (evidenced by the Fe o and Fe d values) and leaching processes (evidenced by pH and SATB decrease), mostly guided by the soil organic matter buildup with time since deglaciation. The decline observed in pH followed the typical pattern of the proglacial chronosequences (Messer, 1988;Frenot et al., 1998;Anderson et al., 2000); it was however much slower than in proglacial areas located at lower elevations (e.g., D'Amico et al., 2014), in which the overall ecosystem productivity was enhanced by higher temperatures and the larger input of plant litter accelerates soil acidification. In our study, soil pH was very high in the recently exposed sediments. As no carbonates were detected, the high pH was probably related to the initial weathering of silicates, which could lead to a release of base cations, able to increase pH values in very coarse-grained soils. In young ice-free till (ca. 5-10 years from deglaciation), not yet vegetated, these processes can lead to very high pH values because of the lack of colloidal buffer surfaces, as deduced by the very low clay content (<2%) and the related scarce cation exchange capacity (max 2 cmol (+) kg −1 in CA horizon in the 5 years old material) (Celi et al., 2013). Moreover, the influence of seasonal glacial runoff cannot be excluded, possibly bringing high Ca and K loads, regardless of bedrock (Anderson et al., 1997(Anderson et al., , 2000 and thus increasing soil pH. High concentrations of K in glacial melt waters, derived from leaching of interlayer cations from biotite, were observed in many works (Drever and Hurcomb, 1986;Stallard, 1995). However, soil acidification occurred over time, likely induced by the increasing vegetation cover and the related release of carboxylic acids through root exudates and/or acidic compounds formed during organic matter decomposition, as well as by organic acids released by bacteria and fungi . The leaching of base cations, due to the lack of exchanging surfaces, could contribute as well to the pH decrease (Burt and Alexander, 1996;Darmody et al., 2005;Bernasconi et al., 2011).
In addition to climatic stress such as drought, frost and cryoturbation, nutrient limitation plays a fundamental role in controlling plant establishment and emphasizes the role of safesites (sensu Harper et al., 1961), where water fluxes increase FIGURE 7 | RDA ordination bi-plot showing the distribution of the six chronosequence stages (represented by the centroids of the three replicates per stage) and the relationships between soil properties (solid lines) and vegetation covers (dashed lines). TOC, total organic carbon; TN, total nitrogen; C micr , microbial carbon; N micr , microbial nitrogen; TP, total phosphorus; P av , bicarbonate extractable (available) phosphorus. The variance explained by each axes is provided in brackets.
Frontiers in Environmental Science | www.frontiersin.org January 2021 | Volume 8 | Article 615499 moisture and nutrient accumulation, enabling the vegetation to survive. Glacial runoff, which represents an important input of nutrients and base cations, could have facilitated the uptake by the early-successional species such as Saxifraga, despite the very harsh substrate of the recently deglaciated debris (Göransson et al., 2016). In the Lauson foreland, we indeed found an optimal base saturation for plant establishment in the first 40 years since deglaciation, with Ca:Mg:K stoichiometric ratios suitable for plant nutrition. Despite the decreasing pH along the terrain age gradient, cation-exchange capacity (CEC) raised, likely as a consequence of organic matter accumulation rather than formation of secondary minerals having higher CEC compared to the unweathered parent material (Bernasconi et al., 2011). Total nitrogen levels were extremely low and were detectable only in the more organic matter-rich and developed upper horizons. This result highlighted that N was one of the major limiting factors of the Lauson soil-vegetation system due to the lack of N-fixing plants (Körner, 2003;Göransson et al., 2016) and to a possible, moderate contribution of N-fixing free-living microbes (Duc et al., 2009). Göransson et al. (2016), in the wetter Damma proglacial chronosequence in Switzerland, found an exponential increasing trend of TN during the first 137 years of soil development and argued that the main N input was probably due to atmospheric deposition (assumed to be approximately 5-10 kg·ha −1 ·yr −1 ). Although in nearby areas (Indren Glacier, Lys Valley, Monte Rosa Group), the atmospheric N deposition was similar, being 6.6 kg·ha −1 ·yr −1 in the last few years (Colombo et al., 2019), the smaller N inputs in our soils were associated with lower mean yearly precipitation because of the isolated position of the Lauson proglacial area, surrounded by high mountains. Thus, our findings suggest that in the Lauson foreland, similarly to areas with low N-deposition like Glacier Bay in Alaska or Franz Josef Glacier in New Zealand (Chapin et al., 1994;Menge and Hedin, 2009), the role of atmospheric deposition was less important in N cycling. Furthermore, Göransson et al. (2016) found, particularly at early stages, that N inputs exceeded plant uptake, and low P levels limited plant colonization more than low N levels. In the entire Lauson chronosequence, soil TP contents ranged between 600 and 700 mg·kg −1 , which represent a favourable range for P cycling and plant establishment (Yang et al., 2013). With time, the release of protons and carboxylic acids caused the dissolution of P-bearing minerals, slightly contributing to feed the available P pool and then plant uptake (Celi et al., 2013). In the older, more structured, and more functional soils, with increased cover of alpine grassland species, P biocycling became particularly evident. Particularly at the fourth and sixth stages of the chronosequence, the upper soil horizons became enriched in TP content. This depth trend can be attributed to plant root exploration and uplift of P from the deeper soil horizons to the surface, through deposition of plant residues. However, compared to other proglacial chronosequences in the Alps, with similar time span, the pedogenic processes remained at an incipient stage, likely because of the lower biomass inputs (D'Amico et al., 2014): the cold temperatures characterizing the Lauson forefield and the associated processes (i.e. short growing season, cryoturbation, root damages and slow decomposition rates) slow down plant colonization, which in turn reduces TOC accumulation and N availability; the low N availability, in a low-N deposition site, in turn, slows down plant colonization, creating a negative feedback which inhibits ecosystem development.
The present study provided novel information concerning the differing function of pioneer and late-successional herbaceous species in affecting soil biogeochemistry in recently deglaciated debris. Indeed, the comparison among soils sampled under plants (SaxS and SilS) and non-vegetated soils (BareS) highlighted the pivotal role of vegetation in soil development and nutrient mobilization, along the whole chronosequence, from the very early stages to the latest ones. Silene performed a greater soil conditioning than Saxifraga in terms of TOC and TN accumulation and development of microbial biomass even within the early stage time points of soil development where both plant species occurred. This may be due to the different growth forms of the two target species, as the cushion species Silene may exert more enhancing effects on soil formation than the prostrate species Saxifraga (Bonanomi et al., 2016). The drastic decline of Saxifraga after 65 years in the chronosequence, likely induced by the competition with more nutrient-demanding species, suggested a facilitation process for the establishment of other plants, linked to the several microsites enriched in organic matter and nutrients. After 140 years, soil development and functionality are proved by higher TOC contents in BareS, as a result of organic matter translocation and mixing processes, which led to even higher concentrations than those in SaxS after 65 years. Microbial and fungal processes could have contributed to this TOC and N enrichment as well. Nitrogen content increased significantly in BareS with time, as observed also in soil pit horizons. Nonetheless, unlike processes involving organic matter, TN resulted in very low concentrations even in the oldest stages, representing the major limiting factor for more nutrient-demanding plant species. Where measurable, the TOC:TN ratio was between 8 and 17.5, thus in the normal range for high-altitude grassland and tundra soils (e.g., D'Amico et al., 2014D'Amico et al., , 2015. The increase in both microbial C and N associated to the target species (and more generally to vascular species cover) along the chronosequence proved again the strong biotic impact in nutrient accumulation compared to abiotic processes, such as atmospheric deposition. In particular, the enhanced microbial biomass and activity associated to the increased concentration of C and N under the target species may have exerted a positive impact on enzyme activity and therefore on microbial N fixation. Similarly, available P content in BareS was very low, but both target species enhanced its concentration along the chronosequence confirming the priming effect of vegetation on P biocycling.
In the youngest stages (ca. 5-10 years since deglaciation) the observed accumulation of P in Saxifraga tissues was probably a consequence of N limitation: under drastic limiting conditions the accumulation of nutrients in tissues has been demonstrated as a possible plant physiological response (Simon et al., 2017). This may be further confirmed by the quite high C:N ratio in the herbaceous tissues, especially for Silene (Hobbie et al., 1998). In older stages, where organic matter and N became more abundant, Frontiers in Environmental Science | www.frontiersin.org January 2021 | Volume 8 | Article 615499 P content in Saxifraga tissues fell to more typical values (Richardson et al., 2004). However, N:P in plant tissues was very low, along the whole chronosequence, and far from the value of 12, a threshold beyond which plant development is P-limited, such as in highly developed soils (Wardle, 2004) or in particular P-limited systems. For instance, in the Damma proglacial chronosequence (N:P ratios in the aboveground vegetation >19, Göransson et al., 2016) the main P-containing mineral was the slowly weathering fluoroapatite (Bernasconi et al., 2011), while the Verra Grande forefield (D'Amico et al., 2020b) developed on serpentinite and almost devoid of P-bearing minerals.

CONCLUSION
This work provided novel results on individual plant species trajectories and related soil development through successional stages of a proglacial chronosequence. Despite the extreme conditions for plant establishment, the Lauson proglacial foreland allowed a remarkable early (5-10 years) and speciesrich colonization, attained by pioneer species in the available safe micro-sites. Carbon seemed to be the major constraint to microbial development together with N, which also hampered plant colonization. However, nutrient requirements for pioneer plants were met, likely thanks to the water fluxes funnelled between rocks. The alpine grassland species colonization started indeed immediately, although species number and cover pronouncedly increased only when the soil reached sufficient nutrient supply and functionality. Then, after about one century of ecosystem development since deglaciation, soil conditions appear suitable for the establishment of more competitive vegetation, typical of alpine grassland communities. Nonetheless, the improving effect on soil triggered by S. oppositifolia was already evident on quasi-inert substrates, free from ice for a limited time, highlighting the relevant role of biotic factors in stimulating the initial phases of soil formation. However, more biomass-productive species, such as S. acaulis, performed a greater soil conditioning than S. oppositifolia, even within the early stage time points where both plants were present, becoming progressively more expressed with time. Therefore, C:N:P stoichiometry in proglacial areas is plantinfluenced, and the impact on nutrient dynamics is related to functional species pool, with a major role played by latesuccessional grassland species.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: 10.5281/zenodo. 4067891.

FUNDING
The project was funded by Gran Paradiso National Park.