Anthropogenic activities in the páramo trigger ecological shifts in Tropical Andean lakes

Abstract

A robust understanding of the impact of anthropogenic activities on high-altitude tropical aquatic ecosystems is key for the conservation and protection of the Tropical Andean biodiversity hot spot. We present the results of a multiproxy study of lake sediments from the high Andean páramo of El Cajas National Park, a UNESCO biosphere reserve in Ecuador. The main site, Laguna Pallcacocha, is well known for recording El Niño–driven clastic flood layers that are triggered by high-intensity rainfall anomalies from the eastern Pacific. The second site, Laguna El Ocho, does not contain clastic laminations, providing a control. The records show abrupt shifts in diatom assemblages ca. AD 1991 in both high-elevation Andean lakes accompanied by local changes in páramo composition that suggest a sudden nutrient enrichment of the environment. The diatom assemblages from Laguna Pallcacocha, in relation to the clastic input events, are remarkably stable and do not show evident El Niño signals at the analysed resolution. Based on comparison with the nonlaminated El Ocho record, we deduce the main source of this nutrient enrichment to be the construction of a heavily transited road that runs through the park, while climate warming played secondary role by amplifying its effects.

INTRODUCTION

Climate change and increased anthropogenic pressure over the past century have resulted in ecological regime shifts and species extinctions in many ecosystems worldwide (Ceballos et al., 2015; IPBES, 2019). Understanding these ecological shifts, especially in biodiversity hot spots, is of importance, as they often result in profound changes in ecosystem services that are costly, complex, or even impossible to reverse (Hughes et al., 2013). The Tropical Andean biodiversity hotspot (Myers et al., 2000), where recent climate warming may have forced high-elevation lakes towards new ecological and physical states (Michelutti et al., 2015), is a well-known example. These lakes are situated in the UNESCO protected region of the El Cajas National Park (CNP), Ecuador, and are considered an important drinking water source for the nearby city of Cuenca (Ministerio del Ambiente del Ecuador, 2018).

Climate and sedimentation in El Cajas, Ecuador

In the absence of long-term observational and monitoring data from the CNP, sediments from these lakes have been used as natural archives to gain insights into aquatic ecological shifts since the Industrial Revolution (Michelutti et al., 2015; Benito et al., 2022). Fossil diatom assemblages obtained from these lacustrine sediments have revealed abrupt replacements of species in the past ~50 yr. Based on the assumption that the lakes are pristine and isolated from direct human activities in their catchments, it was suggested that rising temperatures and longer periods of thermal stratification in the lakes are the main mechanisms behind elevated abundances of certain diatom taxa (Michelutti et al., 2015, 2016). If these algal communities are indeed strongly temperature driven, they could present a sensitive proxy for reconstructing rapid temperature shifts further back in time in the Tropical Andes, a region where local climate variability on interannual timescales is primarily driven by the El Niño Southern Oscillation (ENSO) (Kiefer and Karamperidou, 2019). For example, the lithostratigraphy of sediments from Laguna Pallcacocha, a small lake situated in the páramo of the CNP, is a key record for reconstructions of past ENSO dynamics (Moy et al., 2002; Hagemans et al., 2021; Mark et al., 2022). Current palaeoclimatological interpretations from the lake archive are based on depositional patterns and elemental composition of clastic layers in the lake record (Moy et al., 2002; Mark et al., 2022), as well as pollen assemblages registering concurrent landscape disturbance (Hagemans et al., 2022). Recent debate on the interpretation of the sedimentological proxy for ENSO reconstructions (Schneider et al., 2018) has led to renewed study of the exact sediment cascade and meteorological conditions that lead to the layered organic–clastic pattern in the lake sediment record. The resulting new chronological data from both subbasins of Laguna Pallcacocha and the alluvial fans and improved high-resolution meteorological data show a tight link between high-intensity rainfall anomalies, debris flow activity, and clastic sedimentation patterns in the lake record and confirm it as a valid ENSO proxy (Hagemans et al., 2021).

Anthropogenic-driven change

In addition to the need to understand climatic impacts, the recent ecological changes in the lake basin and the catchment are unknown, while knowledge of the ecological responses to potential human interference is important for development of other ecological and climatological proxies in the lake. ENSO-driven sedimentation patterns have potentially altered nutrient input to the site, and the associated rainfall anomalies likely impact the vegetation. A broad assessment of recent drivers and impacts is therefore important for understanding the system's response to natural and human-made disturbances, both for its significance in paleoclimatology and biodiversity research.

Despite their protected status, the lakes have been directly or indirectly subjected to increased anthropogenic pressure in the past century, which might obscure the observed climate signals. Indeed, the combination of widespread climatic and local anthropogenic disturbances has been shown to lead to dynamically stable lake diatom assemblages across lakes in the Ecuadorean Andes (Benito et al., 2020). Moreover, Giles et al. (2018) demonstrated a differential response to climate change between the shallow and deep lakes in the páramo, with the shallow ones showing only minor climate-related changes in diatom assemblage composition. Diatom-inferred regime shifts in the Ecuadorian lakes have been related to the combination of intensified agricultural use and human-induced climate change (Benito et al., 2022), as well as to nutrient increases due to alteration of the páramo landscape and increased cropland area (Giles et al., 2018). Intensified cattle farming in the páramo, for example, has been linked to changes in phytoplankton community composition in the lakes (Colen et al., 2018). Furthermore, intentional burning of the páramo and chemical pollution of the lakes from road transport are considered primary threats to the sensitive ecosystems of the CNP (Bush et al., 2017).

Such anthropogenic activities might overprint or even amplify climate signals in biological proxies obtained from lake sediments in the CNP park. In the western Amazonian lowlands, for example, a subannually resolved sedimentary record from Lake Sauce showed that the El Niño climate signal only became detectable in the record after people increased land-clearing activities, such as setting fires and growing crops, in the past ~2000 yr (Bush et al., 2017). Hence, the consideration of the signature of people on Tropical Andean ecosystems is essential for the adequate interpretation of paleoclimate proxies. Although the impacts of recent climate warming and the anthropogenic footprint on the ecology of the CNP lakes have been addressed in earlier studies (Michelutti et al., 2015, 2016; Labaj et al., 2017; Benito et al., 2019, 2022), we need additional insight into the potential effects of the various anthropogenic activities. As an example, the impacts of intentional burning of the páramo and recent road construction have not been assessed.

Multiproxy lake records

Here, we study aquatic community composition in sediment surface samples along an altitudinal transect in high-elevation Andean lakes in the CNP. By building upon our earlier studies (Hagemans et al., 2019, 2021), we hypothesize that shifts in the lake microfossil assemblages will reveal an ENSO-driven signal, as the precipitation anomalies that trigger clastic sedimentation (Kiefer and Karamperidou, 2019) might increase nutrient and/or silica input. Together with precipitation anomalies, temperature increases during El Niño events (Vuille et al., 2000) in turn might have effects on biota at high elevations. To test our hypothesis, we perform detailed multiproxy (fire, land use, and botanic) analyses of sediment records recovered from two contrasting lake types, the deep Laguna Pallcacocha and the shallow Laguna El Ocho and ask: (1) Is temperature the main driver of changes in aquatic community composition across various lakes at different altitudes in the CNP? (2) Is there an ENSO-driven signal in the diatom records of both lakes, Pallcacocha and El Ocho? (3) Will our sediment archives reveal a signal of nutrient increase in the studied lakes, related to the postindustrial human impact? (4) Will our multiproxy data reveal other anthropogenic impacts in the CNP?

To test our assumption and answer the research questions, we combine high-resolution diatom data, microcharcoal particles, and pollen analyses for a comprehensive reconstruction of the in-lake and catchment processes over the past century. We link ecological data from multiple proxies with various sensitivities (temperature, nutrients, catchment erosion, element deposition, etc.) to detailed records of local land use and conservation practices for the lakes and the vegetation of the CNP to detect pre-impact conditions and assess causes for (rates of) ecological change.

STUDY REGION

ECP (3100–4500 m above sea level [m asl]; Fig. 1a) is situated on the western Cordillera of the southern Ecuadorian Andes (2°46′46″S, 79°13′27″W), ~25 km west of Cuenca in the province of Azuay (Rodbell et al., 2002; Hansen et al., 2003). Bedrock in the CNP consists mainly of thick Quaternary silicic ignimbrite, rhyolite, and andesite of the Tarqui Formation (Rodbell et al., 2002). The vegetation between 3100 and 3500 m asl is an upper montane Andean forest transgressing into páramo, which occurs from 3500 to 4500 m asl. The vegetation in the CNP consists of ~500 species of vascular plants, of which 71 are endemic to Ecuador and 16 species are unique to the CNP (Ministerio del Ambiente del Ecuador, 2018). Average air temperature in the region is 6.5°C with no significant seasonality but a strong diurnal variation: from −1.8°C to 17.8°C on a daily average (Carrillo-Rojas et al., 2016). Precipitation ranges from 800 to 1500 mm per year and varies with the occurrence of warm El Niño events (Vuille et al., 2000; Moy et al., 2002; Kiefer and Karamperidou, 2019). The park has been managed by the local utility company ETAPA (Empresa Pública Municipal de Telecomunicaciones, Agua Potable y Alcantarillado y Saneamiento de Cuenca) since AD 2001 (Ministerio del Ambiente del Ecuador, 2018). Human-ignited fires in the montane forest and páramo occur regularly and form a major threat to the conservation of the natural environment in the CNP. Traditional livestock keeping consists of letting animals such as llamas, cattle, and horses graze in the páramo of the CNP. However, due to cattle theft in the mid-1990s farmers now keep their cattle closer to their farms, reducing the number of animals grazing in the páramo (Ministerio del Ambiente del Ecuador, 2018). Fish farming occurs in some of the lakes, for example, Laguna Toreadora (Ministerio del Ambiente del Ecuador, 2018). In AD 1991, a road was constructed in the northern part of the CNP (ETAPA, personal communication, 2019) that is heavily transited by trucks, cars, and buses (Fig. 1d).

Figure 1. Study region and sites. (a) Map showing the different vegetation types of the study region, outline of El Cajas National Park (CNP; grey), the major road (black), and sampling locations for the multiproxy analysis. (b) Photograph of Laguna El Ocho and its catchment. (c) Photograph of Laguna Pallcacocha and its catchment. In band c, a star indicates the coring location in the lake. (d) Lake catchments and their distances to the road. Coring locations are indicated with red stars.

METHODS

Sampling and site characteristics

Surface sediment, plant material, and water samples were collected with a UWITEC gravity corer and subsampled on-site from eight lakes in the CNP to assess changes in modern diatom community composition along a steep air temperature gradient of 7°C, that is, 3100–4200 m asl (Hagemans et al., 2019; Supplementary Table S1). Water samples were additionally used to determine the physical and chemical properties of the lakes, including temperature, pH, alkalinity, conductivity, and concentrations of major elements. Two small lakes situated in the CNP, namely Laguna Pallcacocha and Laguna El Ocho, were selected for obtaining paleoecological data. In August 2015, sediment cores of 50 cm length were retrieved from the deepest point of each lake with an UWITEC gravity corer. Laguna Pallcacocha is situated at an altitude of 4050 m asl (2°46′13″S, 79°13′59″W) (Fig. 1b, Supplementary Table S1). The lake has a surface area of 0.05 km², maximum water depth of ~8.3 m, and a pH of 7.54. The lake basin is partly filled in by a small river delta, which supplies abundant clastic material when activated by high-intensity rainfall events (Hagemans et al., 2021). Water retention time of the lake is 0.5 yr, with input from catchment runoff, direct precipitation, and minimal groundwater flow and output via river outflow and evaporation (Hagemans et al., 2021). In AD 2015 a bathymetry map was constructed using an echo sounder, and a sediment core was collected from the deepest point in the southern subbasin with an UWITEC gravity corer (PAL IV). The 46-cm-long core was split and photographed, and element ratios were measured with an X-ray florescence scanner at 1 mm resolution (see details in Hagemans et al., 2021). One half was sliced at 0.5 cm (0–16 cm) and 1 cm (16–46 cm) intervals for ²¹⁰Pb dating and diatom analysis. The other half was sliced at 0.2 cm (0–20 cm) and 1 cm (20–46 cm) intervals and used for pollen analysis. The unlaminated top 9.6 cm of the PAL IV core was analysed at very high resolution, as low clastic input likely results in lower sedimentation rates. Laguna El Ocho is a shallow lake located at an altitude of 3990 m asl (2°47′56″S, 79°12′57″W) (Fig. 1c). The lake has a surface area of 0.02 km², a maximum water depth of 4.2 m, and a pH of 7.42 (Supplementary Table S1). Retention time of the lake is estimated at 0.5 yr, based on an average annual precipitation and evapotranspiration of respectively 1000 and 700 mm/yr and a neglected input through groundwater flow. For the calculation of the lake water volume, an average depth of one half the maximum depth has been used. Catchment sizes are derived from the digital elevation model (based on SRTM data). A 50-cm-long gravity core was collected from the northern, deepest subbasin of the lake (OCH) and subsampled on site at 1 cm intervals for ²¹⁰Pb dating and diatom and pollen analysis.

Chronology

The age–depth models of Laguna Pallcacocha and Laguna El Ocho cores are based on ²¹⁰Pb isotopes analysed by gamma spectroscopy at Utrecht University. The cores were analysed at 0.5 cm to 1 cm (PAL IV) and 2 cm (OCH) intervals, respectively. The ²¹⁰Pb inventories were counted on a broad-spectrum germanium detector following standard procedures using freeze-dried samples. A chronology for both cores was established using a constant rate of supply (CRS) (Sanchez-Cabeza and Ruiz-Fernández, 2012), assuming that sedimentation rates may vary throughout the core.

The age–depth model of Laguna El Ocho was established on the basis of the ²¹⁰Pb dating and was interpolated with the cubic-b spline method in Tilia v. 1.7.16 (Grimm, 1991). Negative values in the observed ²¹⁰Pb profile of Laguna El Ocho (Fig. 2c), caused by near-background values that are within the measurement uncertainty of the detector, were corrected by adding the sum of the negative counts to the total ²¹⁰Pb (+35.7 Bq/kg) for calculation of the unsupported ²¹⁰Pb, which renders the base of the record ~20 yr younger, but mostly affects the levels with below-background signals (26–50 cm). The age model of Laguna El Ocho should therefore be interpreted with some caution. The effects of the minimum correction (+18.4; the mean of negative values) and maximum correction (+46; value leading to no negative supported or unsupported ²¹⁰Pb values) age scenarios on the stratigraphic signal are shown in Supplementary Figure S6.

Figure 2. (a) ²¹⁰Pb activity throughout the core from Laguna Pallcacocha. (b) Constructed age model for Laguna Pallcacocha showing the fluviolacustrine model and the cubic-b spline model. (c) ²¹⁰Pb activity throughout the core from Laguna El Ocho. (d) Constructed cubic-B spline age model for Laguna El Ocho.

The age–depth model of Laguna Pallcacocha consists of a linear interpolation for the unlaminated and consolidated top 9.5 cm and a model based on sedimentation rates for the laminated part of the sediment core, already reported in Hagemans et al. (2021). For the highly laminated record of Laguna Pallcacocha, it is necessary to account for short-term fluctuations in sedimentation rates that are the result of clastic sediment input during high-intensity rainfall events (Rodbell et al., 1999; Moy et al., 2002; Hagemans et al., 2021). We therefore applied an age–depth interpolation model developed to account for variable sedimentation rates in fluviolacustrine deposits (Minderhoud et al., 2016). The fluviolacustrine age–depth model (Fig. 2) corrects the CRS ²¹⁰Pb model for short-term sedimentation rate changes using a continuously sampled red colour intensity (RCI) record. The RCI is a measure of variation between rapid clastic sediment input from the Pallcacocha catchment and the background autochthonous organic deposition in the lake itself (Rodbell et al., 1999; Moy et al., 2002). The laminae show a characteristic pattern that can be correlated across the southern basin (Hagemans et al., 2021; Mark et al., 2022). The sedimentation rates, initially constrained by the high-resolution ²¹⁰Pb dates, are varied as a function of the RCI. We preferred this method over eliminating clastic layers (as presented in Mark et al. [2022] and Schneider et al. [2018]), as it explicitly incorporates changing sedimentation rates that occur frequently in Laguna Pallcacocha and is consistent with the approach followed by Rodbell et al. (1999) and Moy et al. (2002) for the Holocene records.

Water chemistry and diatom analyses

The pH values of the lakes were measured with an Orion BCG ELE 014 electrode. The alkalinity of the lakes was determined in the lab via titration with an acid of 0.010 mM to reach a buffer capacity. The conductivity was measured with an EC electrode, and the oxygen content with a galvanic dissolute oxygen electrode. Major element concentrations were determined by X-ray fluorescence spectroscopy with an ARL Perform'X 4200W high-performance WDXRF instrument at Utrecht University. Between 0.1 and 0.04 g of freeze-dried sediment was processed for diatom analyses at every 0.5 to 1 cm (PAL IV) and at every 2 cm (OCH) down to the base of the core. Additionally, a total of 28 surface sediment, plant and water samples were processed and analysed to assess modern diatom communities along the altitudinal lake transect. Samples for diatom analyses were processed by using cold H₂O₂ and 10% HCl to oxidise the organics and remove the carbonates (Renberg, 1990). To calculate the absolute diatom concentrations (DC; [n] valves × 10⁷/g dry sediment), known aliquots were used for slide preparation using the sedimentation tray method of Battarbee (1986). Concentration values were converted to influx values (DI; n/cm/yr) using the sedimentation rate and sediment density. Permanent diatom slides were prepared with Naphrax® as the mounting medium. Diatom slides were counted to approximately 350 diatom valves per slide. Counting was conducted with oil immersion at 1000× magnification using a Leica DM LB2 microscope. Diatom identification followed Hofmann et al. (2013), Krammer and Lange-Bertalot (1986, 1991a, 1991b, 1997, 2000), and Rumrich et al. (2000). Delineation of the diatom zones in both data sets (including only species at >2% relative abundance) is based on the constrained incremental sum of squares (CONISS) cluster analysis (Grimm, 1987). Multi-proxy stratigraphic diagrams were plotted with the C2 software (Juggins 2007).

Pollen and microcharcoal

Pollen samples were processed from the sediments of Laguna Pallcacocha at 0.2 to 1.0 cm intervals down to 20 cm and Laguna El Ocho at 2 cm resolution down to 50 cm. Volumetric sediment samples (approximately 0.5 ml PAL IV; approximately 0.5 g of freeze-dried sediment for OCH) were spiked with 4 mL of Lycopodium clavatum solution containing approximately 8339 (PAL IV down to 9.6 cm), 3864 (PAL IV 9.6–20cm), and 7942 (OCH) spores. Samples were treated with 10% KOH at 70°C and sieved over 200 μm mesh to remove humic acids, to remove coarse fragments, and to disaggregate the sediments, respectively. Samples were dewatered with 99% acetic acid and acetolysed for 10 min at 100°C in a 1:9 mixture of sulphuric acid and acetic anhydride (Faegri and Iversen, 1989). The sediment surface samples were floated over sodium polytungstate (density = 2.0) to remove heavy minerals. Residues were mounted in glycerol and analysed to a minimum of 300 pollen grains and spores with a Leica DM2500 light microscope at 400× magnification. Pollen types were identified to the lowest taxonomic level possible following descriptions by (Hooghiemstra, 1984) and by comparison with the reference collections at the University of Amsterdam and Utrecht University. Pollen percentages for each taxon were calculated relative to the total pollen number found in the specific sample (except aquatic taxa such as Isoëtes spp.). Percentage diagrams were plotted with C2 (Juggins, 2007).

Initial evaluation of the coarse fraction obtained from the palynological preparation yielded very few macroparticles. Given the limited amount of core material available, charcoal particle analysis was based on the microscopic fraction <120 μm. Microcharcoal particle numbers and surface area and L. clavatum occurrences were recorded using digital image analysis in ImageJ v. 1.52a (Schneider et al., 2012). The analyses were based on tile scans of microscopic slides at 200× magnification produced with a Leica DM6000 B automated stage light microscope. Estimates of total microcharcoal influx were based on L. clavatum added per volume of sediment.

Multivariate analysis

Variation in the data was further explored with ordination analyses to summarize changes in species composition. Diatom and pollen data were resampled to achieve the same resolution and create combined data sets. Percentage data were square-root transformed, rare species (<2%) were deleted, and for Laguna El Ocho we omitted Alnus acuminata from the analyses, due to extraordinarily high levels at ca. AD 1949, to minimize their influence on the outcomes. For both lakes, data are compositional and have a gradient length of <1 SD units, and the recommended method (Šmilauer and Lepš, 2014), a principal components analysis (PCA), was performed in C2 (Juggins, 2007).

RESULTS

Chronology

The ²¹⁰Pb age control points indicate that the homogenous top part of the sediment core (0–9.5 cm) from Laguna Pallcacocha covers ~22 yr, indicating that the first clear laminae at 9.5 cm sediment depth was deposited around AD 1993 (Fig. 2). The age model of Laguna Pallcacocha is corroborated by the charcoal record showing increased fire events in AD 2009–AD 2012, a decline in AD 2013 and AD 2014, and an increase in AD 2015, which is in accordance with fires reported in the CNP by the park managers (Ministerio del Ambiente del Ecuador, 2018). The sediments of Laguna El Ocho do not show any laminations, and the ²¹⁰Pb age control points indicate a lower sedimentation rate compared with Laguna Pallcacocha. The ²¹⁰Pb age control points show that the top 10 cm cover ~27 yr (Fig. 2). The overall Pallcacocha sedimentation rates we observe agree well with the recently published chronology of Mark et al. (2022), which is based on a core collected in AD 1999. The clastic laminae in the PAL IV record have maximally a 1 to 2 yr offset to the ENSO 1+2 index maxima (Hagemans et al., 2021). Similarly, the analysis by Mark et al. (2022) presents very high consistency between the amount of laminae and archive-based occurrences of El Niño. The sediments formed after approximately AD 1999 show more variability due to a large debris flow and rerouting of the main sedimentation pattern into the lake (Hagemans et al., 2021).

Modern aquatic conditions and diatom communities

The eight lakes sampled along the altitudinal transect have comparable limnological properties with circumneutral pH values (Fig. 3, Supplementary Table S1). Alkalinity and conductivity are lowest in Laguna Piñancocha (0.29 mM/L; 22.89 μS) and highest in Laguna Llaviucu (0.93 mM/L; 51.17 μS; Supplementary Table S1). Elemental phosphorus (P) concentrations were very low, below detection limits in most lakes, indicating ultra-oligotrophic conditions. The diatom communities in lakes Llaviucu, Taitachugo, Buríno Grande, Toreadora, and Pallcacocha are dominated by Discostella stelligera (Cleve & Grunow) Houk & Klee, followed by Tabellaria flocculosa (Roth) Kützing, Navicula spp. Nitzschia spp., and Achnanthes spp. In Laguna Piñancocha, the highest and coldest site, Aulacoseira sp.1 dominates the diatom community. The shallow lake El Ocho and the Inka bog are dominated by benthic species, such as Achnanthes spp., Fragilaria spp., Navicula spp., and Cocconeis placentula Ehrenberg (Fig. 3).

Figure 3. Diagram showing the dominant diatom taxa observed in the surface sediment, plant, and water samples collected from various lakes in Cajas National Park (CNP), Ecuador. Lake elevation (meters above sea level), mean air temperature (°C), and measured pH are displayed for each lake.

We identified 169 diatom taxa in the Laguna Pallcacocha sediment core. We divided the record into three zones based on the CONISS cluster analysis (Fig. 4). The first zone, PAL 1 (ca. AD 1960– AD 1993) is characterized by low DI (<100 cm/yr diatom valves). Although Aulacoseira sp.1 is dominant in the assemblage, the small epilimnetic D. stelligera increases its relative abundance towards the top of the zone (Fig. 3, Supplementary Fig. S1). The transition to zone PAL 2a (ca. AD 1993–AD 2006) is characterized by a sharp increase in the relative abundance of D. stelligera at the expense of decreased Aulacoseira sp.1 abundances. The benthic taxa, Achnanthidium spp., increase in abundance, reaching maximum percentage values ca. AD 2005. DI values vary in this zone, with maxima at ca. AD 2002 (1025 cm/yr diatom valves) and between ca. AD 2009–AD 2011 (1100–1260 cm/yr diatom valves). Despite a decrease in the relative abundances, D. stelligera remains dominant (~35–50%) throughout zone PAL 2b (ca. AD 2006–AD 2015).

Figure 4. Multiproxy diagram for Laguna Pallcacocha (core PAL IV) showing red colour intensity; Br/Ti ratio from X-ray florescence (From Hagemans et al., 2021); charcoal influx; diatom influx; the relative abundances of selected diatom, algae, and pollen taxa; total upland pollen influx; and PCA axis 1 and PCA axis 2 scores. Red dashed lines represent clastic input after high rainfall events. Full red line marks the period of road construction.

We identified 123 diatom taxa in total in the sediment core from Laguna El Ocho. The CONISS clustering delineated two distinct diatom assemblages (Fig. 5, Supplementary Fig. S2). No influx values were calculated for Ocho given the greater uncertainty of the age model. The lowermost zone, Ocho 1a (ca. AD 1870–AD 1965) is characterized by relatively high DC (~150–300 × 10⁷/g diatom valves) and dominated by benthic species, including the oligo-mesotrophic species Navicula wildii Lange-Bertalolot (Van Dam et al., 1994) and Encyonopsis subminuta Krammer & Reichardt. DC in Ocho 1b (ca. AD 1965–AD 1992) remain relatively similar to those in in zone Ocho 1a, ranging from ~150 to 250 × 10⁷/g diatom valves. The benthic, meso-eutrophic species Sellaphora pupula (Kützing) Mereschkovsky and Sellaphora laevissima (Kützing) Mann (Van Dam et al., 1994) and Neidium spp. slightly increase in abundance compared with OZ 1a and co-dominate the assemblage. Zone Ocho 2 (ca. AD 1992–AD 2015) marks an increase and the highest relative abundances of the small centric planktonic species D. stelligera (~10–22%). Nitzschia oberheimiana Rumrich & Lange-Bertalot and the meso- to eutrophic Nitzschia palea (Kützing) W.Smith and Nitzschia spp. (Van Dam et al., 1994) are abundant throughout this zone. The DC values are lowest in this part of the record, varying between ~30 and 100 × 10⁷/g diatom valves.

Figure 5. Multiproxy diagram for Laguna El Ocho (core OCH) showing the diatom and total upland pollen concentrations (with and without Alnus), relative abundances of selected diatom and pollen taxa, and PCA axis 1 and PCA axis 2 scores. Note that pollen percentage is based on trees, ferns, and herbs but without Alnus due to overrepresentation. Full red line marks the period of road construction.

Pollen and charcoal

In total, 65 pollen taxa were identified in the sediment samples from Laguna Pallcacocha and Laguna El Ocho. Following the zonation of the diatom assemblages of Laguna Pallcacocha, PAL 1 (ca. AD 1960–AD 1993) shows a strong increase in Poaceae (from 38% to 63%), while Apiaceae and Cyperaceae decline (Fig. 3, Supplementary Fig. S3). Montane forest taxa, such as A. acuminata and Hedyosmum spp. remain stable. Isoëtes novo-granadensis Fuchs, an aquatic submerged plant, increases strongly in abundance in this zone. Charcoal concentrations remain relatively stable, with a peak in AD 1978 (Fig. 4). In PAL 2a (ca. AD 1993–AD 2006) pollen percentages of all taxa remain relatively stable, although an increase in the abundance of Podocarpus spp. and Cyperaceae is observed at the end of this zone. Isoëtes novo-granadensis shows a sharp decline in abundance (from 63% to 13%) in this zone. Botryoccocus sp. and Pediastrum boryanum var. boryanum (Turpin) Meneghini, both planktonic algae, decline. PAL 2a (ca. AD 2006–AD 2015) is characterized by an increase in Poaceae and Podocarpus spp., while other taxa remain relatively stable. Isoëtes novo-granadensis remains stable between ca. AD 2004 and AD 2012 but increases again at ca. AD 2012–AD 2015 (Fig. 4). Charcoal influx values fluctuate, with relatively high values between AD 1963 and AD 1970 and maximum values around AD 1990 followed by lower values. Clastic sediment layers, identified from high RCI and low Br/Ti vales and indicative of high rainfall events (Hagemans et al., 2021), line up with troughs in the charcoal influx. Study of moss pollsters and surface lake sediments have shown that the pollen signal from the montane forest in CNP lakes positioned above the treeline represents sources of up to 40 km distance (Hagemans et al., 2019). Given the comparable particle sizes, the microcharcoal fraction likely represents a regional signal, with significant input from areas beyond the direct lake catchment area. A recent Holocene record from the CNP valley site Laguna Llaviucu (3150 m asl) within the montane forest shows significant macrocharcoal content. Areas in the valley were regularly burned (Nascimento et al., 2020) and were a potential source for the microcharcoal deposited in Pallcacocha. However, in Llaviucu, charcoal yield declined in the post-Columbian period, suggesting either sources changed or the overall burning intensity in the CNP is relatively low in this period. The Holocene Pallcacocha microcharcoal record also shows a recent decline, but only over the last century (Hagemans et al., 2022), confirming the Pallcacocha charcoal signal represents regional signal in the CNP.

The low RCI values after AD 1997 are related to a shift of the contributing stream away from the lake depot centre, reducing overall clastic input (Hagemans et al., 2021). Taxa such as Apiaceae, Polylepis sp., Hedyosmum spp., and Plantago spp. have positive loading on axis 1, and Poaceae, Podocarpus spp., and Ambrosia sp. have negative loadings on the PCA first axis. PCA 1 revealed clear changes in composition of the terrestrial pollen assemblages, with a marked shift around AD 1990 to negative values (Fig. 4). The total pollen influx values show relatively stable values of ~15,000/cm²/yr until AD 1990, followed by lower influx of ~5000/cm²/yr and, from AD 2000, a gradual increase to >25,000/cm²/yr. In contrast to the charcoal record, the pollen composition (PCA 1) and influx values do not show an evident relationship with the clastic sediment input events at this resolution.

Relative to the diatom-based zone boundaries of Laguna El Ocho, the pollen assemblages in zone Ocho 1a (ca. AD 1870–AD 1965) are characterized by a substantial increase in abundance of A. acuminata at ca. AD 1938 with a maximum at ca. AD 1949 and much higher values than in PAL IV, while other montane forest taxa, such as Hedyosmum sp. and Melastomataceae, show short phases of increased abundance (Fig. 5, Supplementary Fig. S4). Towards the end of zone Ocho 1a, Poaceae decline. In Ocho 1b (ca. AD 1965–AD 1992), A. acuminata abundance declines together with Cecropia spp., while Juglans neotropica Diels and Poaceae increase. During zone Ocho 2 (ca. AD 1992–AD 2015), Poaceae and J. neotropica increase to their highest relative abundances in the record, while aquatic taxa such as I. novo-granadensis and Pediastrum boryanum var. boryanum decrease. Alnus acuminata remains low in this zone. We observed consistent but low occurrences of Zea mays Linnaeus pollen throughout the record (Fig. 5).

Multivariate analysis

The PCA results clearly support the stratigraphic zone delineation with main shifts in principal component 1 and 2 values occurring at major zone boundaries (Figs. 3 and 4, Supplementary Fig. S5). For Lake Pallcacocha, a total of 55 diatom and pollen variables were included in the model. Planktonic diatoms, Isoëtes novo-granadensis, Pediastrum boryanum, and macrocharcoal show positive scores on PCA axis 1 (eigenvalue = 0.45), while Aulacoseira spp., Alnus acuminata, and montane forest taxa show negative scores. Discostella stelligera and Poaceae show negative PCA axis 1 and positive PCA axis 2 values (eigenvalue = 0.22). For Laguna El Ocho, 42 variables were included in the PCA model. The small planktonic diatom D. stelligera, benthic, meso-eutrophic Nitzschia species, and Poaceae show positive PCA axis 1 (eigenvalue = 0.39) and negative PCA axis 2 (eigenvalue = 0.17) scores. Isoëtes novo-granadensis, macrocharcoal, and N. wildii have negative PCA axis 1 and 2 scores.

DISCUSSION

The sediment surface samples of the lakes recorded remarkably similar diatom assemblages despite differences in local average annual air temperature (Fig. 3). For the surface samples, temperature alone does not explain the presence, abundance, or distribution of the most abundant diatom species, such as D. stelligera. In fact, D. stelligera is present and abundant in all deeper lakes despite differences in air temperature of up to 7°C. This suggests that diatom assemblages in CNP do not directly reflect a temperature signal but a more complex shift in lake ecosystem functionality that depends on multiple factors.

Overall, the fossil diatom composition in both lakes, Pallcacocha and El Ocho, primarily reflects the physical differences in water depth and lake watersheds. Due to its shallower water depth, El Ocho shows more short-term compositional variability within the benthic assemblage, while planktonic species dominate the diatom assemblages of the deeper Laguna Pallcacocha. Despite these differences, our multiproxy study revealed a similarly timed shift in diatom communities and pollen assemblages in the CNP over the past century, reflecting ecological changes in the lakes and their surroundings. The small epilimnetic species, D. stelligera, replaced Aulacoseira sp.1 in the planktonic communities of both lakes since ca. AD 1992 (Fig. 6). Concurrent with the increased abundance of D. stelligera, we also observe shifts in the local páramo vegetation of the catchment area ca. AD 1991 at both sites. The pollen record indicates clear compositional shifts at ca. AD 1991 that are mostly related to an increase in Poaceae pollen abundance at the expense of Apiaceae and Cyperaceae (Figs. 4–6). The regional pollen signal, consisting of the upper and lower montane forest taxa (Hagemans et al., 2019), is comparable between both sites, although the changes are more pronounced in Laguna El Ocho and it has higher input from upper montane forest pollen, initially Alnus and then Juglans, Myrica, and Melastomataceae. Earlier study of the pollen content of surface samples along the elevational gradient in the CNP showed similar results (Hagemans et al., 2019). In comparison to Laguna Pallcacocha, Laguna El Ocho is smaller and more exposed, and it is facing the forested valley around Laguna Llaviucu on the east side of the CNP, which is a rare area of relatively intact Andean forest (Fig. 1a). The only apparent forest change at Laguna Pallcacocha is a decline in Polylepis spp., a local treeline indicator (Hagemans et al., 2019) at the upper boundary of Zone Pal 1b at ca. AD 1993.

Figure 6. Diagram showing the timing of increase in the relative abundances of Discostella stelligera in sediment records from lakes Pallcacocha, El Ocho (this study), Toreadora, Chorreras, and Llaviucu (data from Michelutti et al. [2015]), alongside changes in Poaceae abundances and PCA axis 1 and 2. Data are compared with the ENSO 1+2 index (NOAA, 2020) and temperature (°C) and precipitation (mm) data obtained from the Cañar meteorological station, situated east of Cajas National Park (CNP) and maintained by the Ecuadorian Meteorological and Hydrological Service (INAMHI). The red line marks the period of road construction.

The shift at ca. AD 1992 is also clearly shown in the ordination analyses, where a change in the PCA axis 1 sample scores summarise the simultaneous increase in the abundance of meso-eutrophic diatoms and grassland pollen in both lakes (Figs. 4–6).

Similar diatom species replacements, as observed in our record, have been reported for other lakes within the CNP, specifically, Toreadora, Llaviucu, and Chorreras (Fig. 6), and related to climate warming (Michelutti et al., 2015). Rising temperatures and reduced wind speeds leading to longer periods of thermal stratification have been suggested as the main mechanisms behind the elevated D. stelligera abundances in the deeper lakes of the park (Michelutti et al., 2015). Since AD 1939, Tropical Andean near-surface temperatures exceeded those of the nineteenth century AD, tripling to 0.34°C per decade since the 1970s (Vuille et al., 2000). Also, the rise of small planktonic diatoms, such as D. stelligera and related centric species from Thallasiosiraceae, has been reported as a widespread phenomenon that occurred over the last two centuries in various lakes across the Northern Hemisphere (Smol et al., 2005; Saros et al., 2012). Although the phenomenon has been generally linked to climate warming, its triggering mechanisms have been widely debated. For example, while some experimental studies show that temperature affects diatom growth rates of D. stelligera (Williamson et al., 2010), others found no effect of temperature on growth rates or cell densities for the same species (Saros et al., 2014). The response of D. stelligera and similar Cyclotella taxa to light, nutrients, and temperature is complex and often controlled by interactive effects between multiple drivers (Malik and Saros, 2016). In situ bottle experiments with natural phytoplankton assemblages from Arctic lakes have shown that temperature alone was insufficient as a driver to explain observed Cyclotella population dynamics (Malik and Saros, 2016). In fact, temperature mediated the effects on water-column stability and, consequently, light and nutrient availability (Malik and Saros, 2016).

Average daily air temperature (maximum and minimum) measured in 19 stations across Ecuador has increased by 0.25°C per decade between AD 1966 and AD 2011 as a result of climate warming (Morán-Tejeda et al., 2016). Because our results show that D. stelligera is present and abundant in lakes with an air temperature difference of up to 7°C, we consider it unlikely that a temperature increase of 0.25°C per decade since the 1960s can solely explain the sudden shifts in diatom communities. Therefore, temperature alone does not explain the presence, abundance, or distribution of D. stelligera and other diatom species in the high Andean lakes of the CNP. Furthermore, a decline of Pediastrum boryanum, a green algae taxon previously shown to rapidly increase in response to climate warming in remote Arctic lakes (Woelders et al., 2018), occurred in the 1990s. Although we acknowledge that climate warming may have provided a background change allowing rapid shifts in diatom assemblages to occur in our studied lakes, we submit that it is unlikely that steadily increasing atmospheric temperatures alone caused the rather sudden shift observed in our records.

The coherent shifts in the diatom communities in the lakes and species composition in the páramo point to an elevated nutrient input in the ecosystems of the CNP. As páramo lakes are oligotrophic and the present-day P concentrations are extremely low (Supplementary Table S1), even a small increase in nutrient loading can have an impact on lake ecology (Colen et al., 2018). The coherent increase in abundance of D. stelligera and Nitzschia taxa over the past three decades suggests elevated nutrient concentrations in the investigated lakes. Discostella stelligera abundance responds strongly to nutrient enrichment, as indicated by in situ experiments in Arctic lakes (Saros et al., 2014). Additionally, high abundances of Nitzschia species, and particularly N. palea, is associated with higher trophic levels (Licursi et al., 2016). Also, increased turbidity in lakes and streams can promote the growth of motile diatoms, such as Nitzschia species (Dickman et al., 2005), and planktonic species that have the ability to resist high turbidity (Rimet and Bouchez, 2012).

Indeed, the shifts in aquatic and terrestrial ecology coincide with the construction of a heavily transited road in AD 1991 that runs through the CNP (ETAPA, personal communication, 2019). Road infrastructure and traffic are an important source of environmental pollutants such as dust, nutrients, heavy metals, and other chemical substances that enter tropical aquatic ecosystems during intense rainfall events and runoff (Trombulak and Frissell, 2000; Laurance et al., 2009). Records of chemical pollution from nearby lakes Fondococha and Llaviucu show the sharp increase after AD 1950 of atmospheric (dust) sources with heavy metal pollution related to industry and traffic (Schneider et al., 2021). As such, nutrient enrichment and chemical pollution can significantly affect species composition in tropical water bodies (Trombulak and Frissell, 2000). Dust is an important source of P to freshwater systems (Brahney et al., 2014), while fuel combustion from vehicles on the road elevates nitrogen (N) deposition in terrestrial and aquatic ecosystems (Coffin, 2007). In addition to road construction, fertilization by cattle grazing in the páramo (Colen et al., 2018) and burning of vegetation (Ministerio del Ambiente del Ecuador, 2018) have likely contributed to enhanced nutrient loads via catchment runoff and atmospheric deposition around the time of the main transition. In the lakes studied here, we observe peaks in diatom concentrations in periods of increased fire occurrences in the CNP area, although this trend is only visible since ca. AD 1991. Additionally, natural and anthropogenic biomass burning in Amazonia can interfere with nutrient cycling in downwind ecosystems such as the Ecuadorian Andes. Elevated deposition of H⁺, NO₃, N, and Mn have been reported in the Ecuadorian Andes during the biomass burning season in Amazonia (Boy et al., 2008). In our research lakes, where the dominant wind direction is from Amazonia (Hagemans et al., 2021 and references therein), acidification and eutrophication from Amazonian fires could also have influenced lake ecology.

Furthermore, industrial development of Andean economies could affect the sustainability of Andean water resources. During the twentieth century, rapid industrialization has led to an increase in emissions of heavy metals and polycyclic aromatic compounds (PAC). Increased deposition of heavy metals (Cooke and Abbott, 2008; Rodbell et al., 2014) and PACs (Bandowe et al., 2018) has already been reported in several previously pristine Andean lakes, potentially negatively impacting lake-water quality.

Our results show that elevated nutrient levels affect not only the aquatic ecosystems, but also the vegetation composition in the páramo. While no regional studies on páramo grassland composition in relation to nutrient deposition are available, in general, eutrophication of grassland communities is known to result in elevated plant productivity (Hautier et al., 2009). Nutrient enrichment increases competition for light between plant species (Hautier et al., 2009), favouring faster-growing or taller species such as the nitrophilous Poaceae. This competition for light can result in substantial plant biodiversity loss following eutrophication in grasslands (Hautier et al., 2009). Indeed, we observe clear shifts in species composition in the páramo, especially a clear expansion of Poaceae abundance at both sites, after the construction of the road at ca. AD 1991 (Fig. 6). The substantial increase of A. acuminata is possibly linked to the occurrences of two major earthquakes in the region in AD 1938 and AD 1949 (NOAA, 2020). Alnus acuminata rapidly colonizes areas recently subjected to mass movements that are triggered by intense rainfall (Clark et al., 2016) and occasional earthquakes (Crausbay and Martin, 2016). The A. acuminata phase is successively replaced by an increase in Melastomataceae and Juglans neotropica, which are secondary forest taxa that confirm regrowth after short-lived landscape disturbance, possibly aided by the creation of the park.

The phases of enhanced clastic input seen in the RCI proxy at Laguna Pallcacocha, which reflect ENSO-driven high-intensity rainfall events (Hagemans et al., 2021; Mark et al. 2022) as well as air temperature anomalies (Vuille et al., 2000), show no clear relation to the percentages of pollen and diatoms. The highest variability is present in the Isoëtes novo-granadensis abundances, but no systematic relation with RCI maxima is evident. An exception is the microcharcoal record. Periods with high charcoal influx seem to correspond to low RCI, suggesting an inverse relation between rainfall intensity and fire intensity likely mediated by ENSO. High sedimentation rates within the clastic layers are corrected by the age model according to the colour intensity (see “Methods”), largely removing the dilution effect of high sedimentation rates in the microcharcoal influx signal, although some bias may remain. A potential relation between clastic input and diatom productivity was not demonstrated, suggesting silica is not a limiting factor and/or the clastic material does not contribute significant nutrients to cause productivity increase. Contrary to our hypothesis, the present set of proxies seem relatively insensitive to high-frequency variability, and the signals are dominated by the assemblage turnover around AD 1992 (Fig. 6). In some intervals, the variable sample resolution of the pollen and diatom data limits rigorous testing of this relationship, and subsequent detrending of the data might reveal more high-frequency variability. Also, a definite answer on this requires further high-resolution analyses, especially in the pre-1995 period.

CONCLUSION

Aquatic diatom communities in the CNP have very comparable composition and, despite differences in lake type, show strong consistent increases in epilimnetic taxa, especially D. stelligera in recent decades. We conclude that atmospheric nutrient enrichment triggered abrupt species replacements and community shifts in the lakes and vegetation in the CNP, while climate warming has provided a background change allowing the rapid shifts to occur but is not a major driver.

Temperature proxies, based on calibration of current aquatic algae compositions, are therefore likely biased by increased nutrient deposition. We deduce the main source of abrupt nutrient enrichment to be the construction in AD 1991 of the heavily transited road that crosses the northern part of the CNP. Other potential sources of nutrient addition to the lakes and vegetation are cattle abundance, páramo burning, and Amazonian forest fires, but exact time-series data on their intensity are not available. These results demonstrate that even in the remote setting of a high-altitude natural reserve, human impact can lead to changes in aquatic conditions and vegetation shifts that are not necessarily related to climate forcing. Grass cover in the páramo shows similar bias from recent nutrient increase, although on decadal timescales the overall vegetation composition seems more influenced by episodes of land clearance and physical disturbance. Comparison of these results with other nutrient-poor aquatic ecosystem changes suggests that distance to intensively inhabited areas is no guarantee for low environmental impact.