# Llamas (Llama glama) enhance proglacial ecosystem development in Cordillera Blanca, Peru

Anaïs Zimmer, Timothy Beach, Sebastián Riva Regalado, Jean Salcedo

Aliaga, Rolando Cruz Encarnación, Fabien Anthelme

► To cite this version:

Anaïs Zimmer, Timothy Beach, Sebastián Riva Regalado, Jean Salcedo Aliaga, Rolando Cruz Encarnación, et al. Llamas (Llama glama) enhance proglacial ecosystem development in Cordillera Blanca, Peru. Scientific Reports, 2023, 13, 10.1038/s41598-023-41458-x. http://dx.doi.org/10.1038/s41598-023-41458-x

HAL Id: hal-04217236

https://hal.inrae.fr/hal-04217236

Submitted on 25 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for updates

OPEN

# Llamas (*Llama glama*) enhance proglacial ecosystem development in Cordillera Blanca, Peru

Anaïs Zimmer<fnref n="1" />✉, Timothy Beach<fnref n="1" />, Sebastián Riva Regalado², Jean Salcedo Aliaga³, Rolando Cruz Encarnación⁴,⁵ & Fabien Anthelme⁶

Worldwide, mountain glaciers are shrinking rapidly. Consequently, large areas are becoming available for the development of novel alpine ecosystems. These harsh environments, however, delay primary succession. In this study with a local community, we conducted an inclusion experiment to investigate whether *Llama glama* influences soils and vegetation primary succession following glacial retreat. At the foot of the Uruashraju glacier in the Cordillera Blanca, Peru (~4680 m.a.s.l.), we established four llama inclusion plots and four control plots that we studied from 2019 to 2022, 24–40 years after deglacierization. After three years, the llama plots had significantly increased soil organic carbon and soil nitrogen. In the llama plots, we found a large, significant increase in vascular plant cover (+57%) between the second and third years of experimentation, and we identified four new species that were not present in 2019. Our results suggest that *Llama glama*, through their latrine behavior and role as a seed disperser, enhances the primary succession and novel ecosystem formation in recently deglacierized landscapes. Our study provides scientific support that rewilding of native Andean camelids may favor adaptation to glacier retreat and inform conservation and management strategies in proglacial landscapes.

Glaciers are retreating rapidly, and 49 to 83% of glaciers worldwide (excluding ice sheets) will disappear by 2100¹. The opening glacier forelands provide an opportunity for novel alpine ecosystems to develop. Here, we investigate to what extent llamas may influence soils and vegetation primary succession following glacial retreat and, thus, influence novel proglacial ecosystem formation on a timescale relevant to the rapid pace of anthropogenic climate change.

Since the beginning of the twentieth century, much work has been done on primary succession in proglacial environments², including recent key advances on biotic colonization and pedogenesis³,⁴. The main ecological constraints to plant succession after glacier retreat include dispersal limitation², deficit in facilitators⁵, and abiotic stress from low temperature, low soil moisture, and low soil fertility⁶. In general, the first colonizers are wind-dispersed, while animal- or water-dispersed individuals are underrepresented⁷. Nevertheless, high dispersal capacity alone does not ensure successful succession. Indeed, the very first colonizers after glacial retreat may be opportunistic species that disappear when water and temperature change as the glacier recedes⁸. Thus, a dispersal lag and an establishment lag constrain primary succession⁹. Nurse species are mostly absent during early succession⁵, which limits facilitation when other species need these positive interactions most to overcome extreme environmental stress (stress-gradient hypothesis¹⁰). Establishing and organizing novel alpine plant communities in proglacial landscapes is therefore complex and slow, and studying primary succession after glacier retreat necessitates long-term observation and experimentation⁴. Manipulative experiments offer great potential for understanding the role of biotic and abiotic factors in driving seed dispersal and plant establishment during primary succession¹¹. Notably, in situ experiments might inform us about how external drivers (e.g., livestock grazing, wildlife re-introduction, and tourism) can transform primary succession trajectories.

Crucial to understanding primary vegetation succession is identifying the drivers of proglacial pedogenesis. First, proglacial pedogenesis is slow: skeletic or lithic Leptosols require decades to form and more distinct soil horizons may take hundreds of years¹². Proglacial soil formation rates vary depending on topography and

<footnote n="1">¹Department of Geography and the Environment, University of Texas at Austin, Austin, TX, USA. ²Laboratorio de Florística, Departamento de Dicotiledóneas, Museo de Historia Natural, Universidad Nacional Mayor de San Marcos, Lima, Peru. ³Departamento de Etnobotánica y Botánica Económica, Museo de Historia Natural. Universidad Nacional Mayor de San Marcos, Lima, Peru. ⁴Área de Evaluación de Glaciares y Lagunas, Autoridad Nacional del Agua, Huaraz, Peru. ⁵Universidad Nacional Santiago Antúnez de Mayolo, Huaraz, Peru. ⁶AMAP, Univ Montpellier, IRD, CIRAD, CNRS, INRAE, Montpellier, France.  email: anais.zimmer@utexas.edu</footnote>
morphodynamics¹³, climate factors¹⁴, parent materials¹⁵, time¹³, and organisms (i.e., plant-soil interactions¹⁶, microbial interactions¹⁷, and biological soil crusts¹⁸). Proglacial soils are usually phosphorus and nitrogen limited⁶. Although microorganism colonization is key during soil development, harsh proglacial microclimates and nutrient limitations impede microbial activity⁶. Plant colonization is tightly linked to many soil characteristics (e.g., texture, structure, depth, organic fraction, and nitrogen¹⁹).

Herbivory in glacier forelands might also affect primary succession. Experimental studies have reported a positive short-term effect (i.e., 7–10 years) of herbivory on plant community assembly and ecosystem functioning on bare landslide surfaces (with *Trichosurus vulpecula*²⁰) and on Taiga river floodplains (with *Alces alces*²¹). First, herbivory affected plant species composition that, in turn, accelerated the successional trajectory. In addition, browsing decreased the biomass of the dominant species and therefore provided greater opportunities for colonization and access to nutrients by other species (e.g., nitrogen-fixing species). The effects of herbivory on plants, however, has a complex periodicity in which biomass increases in seven to ten years but there is a negative impact on the biomass of the colonizing pioneer species from one to three years²⁰.

Second, mammal herbivores enrich soil N and C and soil microbial communities because their wastes increase plant biomass, especially of nitrogen-fixing species²⁰,²¹. Mammal herbivores also feed upon and transport mycorrhizal fungi²². For example, interactions between fungi and large ungulates facilitated primary succession at Mount St Helens²³. Long term (> 30 years) ungulate herbivory can lead to the invasion of unpalatable species, but it can also increase the presence of highly palatable species²⁴. In addition, mammals might play a role in exo- and endo- seed dispersion²⁵. Mammals found in polar and alpine environments, including muskox, lemming, fox, and ermine are potentially significant dispersal vectors in glacier forelands¹⁸. Finally, trampling can also affect vegetation positively or negatively by improving seedling recruitment by pressing seeds into the uppermost soil layer²⁶ or by increasing lichen clonal growth through fragmentation and spreading processes²⁷.

In Andean alpine proglacial environments, Andean camelid dung piles might create more favorable conditions such as resources-rich substrates, sites for seedling establishment, and seed sources, for vegetation establishment (*Vicugna vicugna*²⁸, *Lama guanicoe*²⁹). Reider and Schmidt²⁸ suggested that vicuña dung piles shortcut a 100+ year lag between glacier retreat and primary succession in the Cordillera Volcanota, Peru. Other studies in different ecosystems have highlighted soil nutrient accumulation by alpaca³⁰ or vicuña³¹ dung piles.

Besides soil nutrient enrichment, trampling by large ungulates might also influence soil formation positively or negatively by creating biogeomorphic changes. First, trampling intensity affects mycorrhizae diversity, which affects soil fertility³². Second, trampling can favor pedogenesis by incorporating litter into the soil, enhancing microbial processes and likely influencing the fate of litter carbon, soil organic carbon, and nitrogen mineralization³³. Third, by causing bioturbation, trampling can expose more mineral soil to erosion and generate finer particles³⁴. Moreover, herbivore trampling can compact soils, which might affect soil infiltration, erosion, and nutrient mineralization (e.g., in arid environment³⁵). Nevertheless, camelid-versus cattle-trampling creates less compaction³⁶ and can improve ecosystem functioning through water redistribution (e.g., in hilly dryland³⁷).

Globally, experimental approaches to study pedogenesis and vegetation primary succession after glacier retreat are crucially missing. A few studies implemented experimental research designs in glacier forelands, which mainly included seed burial³⁸ and experimental warming²⁹. Moreover, all studies took place in the European Alps, while pressing ecological challenges related to deglacierization exist globally⁴⁰. While glaciers melt worldwide, emerging proglacial ecosystems are becoming important components of the freshwater and carbon cycles. They might also act as sediment sinks that buffer upstream hazards and be refugia for cold-adapted species⁴¹. The establishment and organization of novel proglacial systems is a key concern for adaptation to glacier retreat. Thus far, we have no study that includes experiments with large herbivores, including camelids, in deglacierized landscapes. However, experimental approaches are necessary to understand more precisely the modalities of proglacial ecosystem development and favor adaptation⁴.

In this study, we used an inclusion experiment to understand if the presence (i.e., herbivory and trampling effects) of llamas (*Llama glama*) in glacier forelands can enhance primary succession after glacier retreat. We evaluated the effects of llamas on (1) soil pedogenesis and (2) primary vegetation succession in a proglacial ecosystem of the Tropical Peruvian Andes, where glaciers are exceptionally vulnerable to global warming⁴². First, we hypothesized that the llama dung piles would drive soil enrichment, catalyzing soil formation processes; soils would have higher nitrogen and carbon contents in the presence of llamas. Second, we hypothesized that the presence of llamas enhances the primary succession of the vegetation. With llamas, we expected an increase in plant biomass, as well as higher nitrogen and phosphorus foliar contents because of soil enrichment. Finally, we hypothesized that the llamas would function as seed transporters between lowlands, highlands, and valleys, by carrying seeds in their hooves, wool (epizoochory), or stomach, and stimulate seed germination via endozoochory. Our study aims to replicate a natural phenomenon to enhance ecosystem formation and inform conservation and/or land management strategies in proglacial landscapes.

# Materials and methods

## Study area

We conducted our experiment in the Uruashraju glacier foreland located in the Quebrada Pumahuacanca within the National Park of Huascarán (NPH), Cordillera Blanca, northwestern Peru (Fig. 1). The Cordillera Blanca is the world's most extensively glacier-covered tropical mountain range⁴³. In the 1960s, the imminent extinction of the overhunted and emblematic vicuña, a native Andean camelid, prompted the creation of the NPH. This area is semi-arid and has highly seasonal precipitation with 80% of the 700 mm (~3500 m.a.s.l.) to 1000 mm (~4550 m.a.s.l.) per year falling between October and May⁴⁴. The mean annual air temperature at 3450 m.a.s.l. is approximately 12–14 °C with strong diurnal variability⁴⁵ but little seasonal variability due to its tropical location
<figure><img src="image_3.png" /></figure>

Source: Base layers originated from © Esri, © NASA, © NGA, © USGS, and © FAO; all other layers were created and edited by the authors.

<caption>**Figure 1.** Location and study site set-up. Map of location with respect to the Santa River watershed (a), and Río Negro sub-watershed (b). Map of the experiment within the Uruashraju glacier foreland (c). The glacier retreat outlines were produced and provided by the ANA (Área de Evaluación de Glaciares y Lagunas, Autoridad Nacional del Agua, Huaraz) based on topographic field surveys of the glacier fronts since 1948, and analysis of photographs. Maps generated by authors with licensed software ArcGIS Pro 3.0.2 (https://www.esri.com/en-us/arcgis/products/arcgis-pro/).</caption>

(9°S). The geology of the Quebrada Pumahuacanca is dominated by granodiorite and tonalite, with outcrops of the meta-sedimentary Jurassic Chicama formation, that consists of shale, pyritic siltstone, and quartzite with volcanics in some areas⁴⁶. Sulfide-rich lithologies occur within the glacier foreland and their oxidation after the retreat of the ice is linked to low pH values⁴⁷.

The Pumahuacanca Quebrada (from 4050 to 5605 m.a.s.l.) has patches of high Andean woodlands (e.g., *Polylepis spp.* and *Gynoxis spp.*), wetland ‘bofedales’ dominated by cushion-forming plant communities (e.g., *Distichia muscoides*), and grasslands dominated by tussock grass such as *Cinnagrostis rigida*⁴⁸. Certain species like *Werneria nubigena* indicates the presence of non-native livestock and overgrazing. Within the NPH, Pasture User Committees (CUP) use these high pasture ecosystems to graze their livestock. This is the case of the CUP of Arwaycancha in the Quebrada Pumahuanca. In recent years, the uncontrolled introduction of exotic livestock for agropastoralism and tourism activity has caused degradation of pastures and the disappearance of some native species within the Park⁴⁸. Due to their soft footsteps and prehensile split upper lips Andean camelid minimize impact on the groundcover compared with introduced livestock⁴⁹,⁵⁰. The NPH has thus supported camelid breeding initiatives as an ecosystem and biodiversity conservation strategy⁵¹.

The study site (77° 19' 20" W, 9° 35' 47" S) lies between 4680 and 4700 m.a.s.l. in a terrain deglacierized between 1979 and 1995² (Fig. 1). This research was in collaboration with the NPH and the *Llama 2000 Asociación*, a local community of farmers, who own the llamas. The farmers' own village of Canrey Chico, Olleros has contaminated water due to acid rock drainage (ARD) that threatens their agriculture⁵³. As an adaptation strategy to glacier retreat, the farmers aim to find solutions to increase the vegetation cover in recently deglacierized areas in order to reduce ARD. Previous floristic plot surveys carried on in the Uruashraju glacier foreland identified 47 vascular plant species (Zimmer et al., unpublished data; Table S3). We selected the Uruashraju foreland based on the *Llama 2000 Asociación*'s shared interest in our research questions, the presence of llamas in the adjacent Quebrada, data availability for deglacierization chronology, floristic comparisons within the overall glacier foreland (Table S3), and the ecological representativeness of the site. Considering the socio-geoecological aspects of the region, we expect our results to be generalizable to the overall Cordillera Blanca and, possibly, to other mountain ranges.
## Context and experimental design

The *Llama 2000 Asociación* launched the *Llama 2000* Project to increase the value of local llama breeding and conserve an essential part of their Inca biocultural heritage. In 2019, the community aimed to strengthen its activity promoting the sustainable management of llamas and vicuñas to develop community-based tourism along the Qhapaq Nan trail, enhance local economy, and develop climate change adaptation strategies. We collaborated with the *Llama 2000* community and the NPH to test a novel adaptation strategy to glacier retreat by introducing llamas in the Uruashraju glacier foreland. Manipulative experiments better predict the causal impacts of short-term ecological change and allow for a better control of confounding factors⁵⁴, which here included livestock disturbance, llama stocking rate and visit frequency to the site, and camelid movement between Quebradas. Here, we are testing a novel adaptation strategy to glacier retreat by introducing llamas in a glacier foreland in a region where the density of Andean camelids has strongly decreased over the last century and where camelid breeding initiatives have been strengthened by National Entities (i.e., NPH and SERNANP—Servicio Nacional de Áreas Naturales Protegidas por el Estado⁵¹).

In June 2019 we set up a complete block design with four monitoring blocks—i.e., four llama plots and four control plots (Fig. 2a, b)—across a chronosequence since deglacierization: the south corner of the experiment was ice-free in ~1979 and the north corner in 1995 (Fig. 1c). We built eight plots (925 m² each) over 12,500 m² of low to moderate slopes. Within each plot, we randomly established eight (1m²) permanent subplots (Fig. 2c). Thanks to personal communications with the *Llama 2000* farmers, specialists of the NPH, and a local zootechnician, we determined the animal load and grazing intensity within the plots based on (1) the natural density of vicuña observed in proglacial landscapes in the Cordillera Blanca (NPH, *Unpublished data*), (2) llama and vicuña behavior (herd grouping and displacement between altitudinal range), and (3) vegetation cover of the study site and llama nutritional requirements. These considerations allowed us to replicate a natural phenomenon (i.e., several proglacial landscapes have been identified as vicuña habitat within the NPH). From June 2019 and to 2022, in each fenced plot, three llamas grazed for three days each month (cf. Appendix 1—Experimental design).

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Experimental design and in situ surveys. Design of the experiment (a), *Llama glama* within a llama plot (b), 1m² vegetation subplot (c), seedling germinated from llama feces found within the experiment in June 2022 (d).</figcaption></figure>
## Field data collection

To test the effects of llamas—grazing, trampling, and latrine behavior—on proglacial ecosystem development following our three hypotheses (i.e., soil enrichment, vegetation succession, and seed dispersal), we collected: (1) pedological, (2) plant diversity and productivity data, and (3) llama dung pile samples. The data collection took place at different times (Table 1).

First, in May 2019, before introducing the llamas, we conducted the first floristic and geomorphic evaluations of the 64 subplots (i.e., Hypothesis 2). At each subplot, we identified all the vascular species (Table S4). We estimated visually their relative surface cover, density, necromass, and their average height. As well, we recorded the presence of biological soil crusts (BSC) measuring their heights and relative covers. In each subplot, we also estimated soil granulometry —percentage cover of sand (<2 mm), gravel (<2 cm), rock (<10 cm), and block (>10 cm). Every year from 2019 to 2022, we replicated the same biotic and abiotic data collection for each of the 64 subplots to study the effects of the presence of llamas over time. In addition, in 2020, 2021, and 2022 we visually assigned a greenness index (0, 1, or 2) to each subplot based on the proportion of green vegetation within each subplot. During the field evaluation in 2020 and 2021 we observed the presence of eaten pasture within the llama plots, which could have biased our estimations of plant cover and height. To overcome this concern, in 2022 we waited three weeks between the last grazing treatment and our field evaluation, giving time for the plants to grow again after the grazing event. We also reported eaten pasture within the control subplots in 2020 and 2021, indicating either that the llamas escaped at least once from the fences, despite attentive care, or the llamas grazed within the control area during their arrival on site while the farmers were distributing the animals between plots.

In addition, we measured three plant functional traits: leaf dry matter content (LDMC), total nitrogen (N), and phosphorus (P) to study plant productivity (i.e., Hypothesis 2). Plant leaf N is directly related to photosynthesis and respiration and reflects the trade-off between greater photosynthetic capacity and potentially suffering more herbivory⁵⁵. LDMC correlates negatively with potential relative growth rate and positively with leaf lifespan⁵⁶. Leaves with high LDMC tend to be relatively tough and are thus assumed to be more resistant to physical hazards (e.g., herbivory) than are leaves with low LDMC⁵⁵. Therefore, LDMC is an indicator of a plant's resource use strategy: usually, stress tolerant and slow growing species found in stressed environments display high LDMC and low N, whereas low LDMC and N-rich leaf are characteristic of competitors and ruderals with fast resource acquisition strategies⁵⁵. For the trait measurements, we selected four species based on (1) their presence in the 8 plots (Table S4), (2) their potential as ecosystem engineers (i.e., *Pernettya prostrata*: slope stabilization; *pers. obs.*, *Senecio sublutescens*: high litter production and micro-habitat creation; *pers. obs.*), and (3) their palatability (*Agrostis. tolucensis*; *P. prostrata*; *pers. obs.*).

In June 2022, we collected 10 leaf samples from different individuals of the same species within each plot. Each sample is a composite of a minimum of three individuals, summing to approximately 15 to 40 leaves per sample depending on the species. We used stratified sampling within the four llama plots, collecting half of the 10 samples from individuals located inside or down-slope of a dung pile and the other five samples from individuals located three meters outside of dung piles (cf. Appendix 1 – Field data collection). The collection of plant material and soil samples was conducted in accordance with relevant recommendations and authorizations of the SERNANP (licenses N°008–2019-SERNANP-JEF and N°005–2022-SERNANP-JEF). The soil samples were exported to the USA under the USDA (United States Department of Agriculture) permit to receive soil P330-17-00,299 and P330-20-00161_20200706. All methods were carried out in accordance with relevant guidelines. As importantly, because we measured and counted the same plant specimens present within the permanent subplots in the control and the llama plots every year, we did not collect the full plant specimens; the collection of plant specimen would have been destructive for the experiment and biased the floristic evaluations results for this study but also for the long-term monitoring of the experiment. Therefore, the plant species identified and studied have not been deposited in a public herbarium, since the registration process requires the full specimen (including vegetative material and root). Accordingly, the plant identification was done in the field by botanists from the Museum of Natural History of the National University of San Marcos (UNMSM), Lima, Peru (Sebastián Riva and Jean Salcedo). Regarding the measurement of the plant traits (nitrogen, LDMC, phosphorus), we only collected the minimum amount necessary for the analysis, collecting only the leaves directly from the plants in the field, and trying to be as minimally destructive as possible for the experiment.

Then, to infer the effect of the llamas on soil development (i.e., Hypothesis 1) in May 2019 (before introducing the llamas) and June 2022 (three years after initiation), we collected three soil samples distributed in the upper, middle, and lower slopes of each llama plot. Each sample is itself a composite of three subsamples collected within 1m², between 3 to 15 cm deep. We excluded from the sampling any pure organic surface layer, litter, and vegetation (cf. Appendix 1 – Field data collection).

Additionally, we collected data of covariate variables that might potentially affect our response variables: subplot slope degree and soil temperature. To measure soil temperature, we buried two HOBO 8 K Pendant data

<table><caption>**Table 1.** Data collected during the field evaluations of the experiment.</caption><thead><tr><td>Data collected</td><td>May 2019 (dry season)</td><td>December 2020 (wet season)</td><td>May 2021 (dry season)</td><td>June 2022 (dry season)</td></tr></thead><tbody><tr><td>Floristic & geomorphic surveys</td><td>X</td><td>X</td><td>X</td><td>X</td></tr><tr><td>Plant functional traits</td><td>X</td><td></td><td></td><td>X</td></tr><tr><td>Soil sampling</td><td>X</td><td></td><td></td><td>X</td></tr><tr><td>Llama dung pile sampling</td><td></td><td></td><td></td><td>X</td></tr></tbody></table>
loggers (five to 10 cm depth) per plot within the upper and lower parts of the plot. Since June 2019, the 24 data loggers have recorded temperature every four hours. The slope and soil temperature measurements allowed us to account for confounding factors such as a slope effect on soil composition (i.e., Hypothesis 1) or vegetation cover (i.e., Hypothesis 2) or an effect of the proximity of the glacier on soil temperature. We used the available precipitation and air temperature data produced by the ANA (Autoridad Nacional del Agua, Huaraz, Peru) from a weather station installed within the Uruashraju glacier foreland) to corroborate our soil temperature analysis and infer the role of seasonality on the response variables (Figure S2).

Finally, to estimate the fertilization effects of the llamas on the soil composition and their role as seed dispersers (i.e., Hypothesis 3; Fig. 2d), we collected three fecal samples from three different llama dung piles in each llama plot. The Environmental Quality Laboratory (EQL) in Huaraz, Peru analyzed samples for nutrients, and the Museum of Natural History of the UNMSM in Lima, Peru analyzed seeds in the samples (cf. Appendix 1 – Field data collection).

## Laboratory and herbarium method

### Laboratory analysis

To infer the role of llamas on pedogenesis, we analyzed soil pH, particle size distribution, soil organic (SOC) and inorganic (IC) carbon, soil nitrogen (N), and carbon isotopes (δ¹³C). We analyzed the 48 soil samples at the University of Texas at Austin, USA, mostly at the Soils and Geoarchaeology Laboratory. Soil pH was determined using a soil:solution ratio of 1:2. Because estimating SOC in proglacial soil is complex (due to low SOC and carbonate or clay contents), we used two different methods: sequential Loss On Ignition (LOI) and Elemental Analyzer (EA). First, we analyzed the 48 samples from 2019 and 2022 using the LOI method for SOC and IC (cf. Appendix 1 – Laboratory analysis). Second, for the 24 samples collected in 2022, we measured SOC and ¹³C by isotope ratio mass spectrometry. Values of δ¹³C provide information about the origin of the SOC. Bardgett et al.⁵⁷ linked high δ¹³C values to the presence of ancient carbon, and Kielland and Bryant²¹ used δ¹³C values to examine the effects of herbivory on soil carbon composition. We measured N using an EA CHNS-O. Finally, at the Laboratory for Environmental Archaeology of UT Austin, we performed Particle Size Analysis (PSA) after deflocculation (cf. Appendix 1 – Laboratory analysis).

The EQL in Huaraz, Peru measured leaf LDMC, N, and P, as well as the total organic C, total N, and total P of the fecal samples. We measured LDMC only for *P. prostrata* and *S. sublutescens* since we did not expect major changes within the tussock grasses⁵⁸. N and P were measured for all species. At UNMSM in Lima, we evaluated the viability of the seeds using the Tetrazolium test (Table S22; cf. Appendix 1 – Laboratory analysis).

## Data analysis

Data analyses were performed using the statistical software R version 4.2.2. To test our first hypothesis (i.e., llamas catalyze soil formation), we first computed summary statistics for all the soil response variables (texture, pH, SOC, IC, N, and δ¹³C) for the control and llama plots. Afterward, we used linear models to test the significance of the variation assignable to the llama treatment. We included as response variables texture, pH, SOC, IC, N, and δ¹³C. We used a fixed blocking factor (4 blocks) to control for time constant unobserved effects within each block (Fig. 2), and “Treatment” and “Years” as fixed factors. We first ran a model based on the experiment’s full three-years (including samples from 2019 and 2022). Second, we ran two separate models, one for the 2019 data and one for the 2022, to remove the variable “Years” from the model to gain degrees of freedom and reduce standard errors. Since the soil nitrogen data from 2019 and 2022 are not comparable because of different drying methods, we did not run the full 3-years model on soil nitrogen. We report the results of all models within the results section and supplementary material.

To test our second hypothesis, we assessed plant community changes based on species richness—count of species per subplot to estimate the taxonomic diversity of the subplots, plant cover, and plant functional traits (i.e., LDMC, N and P, plant height, necromass, and greenness) between the control and llama plots. Plant height and necromass were community weighted for each subplot. As well, we calculated a percentage change in plant cover between 2019 and 2022. First, we computed summary statistics. Then, as for the soil data, we used linear models with a fixed blocking factor to analyze the effect of llamas on plant cover, height, and necromass. The linear models for the vegetation response variables included “Treatment”, “Years”, and “Block” as explanatory variables. In addition, the initial plant cover values reported in 2019, were used in the plant cover analysis as a co-variate to increase the power of the test. We computed the total percentage cover of each species for the overall experiment to select the four species most representative of the plant community. Then we ran distinct models for the four species to test for the effects of the llamas on the species selected. We used Generalized Linear Models with a Poisson distribution to test the effects of the llamas on plant richness and greenness because both variables were discrete, and we hypothesized that species richness patterns were scattered during the first stage of primary succession (i.e., species not aggregated). To evaluate differences between years and between control and llama treatment, when relevant (i.e., Plant cover and Greenness), we computed pairwise comparisons using the *emmeans* package with Bonferroni corrections (*emmeans* function). Afterward, we used linear regressions to assess the effect of the initial subplot cover on the percentage change in cover between 2019 and 2022 and the 2022 plant cover. We log transformed the variable Plant Cover to increase the fit of the model (Shapiro-Wilk normality test: p-value = 0.094, W = 0.988; Breusch-Pagan test: p-value = 0.46; Chi-square = 0.547). We computed all the linear and generalized linear models using the *lm* and *glm* functions of the *stats* library and the model fit diagnostics were checked using the *DHARMa* package.

To analyze variation in soil temperature we used the period with complete data from 23.06.2019 to 22.05.2021. We calculated summary statistics and used Wilcoxon tests to assess the difference between control and llama treatments. To infer differences in air temperature that might have affected our response variables we compared
mean daily temperatures over three periods corresponding to our field evaluations dates: June 2018 to May 2019; June 2019 to May 2021, and June 2021 to May 2022. We used ANOVA and Tukey Honestly Significant Difference tests to assess differences between the three periods. Finally, we tested for significant differences in slope gradient and granulometry of the subplots between the llama and control plots with T-tests and Wilcoxon tests (ggpubr library).

# Results

## Environmental characteristics of the plots

In 2019, prior to llama manipulation, there was no significant difference between the mean slope of the control and llama plots. Similarly, there was no significant difference in sand, gravel, rock, and block contents between the llama and control plots (Table S5). Between June 2019 and May 2021, we observed lower daily minimum (−0.45 °C; *p-value* = 3.8e-10) and higher daily maximum temperatures (+3.88 °C; *p-value* = 0.001) within the llama plots than within the controls (Table S5). The air temperature data (from June 2018 to May 2022) showed significantly lower temperatures between June 2021 and May 2022 than between June 2018 and May 2019 (diff = −0.22, *p-value* = 0.02), and between June 2019 and May 2021 (diff = −0.31, *p-value* < 0.001; Figure S2; Table S2). The mean daily precipitations between 2017 and 2020 displayed differences in precipitation between the dry season (June to September) and the wet season (October to May)⁴⁴.

## Effect of llama presence of soil development

After 3 years of treatment, we observed an effect of the llamas on the silt and sand contents (Figure S3). Both the full three-years model (2019 and 2022 data) and the 2022 model showed an increase in silt within the llama plots (respectively *p-value* = 0.017; Table S7 and *p-value* = 0.059; Table 2). Only the full three-years model showed a significant decrease in sand within the llama plots (*p-value* = 0.019; Table S7). Between 2019 and 2022, SOC and IC obtained by the LOI method and N were affected positively by the llamas (respectively + 0.2% SOC, *p-value* = 0.047; + 0.1% IC, *p-value* = 0.047 and + 0.05% N, *p-value* < 0.001; Fig. 3, Table 2, and Table S6). The SOC data from the EA method also showed a significant effect of the llamas on the SOC in 2022 (+ 0.16% SOC, *p-value* = 0.047). The δ¹³C signature of SOC in 2022 did not differ significantly between the treatments as well as pH (Table 2). The organic C of the dung pile (15.4%) was approximately 9 times higher than the SOC of the control plots in 2022 (1.74%). Similarly, dung pile total N (5.07%) was 36 times higher than the N in the 2022 control plots (0.14%; Table S11). We did not detect a leaching effect from the upslope to the downslope portions of the experimental plots, because there was no significant difference in soil composition between the samples position on slope.

## Primary succession of vegetation

### Vegetation biomass

In 2019, the mean subplot cover of the llama (7.96%) and control (11.38%) plots were not significantly different (*p-value* = 0.12; Table S12). Between 2019 and 2022, the presence of llamas had a significant (*p-value* < 0.001) positive effect on plant cover (Table S13). Between 2020 and 2021, both the llama (*p-value* = 1) and the control (*p-value* = 1) plots did not show significant increases in vegetation cover (Fig. 4; Table S14). In contrast, from 2021 to 2022, the mean cover in the llama plots increased from 8.91% to 13.97% (+57%, *p-value* < 0.001; Fig. 4; Table S12 and S14). Linear regressions to model the relationship between the percentage change in cover between 2019 and 2022 and the 2019 subplot cover showed that the control subplots with higher cover in 2019 have gained less cover during the three years than the subplots with lower cover in June 2019 (*p-value* = 0.004, R-squared = 0.24). The trend is lowered within the llama subplots (*p-value* = 0.018, R-squared = 0.18; Figure S4). We did not detect a significant effect of the llamas on the BSC.

Across the four surveys of the experiment, we observed a total of 14 plant species. Of these, 13 were identified at the species level (Table S4). Results from the generalized linear model showed no effect of the llama treatment

<table><caption>**Table 2.** Significance of the effects of the llamas on soil properties in 2019 and 2022: texture (clay, silt, and sand), pH, organic carbon (SOC) by LOI and EA methods, inorganic carbon (IC), nitrogen (N), and δ13C. Bold values indicate significance (p≤0.05). The full coefficients of the models are reported in Table S8 and S9 in Supplementary material.</caption><thead><tr><th rowspan="2"></th><th colspan="4">2019</th><th colspan="4">2022</th></tr><tr><th>Est</th><th>SE</th><th>t value</th><th>Pr(&gt;t)</th><th>Est</th><th>SE</th><th>t value</th><th>Pr(&gt;t)</th></tr></thead><tbody><tr><td>Clay</td><td>4.498</td><td>2.084</td><td>2.158</td><td><strong>0.044</strong>*</td><td>0.601</td><td>2.511</td><td>0.239</td><td>0.813</td></tr><tr><td>Silt</td><td>1.945</td><td>1.379</td><td>1.411</td><td>0.174</td><td>5.406</td><td>2.696</td><td>2.005</td><td>0.059</td></tr><tr><td>Sand</td><td>-6.441</td><td>2.939</td><td>-2.191</td><td><strong>0.041</strong>*</td><td>-5.995</td><td>4.450</td><td>-1.347</td><td>0.194</td></tr><tr><td>SOC by LOI</td><td>-0.015</td><td>0.094</td><td>-0.155</td><td>0.878</td><td>0.205</td><td>0.096</td><td>2.128</td><td><strong>0.047</strong>*</td></tr><tr><td>SOC by EA</td><td></td><td></td><td></td><td></td><td>0.158</td><td>0.075</td><td>2.122</td><td><strong>0.047</strong>*</td></tr><tr><td>δ13C</td><td></td><td></td><td></td><td></td><td>-0.233</td><td>0.136</td><td>-1.716</td><td>0.102</td></tr><tr><td>N by EA</td><td>-0.009</td><td>0.008</td><td>-1.158</td><td>0.261</td><td>0.048</td><td>0.011</td><td>4.385</td><td><strong>&lt;0.001</strong>***</td></tr><tr><td>pH</td><td>0.093</td><td>0.069</td><td>1.354</td><td>0.192</td><td>0.250</td><td>0.163</td><td>1.535</td><td>0.141</td></tr><tr><td>IC</td><td>-0.007</td><td>0.044</td><td>-0.155</td><td>0.878</td><td>0.097</td><td>0.045</td><td>2.128</td><td><strong>0.047</strong>*</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Effects of llamas on soil organic carbon (SOC by LOI) and nitrogen (N) (95% confidence intervals) in 2019 and 2022. Means not sharing any letter are significantly different by Bonferroni tests ($p \le 0.05$).</figcaption></figure>

on species richness between 2019 and 2022 (*p-value* = 0.534), whereas there was a significant effect of the year 2022 (*p-value* = 0.015; Table S13). Nevertheless, the *Belloa piptolepsis*, *Oritrophium limnophilum*, *Senecio rufescens*, and an unidentified juvenile individual (*sp*. ) appeared only in the llama plots in 2022. Also, *Melpomene peruviana* appeared only within the control plots in 2022 but was already present in 2019 (Table S4). The four species with the highest total cover were *Cinnagrostis rigida*, *Senecio sublutescens*, *Minoides kunthiana*, and *Agrostis toluscensis* (Table S4). Linear models ran separately for the four species showed that the llama presence had a higher effect on the plant covers of *S. sublutescens* (β = 0.45; *p-value* < 0.001) and *M. kunthiana* (β = 0.36; *p-value* = 0.01) than on *A. toluscensis* (β = 0.11; ns) and *C. rigida* β = 0.12; *p-value* = 0.02) (Table S16; Figure S5).

### Plant traits

Plant traits analysis revealed that there is no effect of the llamas on the LDMC of *Pernettya prostrata* (*p-value* = 0.06). However, we observed higher LDMC within the samples collected within the llama plots but far from dung piles (i.e., Llama: no dung pile) than the samples collected in the control plots (+ 28.7 mg.g⁻¹) and samples with dung pile influence (+ 32 mg.g⁻¹; Figure S6; Table S17). There is a positive effect of the llamas on the LDMC when there is no dung pile influence (*p-value* = 0.016; Table S18). We did not observe significant changes in LDMC, N, and P for the other species sampled. Subplot greenness was positively affected by the llamas (Table S13). There were significant decreases in plant greenness between 2020 and 2021 in the control plots (*p-value* ≤ 0.001) and the llama plots (*p-value* ≤ 0.001; Fig. 5; Table S15). Although there was no significant increase in greenness between 2021 and 2022 for both treatments, we observed higher greenness index within the llama plots than within the control plots in 2021 (*p-value* = 0.002) and in 2022 (*p-value* = 0.002; Table S15). During the four floristic surveys we have reported several observations of pasture eaten within the subplots in both the control and llama plots (Table S20). We did not detect an effect of the llamas on plant necromass and height.

### Seed and germination potential

In total, 266 seeds were found within the llama fecal pellets' samples (Fig. 6). Seeds of 12 plant species belonging to seven families were identified within the pellets, from which two families (Poacea and Ericaceae) were also
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Effects of llamas on plant cover (95% confidence intervals). Letters represent the significantly different groups according to the post hoc contrasts (p ≤ 0.05).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Effects of llamas on greenery (95% confidence intervals). Letters represent the significantly different groups according to the post hoc contrasts (p ≤ 0.05). Note: December is the wet season.</figcaption></figure>

present within the llama plots (Table S21 and S22). The viability tests showed that five of the seed species identified still maintain their germinative power between 7.4% and 70% (Fig. 6; Table S21).
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Seed viability per species found within the llama pellets (n = 266). The percentages expressed are the percentages of viable seeds for each species. For the taxa identified at the species or genus level family names are in brackets.</figcaption></figure>

Discussion

Our study investigates the effects of llamas on pedogenesis and vegetation primary succession after glacier retreat. We expected that the presence of llamas and their latrine behavior would catalyze soil development and enhance vegetation succession. Unlike previous studies on camelid herbivory in proglacial landscapes, our results demonstrate that in addition to drive soil enrichment and increase plant cover, the llamas have transported viable seeds via endozoochory to the proglacial landscape, likely facilitating dispersal processes.

The analysis of the soil samples showed that over three years, llamas had a positive impact on proglacial pedogenesis through a modification of the soil texture and increases in SOC and N. First, the very high SOC, N, and P measured in the dung pile samples demonstrate that Andean camelid dung piles act as nutrient hotspots within the proglacial landscape²⁹. Albeit we do not disregard limitations due to the lower temporal resolution of the dung pile analysis, the increase in SOC and N observed in the llama plot soil samples are likely due to a translocation of nutrients via surface runoff from the dung piles to the adjacent areas³⁰. Andean camelids graze in dung pile areas because of their nutritional benefits²⁸,³¹. Therefore, based on our results we hypothesize that in a longer term there would be also an intra-systems redistribution of the nutrients within the proglacial pasture from old to new dung pile sites, in addition to the transfer from downstream to proglacial systems. Besides dung pile nutrient addition, camelid trampling might also facilitate litter and BSC incorporation to the subsurface soil, accelerating its decomposition and favoring C and N enrichment³⁴.

Soil texture can serve as proxy for proglacial pedogenesis. Proglacial chronosequence studies have reported a reduction of the fraction of sand and an increase in silt with increasing age¹²,⁵⁹. We suggest that here the observed increase in silt is due to llama bioturbation processes³⁴ that favor silt incorporation. Second, the increase in plant cover would have reduced erosion processes and increased the capture of recent, carbonate-rich sediment within the denser vegetation, creating feedback effects of the vegetation on soil development (i.e., biogeomorphic interactions⁶⁰). These texture changes reveal an accelerated soil development with the llamas, which would lead to changes in several soil characteristics and microtopography changes, and affect the spatial redistribution of water and soil resources, as well as proglacial ecosystem functioning and productivity³⁷,⁵⁹. Future analysis will investigate nitrogen isotope signature and microbial activity (e.g., mycorrhizae), to improve our understanding of how the llamas affect belowground communities and therefore the overall nutrient cycles in recently deglacierized landscapes.

The δ¹³C signatures of the llama and control plots are consistent with previous observations in glacier forelands and demonstrate the presence of ancient carbon⁵⁷. The slight depletion (but not significant, *p-value* = 0.1) in δ¹³C between the control and llama plots could support an enrichment in modern carbon in line with the observed increase in SOC. Indeed, δ¹³C is higher in low-carbon soils and decreases as SOC increases due to recent organic carbon inputs¹⁴.

After three years, the presence of llamas had a substantial impact on the primary vegetation succession at the Uruashraju glacier foreland, in line with hypothesis 2. The optimization of aboveground net primary production by mammals is reported in grass and shrubland systems²⁴,⁶¹, and few studies also reported positive pathways in arctic tundra⁶²,⁶³. Here, the increased plant cover and overall greenness of the fenced plots imply that the llamas have increased the development of the aboveground vegetation. Although we reported significantly lower minimal and higher maximal soil temperatures within the llama plots in comparison to the control plots, we do not consider that these environmental variations would significantly affect plant growth. The differences are relatively small and occur at temperature ranges that are not critical to alpine plant life¹⁹. Similarly, we recognize a
potential effect of the decreased air temperature between 2018 and 2022, but those differences are relatively small to have a significant effect on vegetation development. We do, however, suggest that differences in precipitations caused the increase in cover observed within the control plots between May 2021 and June 2022. Similarly, we hypothesize that the high greenness value observed in December 2020 in comparison with the May 2021 and June 2022 is due to the highest precipitation occurring during the wet season. Further work will continue to monitor soil temperature and precipitations to infer the role of the llamas on micro climatic characteristics and decipher the role of climatic covariate variables during primary succession.

These changes are likely to have significant effects on the landscape evolution, enhancing vegetation-soil feedback cycles⁵⁹,⁶⁴. For example, *Senecio sublutescens* had the highest cover increase in the llama plots and is very likely to produce significant quantities of litter that will be redistributed to the soil surface through llama trampling and litter decomposition³³. Thus, the effects of the llama latrine behavior on soil pedogenesis will lead to increased plant cover²⁹. The higher llama effects on the plant covers of *S. sublutescens* and *Minoide kunthiana* in comparison to *Agrostis toluscensis* and *Cinnagrostis rigida*, suggests that species abundance might be related to species palatability. Although there is scant literature on these species, *C. rigida* and *A. toluscensis* are known to be palatable for Andean camelid, whereas *Senecio* species are generally not very palatable to livestock⁶⁵, and *M. kunthiana* has a very low stature that might be difficult for llama grazing. Albeit the effects of the llamas on plant cover do remain low after three years, we reported significant changes within the proglacial ecosystem during this short period whereas proglacial ecological dynamics are known to be slow². Indeed, most of experimental studies are carried out for longer periods (i.e., 10 to 35 + years)²⁰,²¹. Thus, our study indicates that biotic and abiotic modifications caused by the presence of the llamas can quickly have significant effects on pedogenesis and primary succession after glacier retreat. Our short-term results and the low temporal resolution of our sampling underline the need to study both short term and long-term effects of camelid herbivory on novel proglacial ecosystems.

Llamas have modified the functional traits of the proglacial plant community. Studies show that plant greenness, leaf thickness, and LDMC are tightly linked to plant growth, chlorophyll production, and correlated with leaf N⁶⁶,⁶⁷. The increase in subplot greenness within the llama plots suggests higher photosynthetic activity and plant growth, which is highly likely connected to the dung pile soil nutrient enrichment. In addition, the higher LDMC values of *P. prostrata* observed within the “llama—no dung pile influence” group in comparison to the “controls” suggests a defense mechanism of *P. prostrata* responding to a grazing effect⁶⁸. The individuals of *P. prostrata* suffering from grazing and not beneficiating from the enriched soil may have lowered their growth rate and invested more resources in structural protection of the photosynthetic tissues, as an adaptation to avoid the effect of defoliation (i.e., avoidance strategy⁶⁸). While the lower values found within the llama plots with “dung pile influence” suggest a nutrient enrichment effect from the dung piles and therefore a reverse response of *P. prostrata* LDMC, lowering the negative consequences of the llama grazing (i.e., tolerance strategy). The llama latrine behavior has locally (i.e., a few meters) transformed a low-nutrient system into a more productive system, and individuals benefiting from the enriched substrate have faster regrowth traits⁶⁹. We hypothesize that variation in climatic conditions (precipitation during the wet versus dry season) explain the higher greenness index observed in December 2020 and the strong decrease in greenness between December 2020 and May 2021.

Our results also suggest a decrease in competition between species within the llama plots. The linear model positive coefficient of the 2019 plant cover means that plots with higher cover in 2019 displayed stronger positive effects from the llamas. In the control plots the 2019 initial cover explained a fourth of the variation of the percentage change in cover between 2019 and 2022; thus we hypothesize that in the absence of llamas, the subplots with low cover are more likely to experience higher rates of plant cover increase because of a potentially reduced plant-plant competition. Whereas in the llama plots the strength of the relationship is lower, indicating that other factors (i.e., the llamas) might have reduced competition through the physical effect of grazing⁶². However, we do not exclude the possibility that abiotic factors, such as precipitation or temperature, might also have influenced plant-plant interactions within the experiment. We also hypothesize that the llamas would have an effect on abiotic facilitation though trampling (e.g., alteration of nurse features such as rock), which are keys associations during the first stage of primary succession⁵. The presence of llamas in proglacial landscapes may modify primary succession trajectories and allow plant communities to overcome the harsh environmental conditions.

We did not observe a notable change in plant richness within the experiment. Although, we registered four new species that appeared only within the grazed areas. Similarly, we noted the disappearance of one species within the llama plots in 2022 that was present in 2019. The disappearance of *C. ovata* might be due to the temporal variability in moisture inputs within the proglacial systems, since *C. ovata* is mainly found in wet terrains²⁰. The precipitation data available and analyzed did not allow us to test differences in precipitation between years. Long-term monitoring is necessary to fully evaluate the effects of the llamas on plant richness and plant functional traits. Finally, future research will need to investigate the feedbacks between plant defoliation from herbivory and its effect on soil biota and all factors on the interactions of mycorrhizae⁷¹.

Ungulates are known to act as dispersal agents through endozoochory, carrying seeds on their coats, or between their hooves (epizoochory), or simply spitting out seeds after mastication or rumination²⁷,²⁷. Here, we show that llamas have the potential to disperse seeds of different plant families and types from lower elevation or neighboring valleys to the proglacial habitat (Table S18). Within the llama dung we found viable seeds from five distinct species: *Alchemilla sp.*, *Luzula sp.*, *Pennisetum sp.*, *P. prostrata*, and *Sporobulus sp.*, of which only the *Pernettya* and *Luzula* were reported during the floristic surveys of the experiment, indicating an external seed input. The seeds of *Pernettya*, *Sporobulus*, and *Luzula*, being viable and endozoochore, evidence the role of Andean camelids as seed transporters by endozoochory⁷⁴. Notably, the *Luzula racemosa*, which was not present in 2019 within the experiment but found in 2022 in the llama plots, is present at lower elevation in the valley (57 years after deglacierization; Table S3) and therefore might have been brought to the experiment by the llamas. Since llamas move regularly between vegetation communities for foraging and ranging, both endo- and exozoochorous dispersal might provide an important potential mechanism for colonization. We recognize that different physical
and chemical treatments of the seed (e.g., digestive treatment, feces humidity, and deterioration) might have produced different viability results or also impacted seed dormancy and therefore germination. In addition, we hypothesize that seeds evaluated as “not viable” during our viability test from the llama pellets, could also have been transported via llamas exozoochory from the lower ecosystems, and thus be present in viable condition within the experimental area. For example, species such as *Plantago tubulosa* and *Werneria pygmaea*—peat-forming cushion plants⁴⁸—and *Alchemilla pinnata*⁷⁴ have been found to be part of the Andean camelid diet, and might be transported via llama endo or exozoochory to the more recently deglacierized terrains. Those species are present within the Uruashraju glacier foreland in terrain deglacierized for more than 57 years (Table S3). In addition, we do not exclude the role of human-directed movement of the llamas and anthropogenic seed dispersion. First, the farmers might have influenced the displacement path and the selection of the pasture area during the daily grazing at lower elevation, which in turn could have affected the identity of the species grazed and the seeds present within the fecal pellets, as well as the seed dispersal by exozoochory. Second, we also recognize the role of humans as dispersal agents from lower elevation to the proglacial area. Finally, we do not exclude the role of the llamas as facilitator of the spread of invasive exotic plants by epizoochoory or endozoochory such as *Rumex acetosella*⁷⁵. In the Andes, the study of archaeobotanical remains have given insights on past landscape uses and on the role of camelid in shaping plant communities over long time scales⁷⁶,⁷⁷. Here, our results demonstrating that llama dung can be a source of viable seeds in deglacierized terrains provides promise for the role of Andean camelids in shaping the emerging and future proglacial landscapes.

In the Andes, the Puna ecosystem was a primary center of ungulate domestication⁴⁹. However, since the Spanish conquest in the 1500s, European livestock (e.g., cattle, horse, donkey, and sheep) have almost completely replaced domestic camelids and ancestral herding practices. Introduced livestock —often associated with overgrazing—have several negative ecological impacts in Andean ecosystems⁴⁹. First, exotic animals can alter the structure and composition of the native vegetation⁷⁵, the hydrologic soil functions³⁶, and can increase erosion by soil compaction and by uprooting plants⁴⁹. Here we showed that the presence of Andean camelid favors pedogenesis and vegetation primary succession that may result in enhanced high Andean alpine biodiversity and novel proglacial ecosystem functioning. By reducing runoff, increasing carbon sequestration, and potentially reducing ARD, the reintroduction or rewilding of native Andean camelid may favor adaptation to glacier retreat while benefiting the local economy. Andean camelids are known for their fiber and meat production, and are also involved in tourism activities and Andean folklore³¹. Rewilding interventions as strategies for biodiversity conservation and ecosystem functioning are novel management alternatives that could be applied to the benefit of proglacial landscapes⁷⁸.

## Conclusion

Our study shows that Andean camelids (*Llama glama*) can enhance soil pedogenesis and primary vegetation succession in glacier forelands. In three years, the latrine behavior of the llamas led to soil organic carbon and nitrogen enrichment and changes in soil texture, suggesting enhanced proglacial pedogenesis pathways. Similarly, we reported a significant increase of the plant cover and the appearance of four species not previously documented. We found that llamas can function as seed transporters bringing seeds from lower elevations or neighboring valleys to the proglacial habitats and possibly facilitating colonization processes, albeit more years might be necessary to also observe an increase in plant richness. These results have at least three implications. First, llamas might help overcome some of the constraints of primary succession after glacier retreat like dispersal limitation, low soil fertility, and plant competition, lowering dispersal and establishment lags⁹ and enhancing ecological connectivity between deglacierized landscapes. Second, our experimental simulation of a natural process suggests that other Andean camelids displaying similar behavior might also enhance proglacial ecosystem establishment (i.e., vicuña, alpaca, and guanaco²⁹). Lastly, our findings may have implications for the conservation and management of novel proglacial ecosystems around the world. Locally, our study supports the goal of national and local efforts in re-introducing Andean camelids from local communities. The experimental protocol designed in this study will allow long-term monitoring to continue to assess the changes observed in our three-year study and discover what changes may happen at longer timescales.

## Data availability

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Received: 15 March 2023; Accepted: 27 August 2023

Published online: 24 September 2023

## References

1. Rounce, D. R. et al. Global glacier change in the 21st century: Every increase in temperature matters. *Science* (80-). 379, 78–83 (2023).

2. Matthews, J. A. *The ecology of recently-deglaciated terrain: a geocological approach to glacier forelands*. (1992).

3. Cauvy-Fraunié, S. & Dangles, O. A global synthesis of biodiversity responses to glacier retreat. *Nat. Ecol. Evol.* 3, 1675–1685 (2019).

4. Ficetola, G. F. et al. Dynamics of ecological communities following current retreat of glaciers. *Annu. Rev. Ecol. Evol. Syst.* 52, 405–426 (2021).

5. Zimmer, A. et al. Time lag between glacial retreat and upward migration alters tropical alpine communities. *Perspect. Plant Ecol. Evol. Syst.* 30, 89–102 (2018).

6. Darcy, J. L. et al. Phosphorus, not nitrogen, limits plants and microbial primary producers following glacial retreat. *Sci. Adv.* 4, 1–8 (2018).

7. Erschbamer, B., Niederfriniger Schlag, R. & Winkler, E. Colonization processes on a central Alpine glacier foreland. *J. Veg. Sci.* 19, 855–862 (2008).
8. Anthelme, F., Cauvy-Fraunié, S., Francou, B., Cáceres, B. Living at the Edge: Increasing Stress for Plants 2 – 13 Years After the Retreat of a Tropical Glacier. 9, (2021).

9. Alexander, J. M. et al. Lags in the response of mountain plant communities to climate change. *Glob. Chang. Biol.* **24**, 563–579 (2018).

10. Callaway, R. M. et al. Positive interactions among alpine plants increase with stress. *Nature* **417**, 844–848 (2002).

11. Walker, L. R., Wardle, D. A., Bardgett, R. D. & Clarkson, B. D. The use of chronosequences in studies of ecological succession and soil development. *J. Ecol.* **98**, 725–736 (2010).

12. Egli, M., Mavris, C., Mirabella, A. & Giaccai, D. Soil organic matter formation along a chronosequence in the Morteratsch proglacial area (Upper Engadine, Switzerland). *CATENA* **82**, 61–69 (2010).

13. Temme, A. J. A. M. The Uncalm development of proglacial soils in the European Alps since 1850. in *Geomorphology of Proglacial Systems* 315–326 (Springer, 2019). https://doi.org/10.1007/978-3-319-94184-4_18

14. Khedim, N. et al. Topsoil organic matter build-up in glacier forelands around the world. *Glob. Chang. Biol.* **27**, 1662–1677 (2021).

15. D'Amico, M. E., Freppaz, M., Zanini, E. & Bonifacio, E. Primary vegetation succession and the serpentine syndrome: The proglacial area of the Verra Grande glacier. *North-Western Italian Alps. Plant Soil* **415**, 283–298 (2017).

16. Mainetti, A. et al. Successional herbaceous species affect soil processes in a high-elevation alpine proglacial Chronosequence. *Front. Environ. Sci.* **8**, 1–15 (2021).

17. Brown, S. P. & Jumpponen, A. Contrasting primary successional trajectories of fungi and bacteria in retreating glacier soils. *Mol. Ecol.* **23**, 481–497 (2014).

18. Breen, K. & Lévesque, E. The influence of biological soil crusts on soil characteristics along a high arctic Glacier Foreland, Nunavut Canada. *Arctic Antarct. Alp. Res.* **40**(2), 287–297 (2008).

19. Körner, C. *Alpine plant life: Functional plant ecology of high mountain ecosystems*. (Springer Nature, 2021).

20. Bellingham, P. J. et al. Browsing by an invasive herbivore promotes development of plant and soil communities during primary succession. *J. Ecol.* **104**, 1505–1517. https://doi.org/10.1111/1365-2745.12624 (2016).

21. Kielland, K. & Bryant, J. P. Moose herbivory in taiga: effects on biogeochemistry and vegetation dynamics in primary succession. *Oikos* **377**–383 (1998).

22. Vašutová, M. et al. Taxi drivers: the role of animals in transporting mycorrhizal fungi. *Mycorrhiza* **29**, 413–434 (2019).

23. Allen, M. F., O'Neill, M. R., Crisafulli, C. M. & MacMahon, J. A. Succession and mycorrhizae on Mount St. Helens. *Ecol. responses* *Mt. St. Helens Revisit. 35 years after 1980 Erupt*. 199–215 (2018).

24. Augustine, D. J. & McNaughton, S. J. Ungulate Effects on the Functional Species Composition of Plant Communities: Herbivore Selectivity and Plant Tolerance Author (s): David J. Augustine and Samuel J. McNaughton Published by: Wiley on behalf of the Wildlife Society Stable URL: http://.J. Wildl. Manage. **62**, 1165–1183 (1998).

25. Pakeman, R. J., Digneffe, G. & Small, J. L. Ecological correlates of endozoochory by herbivores. *Funct. Ecol.* **16**, 296–304 (2002).

26. Eichberg, C. & Donath, T. W. Sheep trampling on surface-lying seeds improves seedling recruitment in open sand ecosystems. *Restor. Ecol.* **26**, S211–S219 (2018).

27. Heggenes, J. et al. Herbivore grazing—or trampling? Trampling effects by a large ungulate in cold high-latitude ecosystems. *Ecol. Evol.* **7**, 6423–6431 (2017).

28. Reider, K. E. & Schmidt, S. K. Vicuña dung gardens at the edge of the cryosphere. *Ecology* **102**, 1–3 (2021).

29. Henríquez, J. M. Influencia de los defecaderos de camélidos sobre el desarrollo vegetal y riqueza de especies en morrenas glaciales, Tierra del Fuego. *Rev. Chil. Hist. Nat.* **77**, 501–508 (2004).

30. Mcgregor, B. A. & Brown, A. J. Soil nutrient accumulation in alpaca latrine sites. *Small Rumin. Res.* **94**, 17–24 (2010).

31. Franklin, W. L. Biology, ecology, and relationship to man of the South American camelids. *Mamm. Biol. South Am.* **6**, 457–489 (1982).

32. Li, Q. et al. Arbuscular mycorrhizal fungal community structure following different grazing intensities in an alpine grassland. *Soil Sci. Soc. Am. J.* **85**, 1620–1633 (2021).

33. Wei, Y. et al. Transformation of litter carbon to stable soil organic matter is facilitated by ungulate trampling. *Geoderma* **385**, 114828 (2021).

34. Wilkinson, M. T., Richards, P. J. & Humphreys, G. S. Breaking ground: Pedological, geological, and ecological implications of soil bioturbation. *Earth-Science Rev.* **97**, 257–272 (2009).

35. Castellano, M. J. & Valone, T. J. Livestock, soil compaction and water infiltration rate: Evaluating a potential desertification recovery mechanism. *J. Arid Environ.* **71**, 97–108 (2007).

36. Yalli, T. B., Naupari, J. A. & Flores, E. R. Llama and Cattle Grazing Effects on Hydrological Function in a High-Elevation Mountain Rangeland. (2020).

37. Stavi, I., Ungar, E. D., Lavee, H. & Sarah, P. Livestock modify ground surface microtopography and penetration resistance in a semi-arid shrubland. *Arid L. Res. Manag.* **23**, 237–247 (2009).

38. Schwienbacher, E. & Erschbamer, B. Longevity of seeds in a glacier foreland of the Central Alps-a burial experiment. *Bull. Geobot. Inst. ETH* **68**, 63–71 (2002).

39. Erschbamer, B. Winners and losers of climate change in a Central Alpine Glacier Foreland. *Arctic, Antarct. Alp. Res.* **39**, 237–244 (2007).

40. Luckman, B. H., Wiegandt, E. & Orlove, B. S. *Darkening peaks: glacier retreat, science, and society*. (Univ of California Press, 2008).

41. Bosson, J.-B. et al. Glaciers and postglacial ecosystems: common goods to protect in the Anthropocene. (2023).

42. Rabatel, A. et al. Current state of glaciers in the tropical Andes: A multi-century perspective on glacier evolution and climate change. *Cryosphere* **7**, 81–102 (2013).

43. Vuille, M. et al. Climate change and tropical Andean glaciers: Past, present and future. *Earth-Science Rev.* **89**, 79–96 (2008).

44. Kaser, G., Juen, I., Georges, C., Gómez, J. & Tamayo, W. The impact of glaciers on the runoff and the reconstruction of mass balance history from hydrological data in the tropical Cordillera Bianca. *Perú. J. Hydrol.* **282**, 130–144 (2003).

45. Schauwecker, S. et al. Climate trends and glacier retreat in the Cordillera Blanca, Peru, revisited. *Glob. Planet. Change* **119**, 85–97 (2014).

46. Cobbing, E. J. The Geology of the Western Cordillera of Northern Peru. *Inst. Geol. Sci.* (1981).

47. Fortner, S. K. et al. Elevated stream trace and minor element concentrations in the foreland of receding tropical glaciers. *Appl. Geochemistry* **26**, 1792–1801 (2011).

48. Polk, M. H. "They Are Drying Out": Social-Ecological Consequences of Glacier Recession on Mountain Peatlands in Huascarán National Park, Peru. (2016).

49. Baied, C. A. & Wheeler, J. C. Evolution of high Andean puna ecosystems: environment, climate, and culture change over the last 12,000 years in the Central Andes. *Mt. Res. Dev.* 145–156 (1993).

50. Pfister, J. A. et al. Grazing behaviour of llamas, alpacas and sheep in the Andes of Peru. *Appl. Anim. Behav. Sci.* **23**, 237–246 (1989).

51. SERNANP. 2018. Plan Maestro del Parque Nacional Huascarán: 2017 - 2021. Lima, Peru: Servicio Nacional de Áreas Naturales Protegidas por el Estado (SERNANP).

52. ANA. 2022. Informe Anual del Área de Evaluación de Glaciar es y Lagunas (AEGL) – Huaraz: Evaluación de la Disponibilidad de Recursos Hídricos en Glaciar es en un Contexto de Cambio Climático en las Cordilleras Blanca, Central, Huaytapallana, Vilcanota, Ampato Y Vilcabamba.
53. Zimmer, A., Rodríguez, M. K. B., Oropeza, C. J. A., León, J. W. S. & Bullard, J. R. Implementación de Dos Sistemas de Biorreme-
diación como Estrategia para la Prevención y Mitigación de los Efectos del Drenaje Ácido de Roca en la Cordillera Blanca Perú.
*Rev. Glaciares y Ecosistemas Montaña* **4**, 57-76 (2019).

54. Yuan, Z. Y. et al. Experimental and observational studies find contrasting responses of soil nutrients to climate change. *Elife* **6**, e23255 (2017).

55. Díaz, S. et al. The global spectrum of plant form and function. *Nature* **529**, 167-171 (2016).

56. Smart, S. M. et al. Leaf dry matter content is better at predicting above-ground net primary production than specific leaf area.
*Funct. Ecol.* **31**, 1336-1344 (2017).

57. Bardgett, R. D. et al. Heterotrophic microbial communities use ancient carbon following glacial retreat. *Biol. Lett.* **3**, 487-490 (2007).

58. Cruz-Maldonado, N. et al. Aboveground-trait variations in 11 (sub)alpine plants along a 1000-m elevation gradient in tropical Mexico. *Alp. Bot.* **131**, 187-200 (2021).

59. Hartmann, A., Weiler, M. & Blume, T. The impact of landscape evolution on soil physics: Evolution of soil physical and hydraulic properties along two chronosequences of proglacial moraines. *Earth Syst. Sci. Data* **12**, 3189-3204 (2020).

60. Eichel. Vegetation succession and biogeomorphic interactions in glacier forelands. in *Geomorphology of proglacial systems* (ed. Heckmann, Tobias; Morch, D.) 327-349 (Springer, 2019).

61. Schoenecker, K. A., Zeigenfuss, L. C. & Augustine, D. J. Can grazing by elk and bison stimulate herbaceous plant productivity in semiarid ecosystems? *Ecosphere* **13**, (2022).

62. Olofsson, J., Moen, J. & Oksanen, L. Effects of herbivory on competition intensity in two arctic-alpine tundra communities with different productivity. *Oikos* **96**, 265-272 (2002).

63. Van Der Wal, R. & Brooker, R. W. Mosses mediate grazer impacts on grass abundance in arctic ecosystems. *Funct. Ecol.* **18**, 77-86 (2004).

64. Bardgett, R. D. & Wardle, D. A. Herbivore-mediated linkages between aboveground and belowground communities. *Ecology* **84**, 2258-2268 (2003).

65. Johnson, A. E., Molyneux, R. J. & Ralphs, M. H. Senecio: a dangerous plant for man and beast. *Rangelands Arch.* **11**, 261-262 (1989).

66. Cornelissen, J. H. C. A. et al. A handbook of protocols for standardised and easy measurement of plant functional traits worldwide.
335-380 (2003).

67. Xiong, D. et al. SPAD-based leaf nitrogen estimation is impacted by environmental factors and crop leaf characteristics. *Sci. Rep.* **5**, 1-12 (2015).

68. Díaz, S. et al. Plant trait responses to grazing - A global synthesis. *Glob. Chang. Biol.* **13**, 313-341 (2007).

69. Coley, P. D., Bryant, J. P. & Chapin, F. S. III. Resource availability and plant antiherbivore defense. *Science* (80-). **230**, 895-899 (1985).

70. Seibert, P. The Vegetation of the Settlement Area of the Callawaya People and the Ulla-Ulla Highlands in the Bolivian Andes.
*Mountain Research and Development.* **14**, 189-211 (1994).

71. Bardgett, R. D., & Wardle, D. A. Aboveground-belowground linkages: biotic interactions, ecosystem processes, and global change.
(Oxford University Press, 2010).

72. Janzen, D. H. Dispersal of Small Seeds by Big Herbivores: Foliage is the Fruit Author (s): Daniel H. Janzen Source: The American Naturalist, Vol. 123, No. 3 (Mar., 1984), pp. 338-353 Published by The University of Chicago Press for The American Society o. 123, 338-353 (2016).

73. Bodmer, R. E. Strategies of seed dispersal and seed predation in Amazonian ungulates. *Biotropica* **255-261** (1991).

74. Bryant, F. C. & Farfan, R. D. Dry season forage selection by alpaca [Lama pacos] in Southern Peru. *J. Range Manag.* **37**, 330 (1984).

75. Veblen, T. T., Mermoz, M., Martin, C. & Kitzberger, T. Ecological impacts of introduced animals in Nahuel Huapi national park. *Argentina. Conserv. Biol.* **6**, 71-83 (1992).

76. Bruno, M. C. & Hastorf, C. A. Gifts from the camelids: Archaeobotanical insights into camelid pastoralism through the study of dung.
*The Archaeology of Andean Pastoralism* (2016).

77. Hastorf, C. A. & Wright, M. F. Interpreting wild seeds from archaeological sites: a dung charring experiment from the Andes. *J. Ethnobiol.* **18**, 211 (1998).

78. Root-Bernstein, M. & Svenning, J. C. Prospects for rewilding with camelids. *J. Arid Environ.* **130**, 54-61 (2016).

Acknowledgements

This work was funded by the United States National Science Foundation under the Human-Environment and Geographical Sciences Program (Award No: 2105826) and National Geographic Society (Grant No. EC-58123R-19). This study was carried out under the scientific research authorization from the SERNANP (licenses No 008-2019-SERNANP-JEF and N°005-2022-SERNANP-JEF) authorizing the collection of plant material. We would like to express our deep gratitude to the Community *Llama 2000* (in particular to Calixto Huerta Trujillo, Alejandro Cruz Salvador, and Jorge Martel Alvarado) for their continuous involvement and contribution to the experiment. We also thank Luis Lizardo Oscanoa Gamarra for his help designing the experiment. We would like to acknowledge the human and material facilities provided by the National Park of Huascarán. We thank Giovana Patricia Vadillo Galvez from the Plant Physiology Laboratory (UNMSM) for the support in the identification of the seeds, and Pier Cisnero Jurado, Joshua Lopez Robles, Edu Jara Tarazona, Fernando Valverde Damián, and Inés Sachahuaman Balbin for their unconditional help in the field. Thanks to the Laboratory for Environmental Archaeology of UT Austin for the use of their equipment and to the multiple undergrad students of funding from the Soils and Geoarchaeology Laboratory of UT Austin for helping with the soil analysis. We thank Eugenio Arima for his statistical advice. Finally, the authors are also grateful to seven anonymous reviewers, whose comments on the manuscript greatly improved the final article. No experiments on live animals were carried out and we certify that the study is exempt from IACUC (Institutional Animal Care and Use Committee) review and approval.

Author contributions

A.Z., T.B., and F.A. conceived the ideas and designed the experiment; A.Z., S.R., and J.S. collected the data; A.Z., F.A., and T.B. analyzed the data; A.Z. led the writing of the manuscript. F.A. and T.B. reviewed, added to and edited the manuscript, and R.C. analyzed and provided the deglacierization data. All authors contributed critically to the drafts and gave final approval for publication.
**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-023-41458-x.

**Correspondence** and requests for materials should be addressed to A.Z.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023