Vegetation Assembly, Adaptive Strategies and Positive Interactions During Primary Succession in the Forefield of the Last Venezuelan Glacier

Llambí, Luis D.; Melfo, Alejandra; Gámez, Luis E.; Pelayo, Roxibell C.; Cárdenas, Mariana; Rojas, Cherry; Torres, Jesús E.; Ramírez, Nerio; Huber, Bárbara; Hernández, Jesús

doi:10.3389/fevo.2021.657755

ORIGINAL RESEARCH article

Front. Ecol. Evol., 27 October 2021

Sec. Biogeography and Macroecology

Volume 9 - 2021 | https://doi.org/10.3389/fevo.2021.657755

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

Vegetation Assembly, Adaptive Strategies and Positive Interactions During Primary Succession in the Forefield of the Last Venezuelan Glacier

$\r\nLuis D. Llambí*$ Luis D. Llambí^1*

Nerio Ramírez⁵

Bárbara Huber^1,6

Jesús Hernández⁷

¹Instituto de Ciencias Ambientales y Ecológicas, Universidad de Los Andes, Mérida, Venezuela
²Centro de Física Fundamental, Departamento de Física, Universidad de Los Andes, Mérida, Venezuela
³Laboratorio de Dendrología, Facultad de Ciencias Forestales y Ambientales, Universidad de Los Andes, Mérida, Venezuela
⁴Instituto de Zoología y Ecología Tropical, Facultad de Ciencias, Universidad Central de Venezuela, Caracas, Venezuela
⁵Instituto de Protección Civil y Administración de Desastres, Mérida, Venezuela
⁶Laboratorio de Microbiología Molecular y Biotecnología, Facultad de Ciencias, Universidad de los Andes, Mérida, Venezuela
⁷Instituto Experimental Jardín Botánico “Dr. Tobías Lasser,” Universidad Central de Venezuela, Caracas, Venezuela

Glaciers are receding at unprecedented rates in the alpine tropics, opening-up new areas for ecosystem assembly. However, little is known about the patterns/mechanisms of primary succession during the last stages of glacier retreat in tropical mountains. Our aim was to analyze soil development and vegetation assembly during primary succession, and the role of changing adaptive strategies and facilitative interactions on these processes at the forefront of the last Venezuelan glacier (Humboldt Peak, 4,940 m asl). We established a chronosequence of four sites where the glacier retreated between 1910 and 2009. We compared soil organic matter (SOM), nutrients and temperatures inside vs. outside biological soil crusts (BSCs) at each site, estimated the cover of lichen, bryophyte and vascular plant species present, and analyzed changes in their growth-form abundance and species/functional turnover. We also evaluated local spatial associations between lichens/bryophytes and the dominant ruderal vascular plant (the grass Poa petrosa). We found a progressive increase in SOM during the first century of succession, while BSCs only had a positive buffering effect on superficial soil temperatures. Early seral stages were dominated by lichens and bryophytes, while vascular plant cover remained low during the first six decades, and was almost exclusively represented by wind dispersed/pollinated grasses. There was a general increase in species richness along the chronosequence, but it declined in late succession for lichens. Lichen and bryophyte communities exhibited a higher species turnover than vascular plants, resulting in the loss of some pioneer specialists as succession progressed. Lichen and bryophyte species were positively associated with safe-sites for the colonization of the dominant ruderal grass, suggesting a possible facilitation effect. Our results indicate that lichens and bryophytes play a key role as pioneers in these high tropical alpine environments. The limited initial colonization of vascular plants and the progressive accumulation of species and growth-forms (i.e., direct succession) could be linked to a combination of severe environmental filtering during early seral stages and limitations for zoochoric seed dispersal and entomophilic/ornithophilic pollination. This could potentially result in a slow successional response of these ecosystems to accelerated glacier loss and climate change.

Introduction

The process of plant community structuring during succession can be interpreted as being the result of several interacting filters which determine: (a) the species that can reach a site (biogeographic/dispersal filters); (b) the species that can tolerate conditions at the site (adaptive/functional filters); and (c) the species that are able to coexist at the site as a result of both positive and negative biotic interactions (Lortie et al., 2004). Some generalizations have been proposed on the changing importance of these processes for vegetation dynamics, including the key role played by effective seed dispersal, local micro-site availability and autogenic habitat amelioration via facilitation during the first stages of primary succession (Walker and Chapin, 1987; del Moral and Wood, 1993; Jones and del Moral, 2009). However, an integrated understanding of how these and other drivers operate across diverse taxa, succession types and biomes is still being developed (Walker and del Moral, 2003; Cauvy-Fraunié and Dangles, 2019; Pratch and Walker, 2020).

A significant contribution to the study of primary succession under climate change has come from research in glacier forefields, particularly following the worldwide tendency for glaciers to retreat after the Little Ice Age (LIA, Matthews, 1992; Fastie, 1995; Cannone et al., 2008). Glacier retreat rates have accelerated in the last five decades (Zemp et al., 2019), while there is increasing evidence that alpine species distributions are shifting upward (Lenoir and Svenning, 2015). This offers a unique opportunity to study climate change impacts in high mountain environments and opens up important questions regarding the development of novel ecosystems and the ability of high elevation specialists to become established and to maintain viable populations in newly exposed areas (Cannone et al., 2008; Zimmer et al., 2018; Cuesta et al., 2019; Anthelme et al., 2021).

Moreover, we are starting to witness the complete extinction of glaciers in many mountain regions, including Mediterranean and tropical alpine ecosystems, where most glaciers are likely to disappear during this century (Rabatel et al., 2013; Huss et al., 2017; Zemp et al., 2019). Alpine areas where glaciers will soon go extinct constitute particularly limiting environments because of steep slopes and dispersal barriers linked with island-like conditions near mountain peaks, high incident radiation but low partial pressures of CO₂ and O₂, incipient soil development and reduced water supply from glacier melt water (Zimmer et al., 2018; Cauvy-Fraunié and Dangles, 2019). In tropical alpine regions additional limiting factors include reduced plant protection from seasonal snow and particularly fast rates of high-elevation warming and glacier retreat (Rabatel et al., 2013; Wang et al., 2016; Vuille et al., 2018). However, studies of the patterns and mechanisms of primary succession in glacier forefields are scarce in the alpine tropics (see Spence, 1989; Mizuno, 2005; Mizuno and Fujita, 2014; Suárez et al., 2015; Zimmer et al., 2018; Anthelme et al., 2021; Rosero et al., 2021), especially under the harsh conditions that occur in mountains where they are about to become extinct.

The tropical Andes concentrate more than 95% of all tropical glaciers (Veettil et al., 2017). In the few available chronosequence studies of glacier forefields in the region, analyses of spatial association patterns and plant dispersal strategies have suggested that vegetation development could be slowed down by strong limitations in zoochoric seed dispersal and the effectiveness of facilitation mediated by nurse-plants (Zimmer et al., 2018; Anthelme et al., 2021). These studies have focused on the first two to six decades of succession, but how these findings apply for longer successional periods and to areas where remaining glaciers are less extensive or about to go completely extinct, remains largely unexplored. In fact, several essential aspects for understanding primary succession remain little understood in the alpine tropics, including the effect of biological soil crusts (BSCs) on soil development. Moreover, the available studies have focused on vascular plants, while community composition/dynamics and the role as pioneers of lichens and bryophytes have received limited attention [but see Sklenář et al. (2010) for a volcanic sere], despite being key structuring components of Arctic, Antarctic and alpine ecosystems, particularly during early succession (Longton, 1992; Bardgett and Walker, 2004; Cutler et al., 2008; Favero-Longo et al., 2012; Rosero et al., 2021).

In this context, a key unexplored aspect to understand vegetation assembly is to compare successional patterns of change in α and β diversity between lichens, bryophytes and vascular plants. Previous studies in glacial forefields have generally found an increase in vascular plant species richness (α diversity), at least for sites deglaciated after the LIA (Zimmer et al., 2018; Pratch and Walker, 2020). This has been associated with a progressive accumulation of plant species. Hence, in terms of the two main components of β diversity (nestedness and species turnover, see Baselga, 2010), during early primary succession vascular plants tend to show a predominance of nestedness (i.e., species accumulation) over temporal turnover (i.e., species replacement), a pattern called direct succession (Svoboda and Henry, 1987; Walker and del Moral, 2003; Pratch and Walker, 2020). However, incorporating bryophytes and lichens in these analyses could result in different patterns, as they could exhibit a different successional response to vascular plants. A higher species turnover could be expected in this case, as competitive displacement between lichens, bryophytes and vascular plants could become increasingly more important as vegetation develops. This could also result in a decrease of species richness for lichens/bryophytes in more advanced seral stages.

Another key aspect relates to a comparative analysis of changes in functional / growth-form strategies along succession. Previous research in the alpine tropics suggests that ruderal vascular plants tend to be composed of small grasses and non-graminoid herbs with a predominance of anemochoric seed dispersal strategies (Sklenář et al., 2010; Suárez et al., 2015; Zimmer et al., 2018; Anthelme et al., 2021), while a more diverse array of growth forms should be expected later in succession (Sarmiento et al., 2003). The hardy stress tolerant strategies of tropical alpine grasses from high elevation would also support the idea of them being the predominant growth-form under the limiting conditions expected during early succession (Márquez et al., 2006; Rada et al., 2019). In the case of bryophytes we would expect cushion and turfs to be the pioneer growth-forms during early seral stages, as they are well adapted to sun-exposed rocky surfaces and desiccation at high elevations (Vanderpoorten and Goffinet, 2009; Kürschner and Frey, 2012). Later in succession, we would expect the colonization/dominance of other growth forms such as bryophyte mats, better adapted to less exposed sites with more developed soils and higher humidity (Vanderpoorten and Goffinet, 2009). In the case of lichens, crustose growth forms are expected during early seral stages, as their bi-dimensional structure (closer to the substrate) has been found to be more resistant to desiccation produced by high radiation in rocky areas (Büdel and Scheidegger, 2008). In later stages we would expect the colonization/dominance of more complex growth forms like foliose, fruticose and dimorphic lichens (Sklenář et al., 2010; Bässler et al., 2015). Consequently, in terms of functional structure, strong environmental filtering during early succession could result in a progressive accumulation of growth-forms of lichens, bryophytes and vascular plants (i.e., a predominance of nestedness over growth-form turnover). Obviously, competitive displacement as succession progresses could also play a role, but it is more likely to modify the relative abundance of the different growth-forms than to result in the complete displacement of some functional groups.

Our general aim in this study was to analyze successional soil development and the assembly of lichen, bryophyte and vascular plant communities during the last stages of glacier retreat in the high tropical Andes, and to link them to potential underlying processes reflected in their changing functional strategies (i.e., growth-forms, dispersal and pollination syndromes) and local spatial interactions (i.e., facilitation/competition). A detailed historical reconstruction allowed establishing a chronosequence of four sites in the forefront of the last Venezuelan glacier at Humboldt’s Peak (4,945 m asl), in areas where the glacier retreated between 1910 and 2009 (c. 100 years, Ramírez et al., 2020). Our specific objectives were: (1) to analyze the role of BSCs on soil development, documenting successional changes in soil organic matter, nutrients and superficial soil temperatures in areas inside vs. outside their influence along succession; (2) to quantify changes in the composition, richness and abundance of lichen, bryophyte and vascular plant species and compare the degree to which species and growth-form accumulation vs. turn-over (the two components of between-site β-diversity) characterize successional community assembly; (3) to analyze the importance of abiotic vs. animal mediated dispersal and pollination, quantifying changes in the relative abundance of vascular plants with different dispersal/pollination syndromes; and (4) to explore the role of facilitation mediated by pioneer lichens/bryophytes on early vascular plant colonization, analyzing local spatial association between them and favorable sites for the establishment of the dominant ruderal plant (the grass Poa petrosa).

Based on the previous research on primary succession in temperate and tropical alpine regions, and considering the expected limitations for ecosystem development in glacier forefields of vanishing glaciers, we evaluated the following predictions: (1) the presence of BSCs will be associated with an increase in soil organic matter and nutrients; (2) a progressive accumulation of species (nestedness), more than species turnover, should characterize the dynamics of vascular plants under these extreme conditions, resulting in a successional increase in the richness of species and growth forms (i.e., direct succession). However, these patterns should differ if we incorporate lichens and bryophyte species in the analysis, reflecting a possible increase in species turn-over and competitive exclusion between them and with vascular plants as vegetation develops; (3) limitations in zoochoric seed dispersal and pollination by birds and insects in these extreme environments should result in less diverse vascular plant pollination and dispersal syndromes during early succession; and (4) the presence of lichens and bryophytes will be positively associated with favorable microhabitats for the early establishment of ruderal plants suggesting a facilitation role.

Materials and Methods

Study Area

The study was carried out in the North-West flank of the glacier forefront at Humboldt peak (8°33′0.8″ N, 71°59′50.56″ W), in the Sierra Nevada National Park, Venezuela (Figure 1). Venezuela will soon be the first Andean country to lose all of its glaciers. Ramírez et al. (2020) showed that glaciers in the Sierra Nevada de Mérida, where the highest peaks in the country are located, have lost 99% of their area since 1910. There is only one remaining ice mass at Humboldt Peak (4,942 m asl), which has lost 98.7 % of its surface since 1910, with an area for 2019 of only 0.0454 Km² (Ramírez et al., 2020). These glaciers have been retreating since the local expression of the LIA, with moraines at the terminus reaching c. 4,200 m elevation (Mahaney et al., 2000; Ramírez et al., 2020).

FIGURE 1

Figure 1. Study region maps. (A) Study region in the northern South America; (B) Detail of Mérida State, the Sierra Nevada National Park and Humboldt Peak, Venezuela. (C) General view of the North West face of Humboldt and Bonpland Peaks indicating the location of our study transects (red lines) and glacier cover for different moments between 1910 and 2020 [based on an integrated historical reconstruction by Ramírez et al. (2020)]. Glacier cover polygons drawn over a Google Earth image.

Average annual precipitation at the highest elevations in the Sierra Nevada de Mérida range from 1,000 to 1,200 mm, with high inter-annual variability (Andressen, 2007). Both precipitation and temperature decrease with elevation, from an average of 1,811 mm annual precipitation and 7.1°C mean air temperature at La Aguada station of the cable car system (3,446 m asl) to 1,173 mm and −0.4°C at Pico Espejo station (4,766 m asl, the closest station to Humboldt peak). Although there is little monthly variability in average temperature (≤3°C), precipitation is seasonal, with a dry period between December and March (see Andressen, 2007). Regarding the area’s geology, the Iglesias group formation is present at the highest elevations, corresponding to Precambrian high-grade metamorphic rocks including banded gneiss, schist, and amphibolite, with frequent granitic dikes and quartz veins (Schubert, 1975).

Sampling Sites

Field work was organized in three 1-week expeditions that occurred between May 2019 and December 2019. For site selection we used a chronosequence approach, based on a detailed reconstruction of glacial retreat in Humboldt Peak (1:5,000 scale, Ramírez et al., 2020). The reconstruction integrated historical photographs with the analysis of maps from 1910 (Jahn, 1925), high-resolution aerial photographs for 1952 and 1998 and satellite images for 2009 (Spot V, 2.5 m). This allowed the selection of four study sites along a successional/elevation gradient in the NW face of Humboldt’s peak (Figure 1): (a) Site 1: 10 years since glacial retreat, corresponding to the glacier position in 2009 (4,670 m asl); (b) Site 2: 21 years of succession, corresponding to 1998 (4,637 m asl); (c) Site 3: 67 years, for the 1952 line (4,537 m asl); and d) Site 4: 109 years, for the 1910 line (4,365 m asl). Although site replication would have been desirable to allow more robust generalizations, selecting more study sites on other mountain faces, logistic limitations prevented it: this is a remote area, requiring two full days to reach base camp; we worked on the more accessible face of the peak (Northwest), as working on other faces can be considerably more dangerous because of the steepness of the terrain.

The glacier forefront area corresponds to a periglacial landscape dominated by steep slopes (45–70°) formed by hard bedrock polished by glacier action. In each study site, we established one 50 × 20 m study area (1,000 m²) on sections with similar geomorphological and topographic characteristics, to allow a better comparison of vegetation and soil development between sites (Figure 1). All vegetation and soil sampling was carried out within each of these 1,000 m² areas (see details below and Supplementary Figure 1). Although it would have been desirable to nest the soil sampling points within the sampling plots used for studying vegetation cover (eight replicate plots, 5 × 5 m each), time/logistic constraints prevented this.

The 50 × 20 m study areas were placed parallel to the glacier front on slopes ranging between 10° and 30°, coinciding with rocky steps created by glacier erosion / scouring along rock foliation lines. The flatter local topography in the study sites allowed the accumulation of loose glacial till mixed with gelifraction materials (small angular rocks and scree), coming from both lateral slopes of the valley. We decided to work on these flatter rocky-steps as they allowed us to analyze areas less limiting for soil and vegetation development than the steep polished glacier bedrock walls, where plant colonization is almost completely absent.

Soil Analyses

Soil samples (0–3 cm depth) were collected in 12 random points along each study area, six in open soils and six directly under BSCs (see Supplementary Figure 1). Soils were transported to the laboratory, where they were dried and sieved using 4 and 0.5 mm mesh strainers to then determine pH (in water), total soil organic matter (SOM, ignition), total nitrogen (N_t, Kjeldahl), available phosphorus (P, Bray-Kurtz for acid soils), and exchangeable bases (calcium, magnesium, potassium and sodium, see details in Gilabert de Brito et al., 1990).

Soil properties were analyzed using a two-way Permutational Analysis of Variance (PERMANOVA+ for Primer 6.0 software, Anderson et al., 2008), defining the successional time as a fixed factor with four levels (10, 21, 67, and 109 years) and the presence/absence of BSCs as a fixed factor with two levels (inside vs. outside of BSCs) and each of the edaphic response variables taken separately (i.e., univariate tests). We used PERMANOVA instead of fully parametric ANOVA procedures, as our data was not normally distributed. In all PERMANOVA tests, the number of permutations was set to 9,999 and the probability of type I error was set at 95% (a = 0.05), although it was adjusted to 99% for post hoc tests (permutational t-tests).

At each study site, we also placed two thermometer data-loggers (Onset HOBO TidbiT v2), 2 cm below the soil surface, one inside and one outside of a BSC (always less than 2 m apart, Supplementary Figure 1). These recorded temperatures every hour for 6 months (between July 1st and December 10th 2019), allowing us to calculate for each situation average, minimum, and maximum superficial soil temperatures, average temperature amplitudes and the proportion of days with freezing superficial soil temperatures (see Supplementary Table 1). However, because of lack of replication (not enough temperature sensors for this), we were not able to statistically compare these results between sites/seral stages.

Patterns of Community Assembly and Diversity

To assess changes in community composition and species abundances, we first carried out a detailed inventory of total species richness in each site recording all the lichenized fungi, bryophyte and vascular plant species found within each 1,000 m² study area. Then, to quantify the abundance of each species, we visually estimated the cover of vascular plants in eight 5 × 5 m plots randomly placed within each study area (Supplementary Figure 1). We estimated lichen and bryophyte cover in three 0.5 × 0.5 m subplots placed along on the diagonals of each 5 × 5 m plot (and averaged results from the three sub-plots for subsequent analyses). Within each 5 × 5 m plot we also estimated the cover of rocks (≥ 50 cm in diameter), scree, loose plant necromass and bare soil/BSCs (with these superficial substrates summing 100% cover in each plot, without taking into account plant cover.

The visual estimation of cover for vascular plants, lichens and bryophytes was each done by a single person in all studied plots (i.e., the corresponding expert in the group). In the case of bryophytes and lichens, cover estimation was aided by placing a frame subdivided into one hundred 5 × 5 cm cells over the 50 × 50 cm sub-plot (with both experts working together to reach consensus). Confirmation consensus between the three experts was also reached for the cover of vascular plants in each plot. Admittedly, the fact that a different researcher estimated the cover of lichens, mosses and vascular plants could have introduced some unintended bias, but the observed differences in cover between groups are sufficiently large to make this small sampling bias unimportant.

We compared community composition and species abundance between 5 × 5 m plots with a Principal Coordinate Analysis (PCoA, Legendre and Legendre, 2003) based on a Bray Curtis similarity matrix (using Primer v6, Clarke and Gorley, 2006). Because of the high incidence of rare species, the original data was transformed with fourth root and then standardized by species for comparison.

Finally, we analyzed the degree to which the successional process involved species replacement (turnover) or species accumulation (nestedness) between sites, using a partitioned β-diversity analysis, following Baselga (2010). First, we separately established presence/absence matrices for lichens, bryophytes and vascular plants in each site using the full species inventory. We then established a similar matrix for all species of photosynthetic organisms taken together. Then, for each type of organism taken separately and for the whole group of species, we calculated total multisite β-diversity using Sørensen’s index (β_SOR), and decomposed it into two additive components: β-diversity resulting from turnover (using Simpson’s index, β_SIM) and β-diversity resulting from nestedness (β_NES). The same analysis was then repeated using a presence-absence matrix of all the growth-forms of lichens, bryophytes and vascular plants for each of the four study sites (see below for details on the grouping of species into growth-forms). These indices were calculated with the “beta.multi” function of the betapart 1.5.1 package (Baselga and Orme, 2012), using RStudio v. 1.1.463 (RStudio Team, 2016).

Changes in Functional Structure: Growth-Forms, Dispersal and Pollination Strategies

To assess successional changes in key aspects of the functional structure of lichen, bryophyte and vascular plant communities, we characterized the growth-form of each species. For lichens, we used the classical system including three growth forms (crustose, foliose and fruticose, see Büdel and Scheidegger, 2008). We defined dimorphic lichens as a mixed form with a primary thallus usually crustose or squamulose/foliose, and a secondary fruticose thallus. For bryophytes, we used Vanderpoorten and Goffinet’s (2009) system, based on the direction of main shoots. The categories included “cushion” (radiating from a central point); “turfs” (erect habit of acrocarps); and “mats” (creeping habit of pleurocarps). Finally, for vascular plants, we adapted the framework proposed by Hedberg and Hedberg (1979), adding to the five growth forms defined by them (grasses, cushions, shrubs, acaulescent rosettes and caulescent rosettes), two important growth forms in many alpine ecosystems: non-graminoid herbs and ferns.

For vascular plants only, we also defined the pollination and dispersal syndrome of each species present, based on on-site observations of their floral visitors (75 h of onsite field observations) and the morphology of their flowers and diaspores, following classifications by Faegri and van der Pijl (1979) and Tovar et al. (2020), respectively. For pollination syndromes we also relied on previous detailed studies of observed plant-pollinator interactions in the high Venezuelan páramo (Pelayo et al., 2019, 2021). Hence, plants were classified as anemophilous, entomophilous or ornithophilous for pollination and as anemochoric, ballochoric, barochoric or zoochoric in terms of dispersal. Given their importance during early colonization, we sub-divided anemochoric species in two groups, those with minute/dust seeds (<3 mm, in our Poaceae species), and those with additional morphological adaptations for wind dispersal (i.e., Asteraceae with plumose, hairy or winged seeds). The total cover of each growth form or syndrome was calculated as the sum of the cover of each species included in a given category in each sampling plot. These were then averaged in each study site to calculate the relative cover of each growth form or dispersal/pollination syndrome in each site.

Spatial Association Between Lichens/Bryophytes and Ruderal Vascular Plants

The possible role of lichens and bryophytes on the establishment of the dominant early- succession vascular plant species (the ruderal grass Poa petrosa) was studied in Sites 2 and 3 (21 and 67 years), using their local spatial association as a proxy. We did not include Site 1 (10 years) because the density of P. petrosa was too low, while in the case of Site 4 (109 years) a much higher vegetation cover of many different plant species complicated establishing clear local spatial association patterns. Moreover, this ruderal grass was no longer among the most abundant vascular plant species at this site.

In sites 2 and 3, we compared the cover of all bryophyte and lichenized fungi species (visual estimation, always by the same two experts working together), within 10 circular micro-plots (30 cm diameter) centered on mature P. petrosa plants (≥5 cm height) and 10 paired micro-plots placed outside, centered on a randomly chosen point within 1–2 m from the target plant (for a total of 40 microplots, 20 in each study site). The small size of the micro-plots facilitated visual estimation.

We compared the total cover of lichens and bryophytes (taking them as a group and also for each species separately) using the average of the relative interaction index [RII = (C_n-C_a)/(C_nC_a)], where C_n is the cover of the organism type or species near P. petrosa, and C_a the cover of the same organism or species away from P. petrosa (Armas et al., 2014). RII values range from 1 to −1, where significant positive values (i.e., 95% CI does not include zero) would indicate a positive local association of lichens/bryophytes with favorable sites for the growth of the pioneer grass, and significant negative values would suggest negative, inhibitory effects.

Results

Soil Development and Substrate Cover Changes

Soil organic matter (SOM), total soil nitrogen (N_t), calcium (Ca), magnesium (Mg) and available phosphorus (P) showed statistically significant differences between sites along succession (p ≤ 0.0001), while the rest of the variables did not (Table 1). In the case of SOM and N_t there was a clear trend to increase from the youngest to the oldest sites along the chronosequence. However, Ca and Mg showed the opposite trend, with significantly higher concentrations in the two youngest sites in the case of Mg, and in Site 1 (10 years) for Ca. P was significantly higher in Site 2 (21 years) than in the rest of the sites. The presence of BSCs did not show a significant effect for any variable, and there was no significant interaction between time and BSC presence (p ≥ 0.05 for the presence of BSC and the interaction term in all the two-way PERMANOVAs for each variable). However, BSCs consistently reduced superficial soil temperature amplitudes in all study sites. This was mainly the result of a consistent reduction in maximum temperature, but the proportion of days during the study period with freezing temperatures was also consistently lower within BSCs than outside (Supplementary Table 1).

TABLE 1

Table 1. Average (standard error) for soil properties (0–3 cm) studied within biocrusts (Bio) and in paired areas outside (Out) in four sites (S1 to S4) along a glacier chronosequence (years since glacier retreat) in the forefront of the Humboldt glacier, Venezuela.

Average soil/BSCs cover was lowest (1.8 % ± 1.3 SE) in Site 1 (10 years), while it ranged between 10.3% ( ± 3.2) and 8.7 % ( ± 1.2) in Sites 2–4. Rock cover varied between 60.6% ( ± 6.0) and 49.4% ( ± 6.5) in the three earlier successional sites and was lowest in Site 4 (28.1% ± 7.4). Average scree cover ranged between 50.0% ( ± 6.4) in Site 1 and 28.4 % ( ± 4.7) in Site 2, with intermediate values in the other two sites. Finally, plant necromass cover ranged between 0.0 and 1.0% in the three younger sites and was highest in Site 4 (12.1 % ± 2.2).

Patterns of Community Assembly and Diversity

Based on our complete inventory in each site (1,000 m² sampling areas), a total of 47 species of lichenized fungi were identified across the four seral stages (19 genera and 16 families). Of these, 25 species (53%) were new reports for Venezuela (Supplementary Table 2). For bryophytes, we identified a total of 55 species (27 genera and 16 families). Six were new records for Venezuela and 8 endemic to the Andean páramos (15%). Finally, 36 species of vascular plants were registered inside the four transects/sites (13 families and 30 genera). There were no new records for Venezuela, but 24 species were endemic to the Andean páramos (67%) and 8 endemic to Venezuela.

Total species richness increased monotonically along succession for bryophytes and vascular plants, but not for lichens, which showed a maximum in site 3 (67 years, Figure 2). Average richness was always lowest in Site 1 (10 years), and increased in the older sites, but leveled-off in Site 3 for lichens and bryophytes, while it continued to increase for vascular plants, with a maximum in Site 4 (109 years). Interestingly, bryophytes showed higher total and average richness in the two earlier sites than plants and lichens, but lichens showed much higher average cover in these initial stages. Average lichen cover peaked after 21 years (29.2%), and then declined, while bryophytes and vascular plants monotonically increased their cover along succession, from very low values in Site 1 (1.17 and 0.38% respectively) to a maximum in Site 4 (Figure 2).

FIGURE 2

Figure 2. Changes in total species richness and average richness (A) and cover (B) for lichens, bryophytes and vascular plants in four study sites (S1-S4) along a post-glacial chronosequence (10–109 years of succession) at Humboldt Peak, Venezuela. Different letters indicate significant differences across successional stages for each group of organism analyzed separately (i.e., differences between seral stages, permutational t-test, α = 0.01).

The PCoA ordination analysis indicated very clear successional changes in species composition and cover (Figure 3). Sampling plots were organized in terms of successional time along the 1st axis (31.7% of total variance), 10 year plots occupying the left side and 67 / 109 year plots the right-hand side. The two late-succession sites were clearly separated along the 2nd axis (14.7% variance). The two earliest seral stages were dominated by lichens (e.g., Rhizocarpon umbilicatum). Bryophytes showed much lower cover in Site 1 (10 years), with Molendoa peruviana and Racomitrium crispipilum being among the dominant pioneers. In the next three sites Racomitrium subsecundum and R. crispipilum became the dominant moss species (see Supplementary Table 2 for details). Only three vascular plant species were present in Site 1: the grass Poa petrosa, and two Asteraceae species. However, only P. petrosa was present within the 5 × 5 m sampling plots, showing a very low average cover (0.38%). Other species progressively colonized the next two sites, with grasses dominating Sites 2 and 3 (e.g., Festuca fragilis). Finally, in Site 4 (109 years) a more diverse community developed, dominated by the shrubs (e.g., P. imbricatifolia), the grass F. fragilis, erect herbs (e.g., Castilleja fissifolia) and the giant caulescent rosette Espeletia moritziana (see details in Supplementary Table 2).

FIGURE 3

Figure 3. Principal Coordinates Analysis based on a Bray–Curtis similarity matrix of the cover of all species present in 5 × 5 m plots (n = 32) in four study sites (S1-S4) along a post-glacial chronosequence at Humboldt Peak, Venezuela. The vectors correspond to the Spearman Rank correlation (r_s ≥ 0.5) of the dominant species of lichenized fungi (red), bryophytes (blue) and vascular plants (green) with the first two principal axes (full species names in Supplementary Table 2).

The analysis of multisite β-diversity revealed contrasting responses for plants, bryophytes and lichens. We found a highest total β-diversity and a highest contribution of nestedness in the case of vascular plants, indicating a predominance of species accumulation over turnover. The lowest total β-diversity and the highest contribution of temporal turnover was observed for bryophytes, while an intermediate situation was observed for lichens (see Table 2 and Supplementary Figure 2). Interestingly, when the three groups of organisms were analyzed together, the species turnover component showed a slight predominance over the nestdness component (54.7 vs. 45.3%), reflecting the higher relative importance of species turnover shown by bryophytes and lichens compared to vascular plants.

TABLE 2

Table 2. Analysis of multi-site β-diversity for lichens, bryophytes and vascular plants and for the whole vegetation analyzing all species of photosynthetic organisms together.

Changes in Functional Strategies

In the case of lichenized fungi all sites were dominated by crustose species, while fruticose lichens were the sub-dominant growth form (Figure 4). Foliose lichens were present from the earliest stages and increased in relative cover, but remained below 10% in all sites, while dimorphic species were only present in the two oldest sites, with a maximum of 4.5% of relative cover in Site 4 (109 years). For bryophytes, the three growth-forms were present across the four seral stages, but turfs were clearly dominant during early seral stages, while in late succession mats dominated. Small cushions showed a low relative importance in all sites/seral stages. Finally, regarding vascular plants, growth-form richness increased across sites, grasses strongly dominating the three youngest sites with values between 100 and 85% of relative cover. Vegetation showed a much higher growth form richness in Site 4, with shrubs becoming dominant (37.5% of relative cover), followed by grasses and non-graminoid herbs. Caulescent and acaulescent rosettes, cushions and ferns were less abundant, with relative cover values below 6% in all cases (Figure 4).

FIGURE 4

Figure 4. Average relative cover (%) of the different growth forms of lichens (A), bryophytes (B) and vascular plants (C) in four study sites (S1-S4) along a post-glacial chronosequence (10–109 years of succession) at Humboldt Peak, Venezuela.

When multisite β-diversity was analyzed in terms of the presence-absence of growth-forms of lichens, bryophytes and vascular plants taken together, its total value was lower than the value based on the presence/absence of individual species, and there was no observed growth-form turnover. Hence, the observed β-diversity between sites was entirely due to the progressive accumulation of growth forms along succession, with an observed increase in growth-form richness for lichens and vascular plants and no loss from early to late succession, so that the nestedness component comprised 100% of the overall β-diversity in this case.

Regarding plant dispersal and pollination strategies, anemochoric grasses with small seeds and anemophilous pollination dominated vegetation in the first three sites (Figures 5A,B). In contrast, plants with anemochoric seeds with appendages (e.g., plumes/hairs, Asteraceae) became dominant in Site 4 (109 years, with 55.2 % relative abundance), followed by dust seed grasses (27.8%), and a much smaller representation of ballochoric and barochoric species (Figure 5A). We found no species with zoochoric dispersal in any of the four study sites. Pollination syndromes also became more diverse in Site 4, with a dominance of insect pollinated species (entomophily), followed by anemophilous and ornithophilous species (Figure 5B).

FIGURE 5

Figure 5. Average relative cover (%) of vascular plant species with different dispersal (A) and pollination (B) syndromes in four study sites (S1-S4) along a post-glacial chronosequence (10–109 years of succession) at Humboldt Peak, Venezuela.

Spatial Association Between Lichens/Bryophytes and Ruderal Plants

The relative interaction index (RII) indicated that bryophytes were significantly overrepresented in the neighborhood of the ruderal grass Poa petrosa, suggesting a positive association in both sites. In contrast, total lichen cover showed a neutral spatial relation, with the RII not being significantly different from zero (Figure 6). When RII was analyzed at the species level (see Supplementary Figure 3), two bryophyte species showed a significant positive association with the grass in Site 2 (21 years), R. subsecundum and Bryum cf. chryseum (among the most abundant species in this site, being particularly frequent within BSCs). Only one species showed a significant negative RII (in Site 2), the saxicolous lichen Tremolecia atrata. In Site 3 (67 years) four bryophytes showed a significant positive RII, including the abundant moss C. flexuosus var. flexuosus (the 2nd most common species within BSCs). In this site four lichenized fungi species also showed a positive association with the grass; three of them were dominant within BSCs and are known to have atmospheric N-fixing photobionts (e.g., S. verruciferum). No bryophyte or lichen showed significant negative spatial associations with the pioneer grass in S3.

FIGURE 6

Figure 6. Mean relative interaction index (RII) and 95% confidence interval for the total cover of lichens (black circles) and bryophytes (white circles) in paired sampled in rings (n = 40) centered on mature individuals of the pioneer grass Poa petrosa vs. random sites outside, in two study sites (S2-S3) along a post-glacial chronosequence (21 and 67 years) at Humboldt Peak, Venezuela.

Discussion

Climate change is resulting in accelerated glacier retreat in the tropical Andes (Rabatel et al., 2013; Zemp et al., 2019), exposing new areas for ecosystem development at high elevations. However, little is known in the alpine tropics about the dynamics of primary succession under the extreme conditions that occur during the last stages of glacier retreat. Our results in the forefield of the last Venezuelan glacier show a strong dominance of lichens and then bryophytes during the first six decades of succession and a progressive accumulation of vascular plant species and growth forms characteristic of direct succession (Svoboda and Henry, 1987; Pratch and Walker, 2020). Moreover, they suggest that in the alpine tropics lichens and bryophytes could play an important role in modulating early vascular plant colonization, generating favorable micro-sites, particularly within BSCs. However, they also support the idea that animal mediated dispersal and pollination could be limited in these extreme conditions, with a complete absence of zoochoric (animal dispersed) vascular plant species during the first century of succession and a strong dominance of anemophily (wind pollinated) grasses in the first six decades (see also Zimmer et al., 2018). Below, we discuss in more detail the dynamics of soil and vegetation development, and some of the key processes that could help us understand the limitations faced by ecosystem development during the first 100 years of post-glacial succession in the highland tropics.

Soil Development

Against our expectations (prediction 1), BSCs did not show a significant effect on SOM or soil nutrients in any of our study sites. This does not agree with what has been documented in several other glacial seres, and the fact that N-fixing organisms are frequently associated with BSCs (Longton, 1992; Belnap and Lange, 2001; Walker and del Moral, 2003; Bardgett and Walker, 2004; Breen and Lévesque, 2006, 2008). Even so, BSCs appeared to reduce daily superficial soil temperature variability (see also Breen and Lévesque, 2008) and the number of days with below zero temperatures (72–21% lower frequency within BSCs in the three highest sites). This could be important in terms of microhabitat amelioration for colonizing plants, as it could reduce soil instability and soil water limitations associated with daily freeze-thaw cycles (Pérez, 1995; Ramírez et al., 2015).

However, in agreement with what has been frequently reported for glacial primary succession (Matthews, 1992; Walker and del Moral, 2003; Pratch and Walker, 2020; Rosero et al., 2021), we found a clear successional increase in superficial SOM and N_t. Values were relatively high considering the limiting conditions for soil development in the study area (e.g., steep slopes with hard exposed metamorphic bedrock and low vegetation cover), although they remained below typical values reported in the superpáramo belt in Venezuela between 4,000 and 4,400 m asl (Malagón, 1982; Pérez, 1995). The accumulation of soil N_t could be partly due to the presence of several species of lichenized fungi which are known to involve nitrogen fixing cyanobacteria as primary (e.g., Sticta, Peltigera, and Cora) or secondary photobiont (e.g., species of Stereocaulon and Placopsis; Rikkinen, 2002; Lücking et al., 2013).

Interestingly, Ca and Mg showed higher values in the two earliest sites and then sharply declined. These initial high content of Ca and Mg could be linked with the presence of amphibolite rocks rich in these cations in the study area (Schubert, 1975). The successional decline of Ca and Mg could be associated with their accumulation in microbial or plant biomass and/or their loss via leaching (given the typically low cation exchange capacity of these sandy soils with low SOM; Malagón, 1982), as reported in other post-glacial succession studies in temperate and tropical sites (Walker and del Moral, 2003; D’Amico et al., 2014; Whelan and Bach, 2017; Anthelme et al., 2021).

Patterns of Vegetation Assembly

Lichens strongly dominated the two early sites (mostly saxicolous species) and shared dominance with bryophytes in Site 3 (67 years, mostly terricolous species present within BSCs). Vascular plants only became the dominant structural component of vegetation after one century of succession (see Figures 2, 3). A similar pattern was reported by Sklenář et al. (2010) in volcanic sere in Cotopaxi, while Rosero et al. (2021) found a relative dominance of bryophytes over lichens and vascular plants during the first 65 years of post-glacial succession in a vanishing glacier in Carihuairazo (both in Ecuador). In agreement with our 2nd prediction, plant species richness increased during the first one hundred years of succession, which paralleled the observed increase in SOM and plant cover. This has been reported in most primary succession studies in glacier forefields, at least for sites deglaciated after the LIA (Zimmer et al., 2018; Pratch and Walker, 2020; but see Anthelme et al., 2021). However, in the case of lichens, total richness peaked after 67 years, but then declined (see also Sklenář et al., 2010).

The general increase in species richness can be linked with the documented predominance of species accumulation over turnover during the initial stages of post-glacial succession (Svoboda and Henry, 1987; Walker and del Moral, 2003; Pratch and Walker, 2020). This was the pattern observed for vascular plants: species accumulation with almost no species loss; of the 36 species recorded across the four sites, 34 species were present in Site 4, with only two high elevation Draba species missing (see Supplementary Figure 1). This was confirmed by decomposing multisite β-diversity in its nestedness and turnover components, with a limited between-site turnover between sites, corresponding to only 17.5% of total β-diversity.

Interestingly, the response was very different when lichens and bryophytes were considered (as proposed in our 2nd prediction). Multisite β-diversity was lower for lichens and bryophytes, but species turnover accounted for a much higher proportion of it, representing 47% of total β-diversity for lichens and 69% for bryophytes. This also resulted in a higher turnover component when all three groups of organisms were analyzed together than when vascular plants were analyzed separately. Both a higher turnover for bryophytes and lichens, and the decline of total species richness during late succession for lichens, suggest that conditions can become less favorable as succession proceeds for some specialist early-succession cryptogamic species. This could be the result of increased competition (among them, and also with vascular plants), and/or a reduced availability of suitable micro-habitats (e.g., open rocky sites for pioneer saxicolous species). Hence, as glaciers completely disappear and vegetation develops, some specialist species (e.g., mosses like Molendoa peruviana and Grimmia fuscolutea) could be lost in the long-term from these periglacial environments (see review by Cauvy-Fraunié and Dangles, 2019).

Changes in Functional Strategies

In the case of lichens and vascular plants, there was a clear increase in growth-form richness, which was maximum after 109 years (hypothesis 2), resulting from a progressive accumulation of growth-forms. This was clearly reflected by decomposing β-diversity in its nestedness and turnover components for the presence-absence of growth-forms, with nestedness accounting for 100% of total β-diversity in this case. For lichens, only two growth forms were present in Site 1 (10 years), while four different growth forms were found after 109 years. However, all sites were dominated by crustose species, with fruticose lichens being sub-dominant. In contrast, Sklenář et al. (2010) found fruticose lichens to dominate their youngest lahar sites, while foliose lichens were more abundant in intermediate stages, possibly as a result of more developed soils and a lower cover of exposed bedrock in their older volcanic study region. In the case of bryophytes, we observed no overall change in the number of growth forms present, but turfs showed a higher relative abundance during early succession. This growth-form is typical in disturbed, sun-exposed areas with high light intensity and low humidity (Vanderpoorten and Goffinet, 2009). In contrast, mats increased their relative cover in the late succession site, possibly as a result of lower sun exposure associated with much higher vascular plant cover in this site.

For vascular plants, Site 4 showed a much higher functional diversity than the three earlier sites. Grasses constituted by far the most important early-succession growth form, strongly dominating the first three sites, while a few non-graminoid herbs, small shrubs and acaulescent rosettes (all Asteraceae) were able to colonize Sites 1 and 2, but showed extremely low cover values. Grasses, non-graminoid herbs and small shrubs from these two plant families (Poaceae and Asteraceae) have also been found to be among the most important ruderal species in available post-glacial succession studies in the tropical Andes (Suárez et al., 2015; Zimmer et al., 2018; Anthelme et al., 2021; Rosero et al., 2021).

The dominance of grasses during early succession in our high elevation sites could be associated with them including many of the most tolerant species to low/freezing temperatures and seasonal drought in the alpine tropics (Márquez et al., 2006; Rada et al., 2019). In Site 4 vegetation showed higher functional diversity, with seven different growth forms and a more even distribution of relative abundance. Here shrubs, grasses and non-graminoid herbs were the most abundant, while the cover of giant rosettes and cushion plants remained low compared to typical superpáramo sites at these elevations in the Venezuelan Andes (Cáceres et al., 2015; Cuesta et al., 2017; Hupp et al., 2017; Mora et al., 2019).

The strong dominance of ruderal grasses in the first six decades of succession was associated with the prevalence of anemochoric dispersal via small dust seeds and anemophilous pollination syndromes. Plants with specialized structures for seed dispersal (hairs/plumes), like those of many Asteraceae, became dominant in Site 4 (109 years). In this site, plants with entomophilous pollination syndromes became dominant, while ornithophilous species were also present. No zoochoric species were found in any of our study sites, although they can have higher relative importance in alpine communities at lower elevations across the tropical Andes (see Tovar et al., 2020). These results support our 3rd prediction regarding the low prevalence seed dispersal and pollination mediated by insects and birds in early post-glacial succession, which has been reported both in tropical and temperate alpine seres (Walker and del Moral, 2003; Jones and del Moral, 2009; Erschbamer and Mayer, 2011; Zimmer et al., 2018; Anthelme et al., 2021; Rosero et al., 2021).

Spatial Association Between Lichens/Bryophytes and Ruderal Plants

Our plant-centered analysis of favorable micro-sites for the colonization of the dominant ruderal grass (Poa petrosa), indicated a significant spatial association with the local cover of bryophytes, but not with lichens when analyzed as a group. Even so, when their cover was analyzed at the species level, we did find several bryophyte and lichenized fungi species that were significantly associated with P. petrosa, supporting our 4th prediction. Most of these species were dominant within BSCs, and three of the lichens are known to have atmospheric N-fixing photobionts (Rikkinen, 2002; Lücking et al., 2013). Hence, our results support the idea that in extreme tropical alpine areas, facilitation mediated by BSCs, more than by nurse plants, could be an important mechanism for the establishment of ruderal plants (see also Zimmer et al., 2018). In fact, growth forms reported to act as effective nurse plants in superpáramo ecosystems were absent or rare in our three early successional sites (i.e., compact shrubs, cushions and giant rosettes, see Cáceres et al., 2015; Hupp et al., 2017; Mora et al., 2019). However, without further experimental evidence, we cannot fully rule out the alternative explanation of shared local habitat preferences between lichens/bryophytes and P. petrosa, although the paired sampling design used should have helped to reduce the likeliness of this alternative explanation.

Conclusion

The results presented here emphasize the key role played by lichens and bryophytes as pioneer photosynthetic organisms in these high tropical alpine ecosystems, and their potential for ameliorating local conditions within BSCs (buffering superficial soil temperature changes). Their local spatial association with the most abundant grass during early succession also supports the idea that they can facilitate the colonization of ruderal vascular plants (see review in Belnap and Lange, 2001). Moreover, while vascular plants exhibited a response characteristic of a direct succession, with a progressive accumulation of species and limited species relay (Walker and del Moral, 2003), there was a higher species turnover in the case of lichens and bryophytes, indicating that different mechanisms, including competitive exclusion, could drive their long-term successional response. This could imply the loss of some pioneer specialists from periglacial environments where a denser vegetation cover develops over time, such as saxicolous species characteristic of rocky open sites (Cauvy-Fraunié and Dangles, 2019).

More generally, our results indicate that under the extreme conditions associated with the last stages of glacier retreat in the northern Andes, ecosystem development can face important barriers, particularly during early seral stages (see also Zimmer et al., 2018; Anthelme et al., 2021). The limited initial colonization of vascular plants and the progressive accumulation of species and growth-forms during the first 100 years of succession (i.e., a nested pattern characteristic of direct succession) could be linked with a combination of: (a) strong dispersal limitations, with a complete absence of animal dispersed vascular plants; (b) severe environmental filtering, with ruderal vascular plants being strongly dominated by stress tolerant grasses; and (c) limitations in the effectiveness of positive biotic interactions, with a complete absence of facilitation mediated by nurse-plants and of insect/bird pollination during early succession. This could have important implications for understanding climate change impacts in high tropical alpine ecosystems, especially near mountain summits. In particular, these results suggests that the projected colonization and upward shifts in the distribution of high alpine communities and plant species linked with global warming (Lenoir and Svenning, 2015) could face important abiotic and biotic limitations in the forefields of vanishing glaciers, which need to be more explicitly incorporated into prospective/predictive models of vegetation dynamics (e.g., Tovar et al., 2013).

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

LL and AM conceived the project. LL, AM, MC, LG, RP, and CR designed the study. LL, AM, LG, MC, CR, RP, JT, and BH collected field data. NR, AM, and LL mapped the chronosequence. LG, MC, CR, and JH carried out all taxonomic identifications. LL and MC performed statistical analyses and wrote the first draft of the manuscript. All authors contributed to the discussion and interpretation of data, revised and complemented the manuscript and approved the submitted version.

Funding

This study was funded by the National Geographic Society (Grant No. NGS-55170R-19).

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.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We would like to thank the constructive comments of three referees and the editor which greatly improved the manuscript. We would also like to thank the Mérida Cable Car System (Mucumbarí) and the Instituto Nacional de Parques for their assistance in access to the study area in Humboldt Peak. We would also further like to thank URO and UGAM mountaineering organizations for their invaluable assistance and logistic support during field expeditions. This article is dedicated to the incredible spirit of our colleagues at Universidad de Los Andes, for their perseverance during challenging times for scientific research in Venezuela.

Supplementary Material

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

References

Anderson, M., Gorley, R., and Clarke, K. (2008). PERMANOVA+ for PRIMER: Guide to Software and Statistical Methods. Plymouth: PRIMER-E.

Google Scholar

Andressen, R. (2007). “Circulación atmosférica y tipos de climas,” in GeoVenezuela 2. Medio Físico y Recursos Naturales, ed. P. C. Grau (Caracas: Fundación Polar), 238–329.

Google Scholar

Anthelme, F., Cauvy-Fraunié, S., Francou, B., Cáceres, B., and Dangles, O. (2021). Living at the edge: increasing stress for plants 2–13 years after the retreat of a tropical Glacier. Front. Ecol. Evol. 9:584872. doi: 10.3389/fevo.2021.584872