The Consequences of Glacier Retreat Are Uneven Between Plant Species

Losapio, Gianalberto; Cerabolini, Bruno E. L.; Maffioletti, Chiara; Tampucci, Duccio; Gobbi, Mauro; Caccianiga, Marco

doi:10.3389/fevo.2020.616562

ORIGINAL RESEARCH article

Front. Ecol. Evol., 29 January 2021

Sec. Biogeography and Macroecology

Volume 8 - 2020 | https://doi.org/10.3389/fevo.2020.616562

This article is part of the Research Topic Consequences of Climate Change for Plant Biodiversity in High Mountain Ecosystems View all 18 articles

The Consequences of Glacier Retreat Are Uneven Between Plant Species

$\nGianalberto Losapio$ Gianalberto Losapio¹^*

Bruno E. L. Cerabolini²

Chiara Maffioletti³

Duccio Tampucci³

Mauro Gobbi⁴

Marco Caccianiga³

¹Department of Biology, Stanford University, Stanford, CA, United States
²Department of Biotechnologies and Life Sciences, University of Insubria, Varese, Italy
³Department of Biosciences, University of Milan, Milan, Italy
⁴MUSE – Museum of Science, Trento, Italy

Glaciers are retreating worldwide, exposing new terrain to colonization by plants. Recently-deglaciated terrains have been a subject of ecological studies for a long time, as they represent a unique natural model system for examining the effects of global warming associated with glacier retreat on biodiversity and the spatio-temporal dynamic of communities. However, we still have a limited understanding of how physical and biotic factors interactively influence species persistence and community dynamics after glacier retreat and glacier extinction. Using hierarchical joint species distribution models, we integrated data on plant species occurrence at fine spatial scale, spatio-temporal context, environmental conditions, leaf traits, and species-to-species associations in plant communities spanning 0 to c 5,000 years on average after glacier retreat. Our results show that plant diversity initially increases with glacier retreat, but ultimately decreases after glacier extinction. The 22% of plant species non-linearly respond to glacier retreat and will locally disappear with glacier extinction. At the local scale, soil carbon enrichment and reduction of physical (topographic) disturbance positively contribute to distribution patterns in 66% of the species, indicating a strong impact of community-level environmental conditions. Furthermore, positive and negative associations among species play a relevant role (up to 34% of variance) in driving the spatio-temporal dynamic of plant communities. Global warming prompts a shift from facilitation to competition: positive associations prevail among pioneer species, whereas negative associations are relatively more common among late species. This pattern suggests a role of facilitation for enhancing plant diversity in recently ice-free terrains and of competition for decreasing species persistence in late stages. Associated to that, species persisting the most show more “conservative” traits than species of concern. In summary, although plant diversity initially increases with glacier retreat, more than a fifth of plant species are substantially declining and will disappear with glacier extinction. Even for the “winners,” the “victory” is not to be taken for granted due to the negative impact of rising competition. Integrating survey data with hierarchical and network models can help to forecast biodiversity change and anticipate cascading effects of glacier retreat on mountain ecosystems. These effects include the reduction of ecosystem services and benefits to humans, including food production from the pioneer species Artemisia genipi.

1. Introduction

Glaciers are retreating, shrinking, and disappearing worldwide due to global warming (Marzeion et al., 2014; Zemp et al., 2015; Roe et al., 2017). The retreat of glaciers affects landscape configuration, water availability, and natural hazards, impacting downstream ecosystems and society (Fell et al., 2017; Cauvy-Fraunié and Dangles, 2019). Therefore, plants inhabiting proglacial environments are highly sensitive to global warming and inherent glacier retreat (Erschbamer, 2007; Dullinger et al., 2012). Understanding how plant communities are responding to ongoing glacier recession is crucial for projecting biodiversity dynamics and planning sustainable ecosystem management (Walker and del Moral, 2003; Erschbamer and Caccianiga, 2016).

As glaciers retreat, new terrains are exposed to colonization by living organisms. The further away and downstream from the current glacier front, the longer ago the terrain was exposed by glacier retreat, the older the community. Thus, in general, community age increases with increasing the distance from the retreating glacier. With the support of geochronological records, one can reconstruct the past positions of glaciers, thus date the terrain and assign an age to the community (Matthews, 1992). During the last 10,000 years, the Little Ice Age (LIA, ca 1350–1850 CE) was the main episode of re-advance of glaciers in the European Alps (Maslin et al., 2019). Despite the rate of glacier retreat varies regionally, including two short recent periods of glacier advances in the 1920s and 1980s, glaciers are retreating at historically unprecedented rates globally (Zemp et al., 2015). For these reasons, the areas in front of retreating glaciers (hereafter, glacier forelands) are unique natural systems and prominent models for studying how populations change, communities assemble, and ecosystems develop over spacetime in response to global warming (Whittaker, 1993; Chapin et al., 1994; Caccianiga et al., 2006).

Various theories have been proposed to describe the orderly process of ecosystem development involving changes in species composition and communities over time, i.e., succession. The classic Connell and Slatyer models account for facilitation, tolerance, or inhibition between initial colonizers and later ones (Connell and Slatyer, 1977). In the facilitation model, initial colonizers ameliorate environmental constraints on plant colonization, growth and reproduction, facilitating the establishment of later colonizers. In the tolerance model, differential life history traits and the differential ability of species to tolerate initial environmental conditions determine the arrival and growth of species, with later colonizers growing in spite of the presence of initial ones. In the third one, initial colonizers inhibit the growth and establishment of later species. The hierarchical framework proposed by Pickett et al. (1987b) expanded on the importance of initial conditions (site availability) and on the concept of differential species availability and performance (e.g., dispersal, growth rate, life history strategies, environmental stress) to explain mechanisms of succession. The geoecological model (Matthews, 1992; Whittaker, 1993) postulates that the environment exerts more control over community development than organisms in the initial stages, while habitat modification by organisms and species interactions become more important in later development phases. Similarly, complex combinations of negative and positive biotic interactions structure plant communities (Chapin et al., 1994). In particular, facilitation is most common in harsh environments (initial stages), while competition becomes more important in more favorable environments (later stages) (Bertness and Callaway, 1994; Callaway and Lawrence, 1997). Finally, Grime's c-s-r theory (Grime, 2001) built on differences among plants in their competitive ability, tolerance to physical stress, and dependence on colonization over resource and disturbance gradients. It describes a pathway defined by control of disturbance initially that changes toward control of resources by competition and then tolerance to biotic stress (Caccianiga et al., 2006). Recent syntheses (Meiners et al., 2015; Erschbamer and Caccianiga, 2016; Pulsford et al., 2016) highlight that mechanisms underlying changes in biodiversity are often a complex interaction of species traits, physical environment, stochasticity, and effects of neighboring plants. However, the role of ecological networks for shaping biodiversity, driving species distribution and contributing to species response to global warming remains overlooked. Since biodiversity is more than a list of species (Bascompte and Jordano, 2014), understanding the contribution of ecological networks to species' response to global warming is of paramount importance for forecasting biodiversity, developing scenarios, and mitigating species loss.

There is evidence that biodiversity increases with glacier retreat (Raffl et al., 2006; Cauvy-Fraunié and Dangles, 2019), but this holds true only as long as glaciers are still present (Stibal et al., 2020). In fact, we are facing the extinction of glaciers around the globe (Dyurgerov, 2005; Zemp et al., 2015), especially in tropical mountain ranges and mediterranean mountains (Apennine, Sierra Nevada of California), while many glaciers will disappear within the next decade in the Alps, temperate Andes, Caucasus and Himalaya (Huss et al., 2017). Glacier forelands comprise unique habitats that host distinctive organisms and communities. Increasing evidence indicates high turnover rates in mountain environments (Steinbauer et al., 2018), where increase in species richness of lower-altitude plants is followed by decrease of populations and range contraction of high-altitude species (Thuiller et al., 2005; Pauli et al., 2007; Dullinger et al., 2012). Accordingly, it is possible that not all the species would have the possibility to migrate and colonize to suitable habitats nearby (Körner, 2005), so high-altitude plants associated to glacier environments will not persist locally following glacier extinction. However, the consequences of glacier retreat and extinction for plant species inhabiting glacier forelands as well as whether and how ecological networks mediate such impact remain poorly understood and quantified.

Here, we hypothesize that glacier extinction will severely reduce a set of plant species, ultimately causing biodiversity loss. Furthermore, we propose that facilitation will support persistence of initial colonizers while competition will exacerbate extinction risk. We answer the following questions: (i) How do plant species and communities respond to glacier retreat? (ii) Which and in what proportion will species persist or disappear after glacier extinction? (iii) How do plant networks change over spacetime and how they mediate species' response to glacier retreat?

2. Materials and Methods

2.1. Data Collection

We integrated published datasets on plant species distribution and leaf traits with unpublished, original data on environmental conditions. Published datasets on plant community surveys were compiled including identity and abundance of plant species at local spatial scale (meters) at the following sites: (i) Vedretta d'Amola glacier (Adamello-Presanella Alps, Italy, 46°13'16” N, 10°40'41” E), characterized by a debris-covered glacier colonized by plants, glacier tongue position dated to 1994 and 1925, and LIA moraine (Losapio et al., 2016); (ii) Western Trobio glacier foreland (Orobic Alps, Italy, 46°03'36” N, 9°57'36” E), characterized by moraines dated to 1985, 1920, 1900, and LIA (Tampucci et al., 2015); (iii) Rutor glacier foreland (Graian Alps, Italy 45°38'0” N, 7°1'0” E), characterized by moraines dated to 1933, 1873, and LIA (Caccianiga and Andreis, 2004), and (iv) Vedretta di Cedec glacier foreland (Ortler Alps, Italy, 46°27'N, 10°35' E), characterized by three main moraines dated to 1986, 1965, and LIA (Gobbi et al., 2010). Each site comprised grasslands that are adjacent to each glacier foreland and occur in terrains beyond LIA moraines. These communities are dated to the Late Glacial period, from 200 to 10,000 years BP. Given this landscape configuration, we consider: (i) communities occurring in glacier forelands being representative of a scenario of glacier retreat, with glaciers still present and affecting ecosystems locally, and (ii) communities beyond LIA moraines (outside glacier forelands, the “potential natural vegetation”, Caccianiga and Andreis, 2004) being the representative of a scenario of glacier extinction, which are not directly affected by Anthropocene glacier dynamic.

Surveys were replicated three to seven times in 25–100 m² plots in each community of different age along a transect starting from the glacier front. Leaf traits including Specific Leaf Area (SLA) and Leaf Dry Matter Content (LDMC) were downloaded from TRY database (Kattge et al., 2020) and complemented with our own original measurements. SLA and LDMC were measured following standard protocols (Pérez-Harguindeguy et al., 2013). For each species, we sampled from 5 to 15 fully expanded leaves selected randomly from the outer canopy of different adult healthy plants within the same population. Leaf material was stored at 4°C overnight to obtain full turgidity for determination of fresh weight (LFW) and area (LA). LA was determined using a digital leaf area meter. Leaf dry weight (LDW) was then determined following drying for 24 h at 105°C. SLA was calculated as the ratio between LA and LDW. LDMC was calculated as the ratio between LDW and LFW.

In order to harmonize data across different studies for the further analysis of plant diversity dynamic, we considered species occurrence (presence/absence) in each study plot (sample). We then pooled data across sites into a single community matrix Y

\begin{array}{l} Y : [\begin{matrix} y_{1, 1} & y_{1, 2} & \dots & y_{1, j} \\ y_{2, 1} & y_{2, 2} & \dots & y_{2, j} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ y_{i, 1} & y_{i, 2} & \dots & y_{i, j} \end{matrix}], & (1) \end{array}

where plant species j are arranged in columns, plots i in rows, and the entries y_ij represent the species occurrence in the plot, with y_ij equals to 0 or 1 in case of absence or presence of the species, respectively.

Community age x of each plot i was estimated as the average years since glacier retreat between two moraines as $x_{i} = \sqrt{x_{s} - \frac{(x_{o l d} + x_{y o u n g})}{2}}$ , where x_s is the year of sampling, x_old and x_young is the geochronological information (age) of older and younger moraines adjacent to the community, respectively. For instance, communities located between the glacier position in 2013 and the moraine dated 1994 were considered having an average age of 10 years. In the same way, communities surveyed in 2014 between two moraines deposited in 1985 and in 1920 have an estimated average age of 62 years, and plots beyond the LIA moraines have an estimated average age of 5,000 years. This approach allows us quantifying biodiversity patterns and dynamics with greater resolution than before by overcoming the limitations of using categorical data. As any other mean estimate, however, it has the limitation of focusing only on average trends excluding variation around such mean. In fact, notice that there is substantial geochronological uncertainty surrounding the age of communities beyond LIA moraines, since their age can span a range of 170–10,000 years.

In light of this potential shortcoming, we then explored the implications of such uncertainty and the robustness of our results against this assumption. We adopted a second method that relies on relative differences in elevation along the glacier foreland transect. In particular, for the seek of among-site comparison and synthesis, we standardized elevation v for each plot i at each site k as $v_{i} = 1 - \frac{v_{i} - min (v_{k})}{max (v_{k}) - min (v_{k})}$ , where the maximum and minimum values of elevation along the glacier retreat transect were taken for each site.

Then, we proceeded with measuring the following soil parameters in the Vedretta d'Amola glacier and Western Trobio glacier: (i) soil organic matter (SOM); (ii) gravel content; and (iii) pH. Soil samples were taken at the surface for physical and chemical analysis. A sample of about 1 kg was taken at every plot for particle size distribution analysis; a sample of about 200 g was taken for soil pH (in 1:2.5 soil:water), calcium carbonate content (Dietrich–Fruhling calcimeter) and organic matter content (Walkley–Black method) (Sleutel et al., 2007).

Spatio-temporal context and soil conditions data were organized into a matrix X

\begin{array}{l} X : [\begin{matrix} x_{1, 1} & x_{1, 2} & \dots & x_{1, m} \\ x_{2, 1} & x_{2, 2} & \dots & x_{2, m} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ x_{i, 1} & x_{i, 2} & \dots & x_{i, m} \end{matrix}], & (2) \end{array}

where m environmental variables (i.e., community age x_i, standardized glacier retreat v_i, and soil conditions) are arranged in columns, plots i in rows. The integrated datesets for the Vedretta d'Amola glacier foreland and Western Trobio glacier foreland comprised 51 plants species and two soil variables (SOM and gravel) over 54 plots and 54 plant species and three soil variables (SOM, gravel, and pH) over 30 plots, respectively.

2.2. Hierarchical Modeling

We analyzed plant diversity and community dynamics using the Hierarchical Modeling of Species Communities (HMSC) (Ovaskainen et al., 2017; Ovaskainen and Abrego, 2020; Tikhonov et al., 2020). This Bayesian framework accounts for the influence of environmental conditions, species' response traits, random variation in species occurrence, and species-to-species associations for modeling the joint distribution of species. In particular, we combined plant species occurrence data (Y matrix with plant species in columns and plots in rows) with the spatio-temporal context (X matrix defining sites and community age in columns for each plot in rows) and trait data (matrix with SLA and LDMC for each species) accounting for the potential effects of environmental conditions (SOM, gravel, and pH) and biotic interactions (positive and negative species associations) (Figure 1).

FIGURE 1

Figure 1. On the left, schematic representation of glacier retreat dynamic, with “pioneer” (dark green), “early” (light green), and “late” (red) stages of ecosystem development. On the right, framework of the hierarchical joint species distribution modeling for forecasting biodiversity change.

The occurrence y of each species j in each plot i was modeled as y_ij ~ B(L_ij), which follows a binomial distribution B, where L_ij is a linear predictor composed of fixed and random parts. Fixed effects F are modeled as the quadratic regression $L_{i j}^{F} \propto \sum_{i = 1}^{S} x_{i} β_{j}$ , where x_i indicates community age at plot i, and β_j the parameters of the second degree polynomial measuring species response j to glacier retreat dynamic across all S plant species. In addition, we ran the same model with standardized elevation v_i in place of x_i. The parameters β_j are estimated following a multivariate normal distribution and are influenced by species-specific traits (SLA and LDMC). Sites were included as random term. This model was fitted to the integrated dataset (Y matrix) comprising 117 plant species and two leaf traits over 170 plots at four glacier-foreland sites.

In a second model, environmental conditions (SOM, gravel, and pH) were included as fixed effects, in addition to community age, and spatial correlation among plots within each site was included as random term. This allowed us inferring signals of “biotic interactions” by identifying species pairs that co-occur more or less often than expected by chance, after controlling for the response of species to environmental conditions (Ovaskainen et al., 2016). This way, associations between species are reflected by significant co-occurrences that cannot be attributed to the responses of species to the environment (Ovaskainen and Abrego, 2020).

HMSC uses a latent variable approach to estimate species associations, which has much fewer parameters to be estimated than other approaches because the number of latent factors is much smaller than the number of species (Ovaskainen et al., 2017; Tikhonov et al., 2020). At the community level, this approach generates patterns of associations between species embedded in a variance–covariance matrix R. Matrix entries R_j1j2 measure the degree to which species j₁ and j₂ co-occur more or less often than expected by chance net of the influences of environmental conditions. Thus, latent variables can represent the outcome of ecological interactions, with factor loading representing the response of species j₁ to that latent interaction with species j₂ (Ovaskainen et al., 2016). Two separate models were fitted for the Vedretta d'Amola site and Western Trobio site.

Parameters are estimated using Gibbs MCMC with a flat prior and 2,000 samples from the posterior distribution across four chains (thin = 2). The MCMC scheme converged as chains reached a stationary distribution. We evaluated model fit in terms of explanatory power. Analyses were conducted using the “Hmsc” package (Tikhonov et al., 2020) in R version 4.0.2 (R Core Team, 2020).

2.3. Network Analysis

To explore changes in facilitation and competition among plant species across glacier forelands, we adopted a multilink network approach (Losapio et al., 2019). Associations between species significantly differrent from zero (95% CI) $R_{j 1 j 2}^{*}$ were then retained to build species association networks. Species association networks were composed of plant species and associations between species taking the value of 1 or −1 if significant $R_{j 1 j 2}^{*}$ were higher or lower than 0, respectively.

Finally, to identify and classify species' response to glacier retreat, we implemented a species distribution network approach (Burns and Zotz, 2010; Losapio et al., 2019; Marini et al., 2019). We first transformed the plant community survey data into a bipartite network of species distribution over the landscape. In such network, plant species and community age of each plot are linked by species occurrence y_ij, representing the two parts of the network. Then, a fast greedy algorithm that optimizes a modularity score (Clauset et al., 2004) was used to analyze community structure and identify modules of plant species that are distributed across community ages. Network modules are dense subnetworks characterized by high occurrence frequency of a group of species within the same community age and low frequency or no occurrence between other species groups at different community age. This allowed us unveiling large-scale patterns of species response to glacier retreat. Two separate networks were built for the Vedretta d'Amola site and the Western Trobio site.

2.4. Statistics

First, we explored parameter estimates for β_j (species responses to glacier retreat), reporting mean and 95% posterior probability, as well as we generated predictions of plant diversity at the community level over glacier forelands as a function of terrain age, from zero (debris-covered glaciers) to 5,000 years.

Second, we estimated the probability of species persistence p_j after glacier extinction, at 5,000-years communities. Then, we related this probability to community type obtained from modularity analysis of species distribution networks. We used linear mixed-effects models including persistence probability as response, community type as predictor, and species identity nested within site as random effects.

Third, we quantified: (i) the frequency of positive and negative associations per species across network modules, and (ii) the density of positive and negative species associations per plot across the glacier retreat gradient in Vedretta d'Amola and Western Trobio sites. For the former, we used zero-inflated generalized linear models including association frequency as response (zero-inflated poisson distribution), association type (positive or negative), community type, and their statistical interaction as predictor. For the latter, linear models were fitted to association density in response to community age (continuous, square-root transformed), community type, and their statistical interaction. Two separate models were fitted for the two sites.

Finally, we analyzed the statistical association between traits and species distribution. We used a linear model for SLA and LDMC including community type (pioneer, early, late) as predictor.

Model output was evaluated in terms of both variance explained and parameter estimates. Analyses were conducted in R version 4.0.2 (R Core Team, 2020), using the packages “lme4” (Bates et al., 2015), “car” (Fox and Weisser, 2019), and “glmmTMB” (Brooks et al., 2019).

3. Results

Glacier retreat affects the 51% of plants in a species-specific way (Supplementary Table 1), among which 29% of species will expand while 22% will disappear with glacier extinction. Examples of species that will flourish the most are: Minuartia verna, Veronica fruticans, Achillea moschata, Trifolium pallescens, Draba aizoides, Poa alpina, Trisetum spicatum, Gentiana nivalis, Myosotis alpestris, Carex curvula, Antennaria dioica, Leontodon helveticus, Potentilla aurea, Senecio incanus. Species of main concern for local extinction are: Saxifraga bryoides, Artemisia genipi, Cardamine resedifolia, Leucanthemopsis alpina, Gnaphalium supinum, Sedum alpestre, Minuartia sedoides, Sempervivum arachnoideum, Hieracium staticifolium, H. glanduliferum. Species that are ubiquitously distributed and do not show significant response to glacier retreat and environmental conditions are: Anthoxanthum alpinum, Gentiana kochiana, G. punctata, Ligusticum mutellina, Pedicularis kerneri, Phyteuma hemisphaericum. Functional traits contribute to explain the 4% of variation in species distribution.

Overall, plant diversity dynamic following glacier retreat is characterized by a hump-shaped trend (Figure 2A). In particular, our predictions show that diversity will increase immediately after glacier retreat with increasing community age x_i (β = 0.017). After saturating, plant diversity will decline toward poorer communities in the late stages of ecosystem development (β = 0.020). Results are qualitatively similar when using standardized differences in elevation v_i (Figure 2B).

FIGURE 2

Figure 2. Predicted effect of glacier retreat (x-axis) on plant diversity (y-axis). (A) Community age was estimated as the mean between two consecutive moraines. Notice that there is substantial geochronological uncertainty surrounding the age of communities beyond LIA moraines, since their age spans a range of 170–10,000 years (here, on average 5,000 years). (B) Community age was estimated as standardized differences in elevation along the glacier foreland transect.

We then considered plant community dynamics in the two sites (Amola and Trobio glacier forelands) with high resolution, local-scale data on environmental conditions. Glacier retreat, soil organic matter, gravel and species associations explain the 66–92%, 10–16%, 14–17%, 1–24% (95% CI) of the variance observed in species distribution along the Amola glacier foreland, respectively. In the Trobio glacier foreland, glacier retreat, soil organic matter, pH, and species associations explain the 61–94%, 1–11%, 1–12%, 1–34% (95% CI) of the variance observed in species distribution, respectively.

In both sites, species distribution networks were structured in three main modules (communities), which were composed by the 29, 21, and 50% of the species pool in the Amola glacier foreland and by the 42, 18, and 40% of the species pool in the Trobio glacier foreland. Results indicate that these communities are characterized by three distinct time periods of glacier retreat: “pioneer” (less than 100 years), “early” (between 100 and 150 years), and “late” (older than 150 years). We found that community type significantly predicted the persistence probability of species (χ² = 117.2, P < 0.001, Figure 3). In particular, pioneer species were only 11.7 ± 8.0% likely to persist after glacier extinction, while early and late species were marginally (P = 0.079) and significantly (P < 0.001) more likely to persist, having a persistence probability of 21.3 ± 8.4% and 56.1 ± 7.8%, respectively. Variance among species and between sites was low (σ = 0.038, σ = 0.011).

FIGURE 3

Figure 3. (A) Distribution of three “model” species along the glacier foreland: the pioneer Artemisia genipi (dark green), the early Luzula alpino-pilosa (light green), and the late Carex curvula (red). (B) Persistence probability of species after glacier extinction across pioneer, early, and late species.

Having shown that species differed in their persistence, we examined the distribution of network-level associations along the two glacier forelands. The Amola network comprised 150 associations among 51 species (connectance = 0.06), among which 78 positive associations and 72 negative. Pioneer, early, and late species differ in their frequency of positive and negative associations (χ² = 5.2, P = 0.075). In particular, late species were significantly more negatively associated to other species than early species (β = 0.79 ± 0.38, P = 0.037; Figure 4A). Furthermore, the density of positive and negative associations depended on glacier retreat dynamic (χ² = 12.8, P < 0.001), as negative associations became relatively denser than positive with increasing years after deglaciation (β = 0.0007 ± 0.0004, P = 0.011; Figure 4B). The Trobio network comprised 90 associations among 54 species (connectance = 0.03), among which 54 positive associations and 36 negative. The frequency of positive and negative associations was significantly different (χ² = 5.9, P = 0.014) and varied among plant communities (χ² = 22.8, P < 0.001, Figure 4C). In fact, positive associations were significantly more prevalent than negative ones overall (β = 1.13 ± 0.29, P < 0.001). Yet, negative associations were significantly more frequent among early and late species than pioneer ones (β = 4.03 ± 1.19, P < 0.001, and β = 4.48 ± 1.24, P < 0.001, respectively). Moreover, the density of associations marginally decreased with glacier retreat (β = −0.002 ± 0.001, P = 0.088; Figure 4D).

FIGURE 4

Figure 4. Distribution of positive (red) and negative (blue) associations across species (A,C) and along the glacier foreland (B,D) in the Amola site (above) and Trobio site (below). Reported marginal means and 95% CI.

Finally, we found that community type can predict plant traits (P = 0.056 for SLA and P < 0.001 for LDMC). In particular, pioneer and early species have different trait values from late species, having significantly higher SLA (β = 3.34 ± 1.60, P = 0.039, Figure 5A) and significantly lower LDMC (β = 7.62 ± 1.90, P < 0.001, Figure 5B).

FIGURE 5

Figure 5. Association between community type (x-axis) and plant traits (y-axis). (A) Specific Leaf Area. (B) Leaf Dry Matter Content.

4. Discussion

We examined the poorly understood consequences of glacier retreat for biodiversity over a pool of 117 plant species in different alpine glacier forelands. Although glacier retreat triggers and favors the development of communities in the short term, plant diversity will decrease with glacier extinction in the long term: about one fourth of plant species would disappear locally once glaciers vane. Soil organic matter and gravel further increase and decrease local diversity by enhancing and diminishing species ranges, respectively. Besides changes in species richness, we also observed changes in the way species are associated—and potentially interact—with each other. Positive associations tend to become less prevalent than negative ones with glacier retreat, potentially decreasing species persistence probability after glacier extinction. Species persisting the most with glacier extinction show more “conservative” leaf traits than species of concern. Taken together, these patterns reveal the species-specific impact of global warming on biodiversity is mediated by complex ecological interactions, highlighting the need of understanding them in a systemic and context-dependent way. Biodiversity conservation and ecosystem management require focusing on both species of concern and their networks of relationships with other species and the environment to support and ensure plant persistence after glacier extinction.

4.1. Plant Diversity Dynamic

Our results indicate that climate warming favors plant diversity as long as glaciers can withstand rising temperature. Model predictions show that plant diversity rapidly increases for c one thousand years after glacier retreat. Late species like Carex curvula, Leontodon helveticus, Potentilla aurea, and Trifolium pallescens are gaining from higher temperatures and their populations will likely flourish also when glaciers disappear. However, climate warming decreases species persistence once glaciers vane, ultimately prompting local extinctions. After a period of saturation and stasis, plant diversity will start declining substantially, losing half of the species across plant communities at the end of the dynamic and facing more than 20% of extinction from the species pool. In fact, not the whole ecosystem will gain species since pioneer and early plants like Artemisia genipi, Saxifraga bryoides, Cardamine resedifolia, and Minuartia sedoides decline already 100 years after glacier retreat. With a probability of persisting lower than 0.2, these species will disappear locally in the absence of glaciers.

Such hump-shaped curve and idiosyncrasy reflect that the effects of climate warming produce “winners” and “losers” (Grabherr et al., 1994; Erschbamer, 2007). The likely reason is that the response of plants to increasing temperature varied among species and communities (Körner, 2005). Indeed, our results are consistent with experimental studies (Erschbamer, 2007), surveys (Walther et al., 2005; Pauli et al., 2007; Steinbauer et al., 2018), and species distribution modeling (Thuiller et al., 2005; Engler et al., 2011; Dullinger et al., 2012) indicating that although species richness increases with increasing temperature in high-altitude environments, there is also range contraction, replacement and biodiversity loss. Furthermore, they support the “replacement model” of succession (Pickett et al., 1987a; Van Andel et al., 1993; Walker and del Moral, 2003; Cutler et al., 2008), according to which changes in species richness and community composition occur as a result of changes in competitive “winners.” Yet, it is worth noting that the persistence probability of late species after glacier extinction is on average 0.6. That is, the “victory” for the “winners” is not to be taken for granted. It is possible that competition among late plants decreases species persistence and therein biodiversity after glacier extinction.

Since ecosystem development following glacier retreat involves also changes in pollination networks (Losapio et al., 2015, 2016), arthropod community (Hodkinson et al., 2001; Kaufmann, 2001; Gobbi et al., 2010), food-webs (Raso et al., 2014; Sint et al., 2019), and soil microbial communities (Sigler et al., 2002; Jumpponen, 2003; Schütte et al., 2009), it is likely that the consequences of changing plant diversity will extend beyond plant communities. For instance, species of concern like Saxifraga bryoides and Leucanthemopsis alpina play a key role in the assembly and structure of ecological networks, supporting pollinators, herbivores, predators, detritivores, and parasitoids (Losapio et al., 2015). Thus, the local extinction of these plants would not only put other species at risk of extinction but would also severely affect the functioning of alpine ecosystems. The loss of biodiversity is going to reduce ecosystem services and benefits to humans, like food production from the pioneer, species of concern Artemisia genipi.

4.2. Plant Association Networks and Traits

After accounting for the response of species to glacier retreat and the effects of environmental conditions and random variation in species occurrence, we found a strong signal of biotic interactions in species distributions. In fact, species-to-species associations contribute to explain a significant proportion of variance, up to 34%, observed in the way plant species tended to be more or less associated to each other than expected by all other factors. Positive and negative species associations characterize plant networks, suggesting that both facilitation and competition may contribute to driving species distribution along glacier forelands and creating the hump-shaped diversity pattern over spacetime.

In particular, facilitation prevails among initial colonizers in pioneer communities, whereas competition was more prevalent in late communities. Overall, competition becomes relatively more important than facilitation with increasing community age. Particularly, this trend results from two distinct and non-mutually exclusive processes: in one case (Vedretta d'Amola site), the intensity of facilitation stays the same while competition becomes more prevalent, whereas in another case (Western Trobio site), the intensity of competition remains constant and it is the decrease in facilitation that changes through time. Taken together, these results provide evidence in support of both the facilitation model of community structure (Connell and Slatyer, 1977; Van Andel et al., 1993; Bertness and Callaway, 1994) as well as the geoecological model of succession (Matthews, 1992; Whittaker, 1993). Since the deglaciation gradient is also a gradient of decreasing harsher conditions, these results also support the shift from competition to facilitation in harsher environments (Callaway and Lawrence, 1997).

Notably, pioneer and early species have significantly different trait values as compared to late species, having higher SLA and lower LDMC. The remarkable differences we found locate the least and most persisting species at two opposite sides of the “leaf economic spectrum” of plant traits (Wright et al., 2004; Díaz et al., 2016). One one hand, pioneer and early species, adapted to physical disturbance (“ruderal” sensu Grime, 2001), show more “acquisitive” strategies (high SLA) that facilitate their colonization and development in recent ice-free areas (Caccianiga et al., 2006; Gobbi et al., 2010). On the other hand, late species pose more “conservative” strategies (low SLA and high LDMC, “stress tolerant” sensu Grime, 2001) that can allow them dominating over other plants (Caccianiga et al., 2006; Gobbi et al., 2010) and persisting after glacier extinction. Everything considered, it emerges that heavier “conservative” leaves may allow late species excluding pioneer and early plants and dominating once glaciers disappear.

Our findings have threefold implications. First, they suggest that both competition and facilitation are acting together within the same ecosystem. Therefore, studying mechanisms of species coexistence and succession should not be limited to competition only (Callaway and Lawrence, 1997; Losapio et al., 2019), especially in mountain ecosystems. The modeling, data-driven approach we adopted is useful for establishing correlative patterns between species. These inferred associations are putative interactions, thus are helpful for generating testable hypotheses concerning biotic patterns and for developing a deepen understanding of the role of species interactions and geophysical forces in dynamic environments. Addressing interactions with other organisms, such as pollinators, herbivores, seed predators, and soil microbes will also broaden our knowledge of community-structuring processes and ecological succession.

Second, global warming is changing environmental conditions, shifting the relative importance of facilitation and competition, and shuffling species and communities over spacetime. In the long run, the consequences of habitat loss that follows glacier extinction are negatively impacting biodiversity both directly and indirectly. Directly, due to the lack of available habitats: glacier extinction will decrease the number and size of areas in proglacial environments, which will progressively disappear without glaciers. Since the species—area relationship is non-linear, a small reduction in glaciers may cause a much bigger species loss, reducing the occurrence of pioneer species. Indirectly, global warming exacerbates the impact of habitat loss via prompting competitive interactions and favoring fewer, dominant plants. Facilitation among pioneer plants may support species persistence in recently ice-free terrains. Indeed, it is well-established that plants can have beneficial, one-way or two-way, effects on the microenvironment and other plants, supporting community development (Connell and Slatyer, 1977; Bertness and Callaway, 1994; Chapin et al., 1994; Walker and del Moral, 2003). Nevertheless, glacier retreat is creating environmental conditions that disfavor “losers” pioneer species via creating an environment that favors competition and supports the “winners” late colonizers. Notice that the increase in competition with global warming would contribute to the local decline of late species too, beside the local extinction of pioneer plants. As experimental evidence supports that the fitness of high-altitude species is reduced due to competition with lower-altitude plants (Alexander et al., 2015), the same process would operate for the negative effects of competition between pioneer and late species.

Third, we highlight that there is no single factor, model or mechanism accounting neither for the response of species to global warming nor for the development of communities. Not only the effects of glacier retreat are uneven between species, which respond in different ways to global warming in general (Körner, 2005) and glacier retreat in particular (Raffl et al., 2006; Erschbamer, 2007), but also multiple factors interact with each other and a significant amount of uncertainty is associated with the future of glaciers and the development of communities. This includes the potential future (and past) actions of Homo sapiens, which have been overlooked in our estimates here. Interesting enough, almost half of the species pool shows inconsistent trends or weak responses to glacier retreat and environmental conditions, but random variation along the glacier foreland. Random events of colonization and dispersal as well as stochasticity associated with births and deaths may play a remarkable yet underestimated role in supporting plant persistence and driving successional change in community composition (Hubbell, 2001).

As we adopted a conservative scenario of glacier extinction ranging from 180 to 10,000 years (5,000 years on average), our predictions of species trends may be faster since it is likely that some of the glaciers examined here may disappear within the next decades. Although we used only two leaf traits with coarse average resolution, remarkable differences emerged between communities, providing meaningful insights into potential processes underlying species distribution patterns. Moreover, our estimates of plant diversity dynamic aid in overcoming the challenges associated with geochronology and unpredictable events, are important for furthering our understanding of biodiversity change, and could help in directing limited resources to critical species and communities. Although our predictions have a relatively high margin of uncertainty, making it difficult to estimate with accuracy and precision when exactly an extinction may happen, there is confidence on the overall pattern that indicates local species loss and extinction cascades following glacier extinction.

In conclusion, increasing the value of community survey data may help better understanding, conserving, managing and restoring ecosystems with scarce (economic) resources. To do so, addressing which and how species will interact with each other while being impacted by global warming is of paramount importance. Here, we have shown that plant diversity increases shortly after deglaciation, but declines in the long run following glacier extinction. This is particularly the case since pioneers and early colonizers are the most vulnerable to global warming and would not persist after glacier extinction. Despite facilitation can support plant establishment in recently-deglaciated terrains, the increase in competition over spacetime exacerbates habitat loss. Forecasting the fate of biodiversity change requires integrating hierarchical processes accounting for variation in the geophysical environment, biotic interactions and random events as well as embracing the challenge of including evolutionary changes, socio-economic drivers and feedbacks into a robust predictive science. Such framework shall be useful for guiding monitoring schemes and sustainable ecosystem management actions.

Data Availability Statement

All raw original data used in this study are published on Stanford Digital Repository, publicly available at https://purl.stanford.edu/hp416gs2665.

Author Contributions

GL conceived the study. MG and MC designed and coordinated data collection. GL, BC, CM, DT, and MC collected data. GL organized the dataset, performed statistical analyses and wrote the first draft of the manuscript. All authors contributed to the discussion and interpretation of data, revised the manuscript, and approved the submitted version.

Funding

This study was financed by the Swiss National Science Foundation through a grant awarded to GL (P2ZHP3_187938). MG acknowledges the support of the Province Autonomous of Trento (Italy).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank the Parco Naturale Adamello Brenta, Parco delle Orobie Valtellinesi, and Parco Nazionale dello Stelvio for providing sampling permission and accommodation.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2020.616562/full#supplementary-material

References

Alexander, J. M., Diez, J. M., and Levine, J. M. (2015). Novel competitors shape species' responses to climate change. Nature 525, 515–520. doi: 10.1038/nature14952