Abstract
High-mountain lakes constitute key archives of past environmental dynamics and climate change. This study aims to reconstruct Holocene vegetation dynamics, fire regimes and hydroclimatic history of a high-altitude pond in central Mexico using a multiproxy approach. The Nahualac pond is located on the NW flank of the Iztaccíhuatl volcano at 3900 m asl. This multiproxy investigation integrates the analysis of pollen, non-pollen palynomorphs, diatoms, cladocerans, ostracods, chironomids, testate amebae, charcoal, and sediment geochemistry. Three stages of environmental change were identified. The first stage (8.4–6.6 ka) was characterized by a wetland environment with alpine grasslands, and low fire activity. In the second stage (6.6–4.7 ka) planktonic diatoms increased, and there was an upward displacement of the montane forests dominated by Pinus, Quercus, and Alnus and enhanced fire activity. A transitional period (4.7–4.2 ka) showed dry-cold environments, similar to those of the first stage. The third stage (last 4.2 ka) was characterized by wetter conditions and fluctuating water levels. A first interval (4.2–3.5 ka) recorded a shallow pond with planktonic diatoms, cladocerans, and ostracods indicating cold winters, the establishment of Pinus forest, and increased fire activity. From 3.5 ka to the present, a pond ecosystem with complex trophic interactions developed. Three drought events were identified at ~3.2–2.7, ~2.0–1.9, and ~1.2–1.1 ka. Overall, Nahualac provides a robust high-altitude reference showing how interactions among temperature, seasonality, and moisture balance governed ecosystem dynamics and fire regimes in the central Mexican highlands over the last 8.4 ka.
Introduction
High-mountain environments are particularly sensitive to climate change (Diaz et al., 2014). They contain key archives of past environmental dynamics and climate change, including glacial landforms, lake and bog sediments. In tropical mountains these archives often record marked changes from glacial and periglacial in the late Pleistocene to near tropical forest conditions in the Meghalayan. Because of their elevation, they record the global temperature rise of the Holocene but are also sensitive to cooling short-lived episodes like the Little Ice Age, as well as moisture fluctuations related to latitudinal shifts of the Intertropical Convection Zone (Heine, 2024).
In central Mexico the Northgrippian and Meghalayan are largely missing on lake records from low and mid elevation areas due to strong anthropic disturbance (Caballero et al., 2025). In contrast, high mountains are much less disturbed and may provide valuable insights into climate and environmental change during the last few millennia.
During the cold periods of the late Quaternary, extensive mountain glaciers developed in central Mexico at elevations above 3,800 m asl (Vázquez-Selem and Heine, 2011). The most complete glacial history of Mexico, and particularly of the Basin of Mexico (BM), is based on studies carried out on the Iztaccíhuatl volcano, where geomorphological evidence has documented variations in glacier extent over the last 200,000 years (Vázquez-Selem and Heine, 2011; White, 1962). The advance of these glaciers, which descended several hundreds of meters, reshaped the geomorphology of the mountains, creating new landscapes.
Iztaccíhuatl (5215 m asl) is located at ca. 19° N in the southeastern part of the Basin of Mexico, a volcano-tectonic depression in the east-central sector of the Trans-Mexican Volcanic Belt. The BM is a high-altitude endorheic basin (2250 m asl) that developed an extensive lacustrine system, with Lake Chalco to the south and Lake Texcoco at its center. Iztaccíhuatl volcano formed between 0.9 and 0.08 Ma by the accumulation of andesitic and dacitic lavas emitted from several vents aligned along a ~20 km north–south axis (Nixon, 1989). It forms part of a volcanic chain that also includes Popocatépetl (5450 m asl), an active stratovolcano that, repeatedly deposited tephra throughout the late Pleistocene and Holocene, mantling the entire mountain range and large parts of the adjacent Basins of Mexico and Puebla (Miehlich, 1991; Siebe et al., 1996).
Due to the complex topography of central Mexico, the upper mountain regions contain sedimentary archives such as peatlands, alpine lakes, and depressions some of which were carved by former glaciers. Paleoecological research conducted at some of these high-altitude sites (Lozano-García and Vázquez-Selem, 2005; Rodríguez-Pérez et al., 2026; Ruiz-Córdova et al., 2019) suggests that they function as sensitive recorders of climate change, with the advantage of being less affected by human activity than lakes located at lower elevations (Caballero et al., 2025). Moreover, Mexico’s high-mountain lakes represent unique marginal conditions within the country’s environmental envelope and within the Trans-Mexican Volcanic Belt (Avendaño et al., 2023), characterized by high elevation, low temperatures, strong seasonality, reduced atmospheric pressure, and generally oligotrophic, low-conductivity waters. As a result, they offer exceptional opportunities to study ecological processes related to diversity, endemism, and environmental dynamics.
An example is the pair of high-mountain lakes on the Nevado de Toluca volcano (Lakes El Sol and La Luna), which occupy marginal ecological conditions (low temperatures and low pH values) in this environmental space and host endemic species not reported from any other sites (Avendaño et al., 2023; Ruiz-Córdova et al., 2019). These lakes date to approximately 6900–6100 cal year BP (6.9–6.1 ka hereafter), formed under periglacial conditions, and evolved into fully established lacustrine systems around 6.0 ka, with sedimentation continuing uninterrupted to the present. Another example is the Valle Agua el Marrano site, a depression formed during the Milpulco-1 glacial advance (12.5–10 ka), at 3860 m asl on the north slope of Iztaccíhuat and filled with colluvial deposits. Palynological analysis of this sequence has documented Holocene vegetation history, revealing shifts in community composition and recording the upward trend of the treeline migration (Lozano-García and Vázquez-Selem, 2005). These records reveal pronounced altitudinal and climatic changes, as well as vegetation dynamics, between ~11 and 2 ka.
Paleoenvironmental studies of lacustrine sequences from the lowlands (ca. 2,200 m asl) of the BM, such as those from Lake Chalco, have documented vegetation changes, lake-level fluctuations, volcanic activity, and multiproxy evidence spanning the Pleistocene and Holocene (Avendaño et al., 2023; Caballero and Guerrero, 1998; Lozano-García et al., 2022; Lozano-García and Ortega-Guerrero, 1998; Lozano-García et al., 1993; Ortega-Guerrero et al., 2020). However, the human impact is so extensive that the record of the last 6000 years is mostly disturbed (Caballero et al., 2025).
In this study, we present the paleoecological record from the Nahualac pond, at 3900 m asl on the northwestern slope of Iztaccíhuatl. We integrate multiple proxies—including pollen and non-pollen palynomorphs, macroscopic charcoal and its morphometric attributes, diatoms, testate amebae, cladocerans, ostracods, chironomids, and sediment geochemistry—to reconstruct Holocene environmental dynamics. We focus on (i) vegetation change, including treeline fluctuations; (ii) reconstructing fire activity and fuel types through macroscopic charcoal morphometry; and (iii) on the process of pond formation and limnological development. Our specific objectives are to: (1) assess local and regional environmental changes during the last 8.4 ka, particularly through comparison with lower-elevation records such as Lake Chalco; (2) evaluate site responses to climatic variability throughout the Holocene; (3) determine climatic trends and ecological responses during the Northgrippian (~6 ka), a period often described as relatively warm in regional records, identifying climatic trends toward humid or dry conditions in this sector of the BM and comparing with modern ecological trends associated with present-day warming, in a qualitative analog framework; and (4) identify drought events during the Late-Holocene and their potential ecological impacts.
Study site
The Nahualac pond, located at 3900 m asl (19°10′ N; 98°40.7′ W), lies in a depression on the upper surface of a lava flow that originated near the summit area of Iztaccíhuatl (Pecho Peak; Figure 1). The lava flow stands out from ca. 4,100 m asl down to its frontal wall at 3750 m asl It belongs to the younger andesites and dacites described and dated to the Middle Pleistocene (0.58–0.35 Ma; Nixon, 1989). The lava flow was shaped by glacial erosion during at least the last glacial cycle. It was entirely covered by glaciers during the Hueyatlaco–1 and −2 advances of Iztaccíhuatl volcano between 21–14.5 ka (Vázquez-Selem and Heine, 2011; White, 1962). Glaciers sculpted the irregular surface of the lava flow, producing a roche moutonnée morphology with distinctive glacial polish and striations. Glacial erosion likely overdeepened some pre-existing depressions on the lava flow surface, one of which hosts the Nahualac pond.
(a) Mexico (left) and the central-eastern portion of the trans-Mexican volcanic belt (right). NT: Nevado de Toluca; IZ: Iztaccíhuatl; PO: Popocatépetl; MA: Malinche; CP: Cofre de Perote; ch: lake Chalco; red inset: area represented in B, (b)the Sierra Nevada volcanic range formed by Iztaccíhuatl and Popocatépetl volcanoes; red star: Nahualac pond; red inset: area represented in C, (c) northwest side of Iztaccíhuatl (oblique view from west to east); red star: Nahualac pond; yellow line: lower limit of terminal Pleistocene (Milpulco-1) glacier; orange line: lower limit of early Holocene (Milpulco-2) glacier; white line: lower limit of Little Ice Age (Ayoloco) glacier.
A cosmogenic 36Cl exposure age from a glacial polish on the northern edge of the pond yields ages from 14.1 ± 1.0 to 17.0 ± 1.1 ka, depending on the scaling method used for age calculation (CRONUS web calculator, Marrero et al., 2016). These ages represent the timing of deglaciation of the site and provide a maximum age for the onset of sedimentation within the pond depression, once the surface became ice-free. The exposure age is consistent with the deglaciation following the Hueyatlaco-2 advance of Iztaccíhuatl (Vázquez-Selem and Lachniet, 2017). Evidence for three subsequent glacial advances has been identified at higher elevations than the Nahualac pond: Milpulco-1 (13–10.5 ka), Milpulco-2 (9–8 ka), and Ayoloco (Little Ice Age). The Nahualac lava flow was also mantled by tephras from major Plinian eruptions of Popocatépetl volcano, located 18 km, to the south. A fine pumice layer from the Ochre Pumice 5.9 ka; Arana-Salinas et al., 2010) was unequivocally identified in the pond sediments, while ash from the Pink Pumice sequence (ca. 1.1 ka BP; Siebe et al., 1996) mantles the lava flow and the surrounding areas.
At present, the Nahualac pond is a shallow waterbody (0.8–1 m depth) with a maximum area of ca. 1600 m2 that forms during the rainy season (May – October) and dries out completely in the dry season; the length of time it stays inundated depends on how wet the year has been (Hernández-Bautista and Junco-Sánchez, 2023). Its maximum area and depth are constrained by a low sill on its northwestern side. The pond is surrounded by Pinus hartwegii forest, which covers the area from 3700 to 4020 m asl, with tussock grasses forming part of the understory (Figure 2). The modern timberline (4020 ± 50 m asl) is constrained by cold and dry conditions (Lauer, 1978). From the timberline to ~4400 m asl, the dominant alpine communities consist of tussock grasses such as Muhlenbergia quadridentata, Calamagrostis orizabae, and C. tolucensis and Festuca toluscencis along with other species including Eryngium proteiflorum, Cirsium nivale, Lupinus montanus, Arenaria spp., and Juniperus monticola (Calderón de Rzedowski and Rzedowski, 2001; Steinmann et al., 2021). Abies religiosa forests develop between 2700 and 3600 m asl. At lower elevations, between 2300 and 3000 m asl, mixed Pinus–Quercus forests are dominant, with patches of cloud forests restricted to humid ravines between 2500 and 2800 m asl.
Nahualac pond.
The earliest known reference to the Nahualac pond dates back to the year 1747, on an old map with a land distribution record (Hernández-Bautista, 2020). In 1880, the French explorer Désiré Charnay reported on the pond, mentioning the existence of a rectangular monument rose in its center whose foundations are still extant; also, he mentioned smaller monuments like altars rounding the pond (Charnay, 1885). Nahualac was used as a ceremonial and symbolic site by the pre-Hispanic inhabitants of the Basin of Mexico. The architectonical elements reported by Charnay inside the pond were built and visited during pre-Hispanic times, between AD 700 (Epiclassic) and AD 1250 (early Postclassic) with a more emphasis in ritual use interval between AD 900 (final Epiclassic) and AD 950 (beginning early Postclassic) based on recent archeological studies. During these pre-Hispanic periods or later, the watercourse from nearby springs was diverted north of the Nahualac pond, and a drain was dug to the south in order to control the water level of the pond around the architectural elements (temple; Hernández-Bautista, 2020). This archeological context is particularly relevant for the interpretation of the Late-Holocene record, as human modification of the hydrological system may have influenced local water levels, vegetation structure, and fire regimes, and should therefore be considered when evaluating potential anthropogenic signals in the multiproxy dataset.
Materials and methods
Sedimentary sequence
A 95 cm-long sedimentary sequence (Nahualac 18-II; Figure 2) was obtained from the southeastern area of the pond using an Eijkelkamp percussion corer to recover an undisturbed core. Samples were taken continuously at contiguous 2 cm intervals. For radiocarbon dating, two samples were selected (Table 1). The age–depth model was constructed using the program Bacon v.3.3.1 (Blaauw and Christen, 2011), which calibrates radiocarbon dates with IntCal20 (Reimer et al., 2020). In addition, the known age of the Ochre pumice from Popocatépetl volcano (Arana-Salinas et al., 2010), identified between 66 and 69 cm, was used as an independent chronological marker (Table 1).
Radiocarbon date from Nahualac-18-II sequence.
Source: Arana-Salinas et al. (2010).
Geochemical analysis
Elemental concentrations were measured from 1 cm3 of sediment dried at 40°C for 24 h. Thirty-seven dried samples were ground with an agate mortar and analyzed with a portable Thermo-Scientific Niton XL3t XRF analyzer. Portable XRF data are semi-quantitative and therefore interpreted in relative rather than absolute terms. Total carbon and total nitrogen contents were determined on 32 dried samples using a Thermo-Scientific Flash 2000 CN analyzer.
Elemental ratios obtained from XRF analyses were used to infer changes in sediment sources, hydrological conditions, and nutrient dynamics. Titanium (Ti) is a conservative lithogenic element associated with terrigenous input, whereas calcium (Ca) is commonly related to carbonate precipitation or biogenic production; therefore, the Ti/Ca ratio reflects variations between detrital and authigenic inputs in lacustrine sediments (Boyle, 2001; Kylander et al., 2011). The Zr/Fe ratio was used as an indicator of sediment transport and hydrological variability, as zirconium is typically associated with coarse heavy minerals transported during higher-energy conditions (Davies et al., 2015). Phosphorus (P) concentrations were interpreted as a proxy for nutrient availability and aquatic productivity, given their close relationship with organic matter inputs and trophic state changes in lake systems (Boyle, 2001).
Biological analysis
Samples for biological analysis were selected according to the main changes observed in the pollen stratigraphy. Due to differences in material availability, the biological proxies (pollen, charcoal, testate amebae, diatoms, ostracods, cladocerans and chironomids) varied in their analytical resolution.
Pollen: A total of 0.5 g of wet sediment were processed (n = 26 samples). Before chemical treatment, two Lycopodium clavatum marker tablets (Batch 1031) were added. Samples were treated with 10% HCl and 40% HF for 24 h, followed by heavy-liquid flotation using sodium polytungstate (density 2.0). The residue was mounted in glycerine-gelatin for microscopic observation. A minimum of 200 terrestrial pollen grains were counted, excluding Pinus, aquatic taxa, and non-pollen palynomorphs (NPPs). Taxonomic identification was based on specialized literature and the reference collection of the Laboratory of Paleoclimatology, Paleoecology and Climate Change of the Institute of Geology, UNAM. Pollen data are expressed as percentages, concentration and influx and we used Tilia software (Grimm, 1991) for these calculations. Observations were made using a ZEISS AxioStar microscope at 400× and 1000× magnification.
Charcoal: Charcoal particles (>100 μm in size for n = 23 samples) were extracted from 1.0 cm3 sediment subsamples following Clark and Royall (1996). Each sample was deflocculated using hot 10% KOH and 10% Na4P2O7 to disaggregate clay and organic matter. Subsequently, the samples were sieved to selectively collect charcoal particles retained on the 100 μm mesh. The recovered material was photographed under an Olympus SZ2-IL-ST stereomicroscope at 1.6x. These images were then processed using ImageJ, and the total area of the particles was quantified. Based on the quantified area, charcoal concentration (mm2 cm−3) and charcoal accumulation rates (mm2 cm−2 year−1) were calculated. To infer changes in dominant fuel types, the length and width of all particles retained on the 100 μm mesh were measured using ImageJ (Schneider et al., 2012). Length-to-width (L/W) ratios were calculated for each particle and used to classify them into three fuel-type categories: L/W <2.5 (woody fuels), 2.51–3.49 (mixed), and >3.5 (non-woody), following practices commonly applied in lacustrine records and acknowledging recent critical evaluations of this proxy (Vachula et al., 2021). Relative intensity fire activity was inferred using Z-scores standardization of CHAR across the full sequence.
Diatoms: For diatom analysis (n = 23 samples), 0.2 g of dry sediment were treated with 10% HCl, 30% H2O2, and HNO3 at 80°C to complete oxidation. Cleaned material was mounted on permanent slides using Naphrax® as a mounting medium. At least 400 diatom valves per sample were counted under a ZEISS Axio Lab.A1 light microscope at 1000× magnification. Identification followed specialized literature (Fallu et al., 2000; Krammer et al., 1991; Krammer and Lange-Bertalot, 1986, 1988, 1991; Spaulding et al., 2021). Total diatom concentration was estimated using the aliquot settling method (Schrader, 1974) and expressed as valves per gram of dry sediment (v/gds). The identified species were classified according to their ecological habitats into three groups: planktonic–tychoplanktonic, benthic, and aerophilous. The complete list of species is provided in Supplemental Table S1 in the Supplemental material.
Invertebrates: For cladoceran, ostracod, and chironomid remains, 1 cm3 of wet sediment was sieved through 200, 100, and 63 µm meshes (n = 23 samples). Larger fractions were examined under a ZEISS Stemi 508 stereomicroscope, separating, counting, and identifying each group of bioindicators. Ostracod valves were isolated, counted, and identified under a compound microscope at 400× magnification. Taxonomic identification followed Meisch (2000) and Pérez et al. (2010). Ostracod abundance was expressed as the number of valves/cm3. Cladoceran and chironomid remains were processed jointly: 1 cm3 was treated with 10% KOH to deflocculate organic matter, rinsed with distilled water, and sorted in a Bogorov counting chamber under a stereomicroscope (400×). Remains were picked with a fine brush and stored in 1.5 ml Eppendorf tubes, separating cladoceran remains (head shields, carapaces, ephippia) and chironomid head capsules. Cladoceran remains were mounted on semi-permanent slides with glycerine, and chironomid capsules on permanent slides with Hydro-Matrix® mounting medium. Carapaces were the most abundant remains and were used to estimate individual abundance. Taxonomic identification was performed with a ZEISS Axiostar Plus microscope at 400× and 1000× magnifications, consulting specialized literature for cladocerans (Wojewódka et al., 2020a, 2020b) and chironomids (Andersen et al., 2013; Brooks et al., 2007; da Silva et al., 2018; Epler, 2001; Hamerlík et al., 2022).
Testate amebae: For testate amebae analysis (n = 21 samples), 1 cm3 was sieved through a 53 µm screen, and the upper fraction was completely examined. The analyses were performed under an AmScope LED Trinocular stereoscopic microscope up to 180× magnification. All tests found were separated with a very fine brush (20/0), counted, and identified following Kumar and Dalby (1998), which includes the use of ecophenotypes to describe morphologies associated with particular environmental conditions. The best-preserved specimens were photographed with scanning electron microscopy on a JEOL NeoScope JCM-600 microscope, low and high vacuum. To synthesize the testate amebae data, and given their generally low concentrations but high diversity, the taxa were grouped into two genera: Centropyxis (C. aculeata aculeata, C. aculeata discoides, C. constricta constricta) and Difflugia (D. urceolata urceolata, D. globosa, D. oblonga oblonga, D. glans; Supplemental Figure 2 in Supplemental material shows the individual species plots).
Statistical analysis: A Multiple Factor Analysis (MFA) was performed to integrate and simplify the multiproxy data (Bécue-Bertaut and Pagès, 2008; Escofier and Pagès, 1994). Diatoms and algae were first classified according to their ecological preferences. These data, together with pollen, cladocerans, ostracods, chironomids, testate amebae, XRF-derived elements, total organic carbon (TOC), total Nitrogen (TN), phosphorus, and charcoal accumulation rates (CHAR), were included in the MFA. Quantitative variables were standardized according to the procedure recommended by the FactoMineR package for quantitative data. Specifically, all variables were standardized to unit variance, ensuring comparability among datasets with different units and scales (Lê et al., 2008). Analyses were performed with the FactoMineR package (v.2.11; Lê et al., 2008) in R (v.4.3.3; R Core Team, 2009).
Results
Stratigraphy and chronology
The pond infill is colluvial, characterized by silty-loam textures with sandy intercalations and interbedded tephras (ash and pumice falls). Part of the colluvial material is rich in organic matter, as inferred from its very dark coloration, while another portion is predominantly light brown sediment. Between 66 and 69 cm, an ochre pumice-rich tephra was identified, corresponding to that described by Arana-Salinas et al. (2010) from Popocatépetl volcano (Figure 3).
Nahualac-18-II age model and lithostratigraphy.
The age–depth model indicates that the base of the Nahualac-18-II sequence is 8.4 ka, suggesting that sediment accumulation at the pond began in the early Holocene (Figure 3). The age–depth model is constrained by two radiocarbon dates and one tephra layer, which, although providing key chronological control points, results in relatively low temporal resolution for an 8.4 ka record. This limitation introduces uncertainties in the precise timing of environmental changes, particularly at centennial scales. Therefore, the ages assigned to the different phases should be regarded as approximate and interpreted within the confidence limits of the model. Consequently, short-term events, such as Meghalayan dry intervals, should be interpreted with caution, whereas the main trends discussed here are supported by consistent multiproxy signals and stratigraphic coherence throughout the sequence.
Multiple factor analysis
Due to the diversity of the proxies analyzed, stratigraphic zonation was derived from MFA (Figure 4) integrating all proxies. This approach jointly incorporates non-biological (Figure 5) and biological (Figures 6 and 7) proxies, accounting for their ecological particularities. The first two axes of the quantitative variables in the MFA explain 68% of the variance (Dim. 1:34%; Dim. 2:23%; Figure 4a). The axis scores are presented in Supplemental Material Figure 3, illustrating the evolution and establishment of the different environmental stages.
Multiple factor analysis (MFA): (a) circle of correlation of the database including geochemical, total nitrogen, total organic carbon, pollen, NPPs, charcoal accumulation rates, diatoms, and invertebrate assemblage (cladoceran, ostracod and chironomid) and (b) factor map showing sample ages, zonation, and four distinct environmental stages.
Geochemical parameters: Ti/Ca, Zr/Fe ratios, phosphorus (P), total organic carbon (TOC) in black and total nitrogen (TN) percentages in orange.
Summary diagram of Nahualac-18-II sequence with percentage and influx pollen and charcoal data.
Diagram of Nahualac-18-II sequence illustrating subaquatic pollen, non-pollen palynomorphes (NPPs), diatoms and invertebrates from the Nahualac sequence.
Dimension 1 separates the wetland environment from the aquatic one. On the negative side, variables such as aerophilous diatoms, Isoetes, Carex, ferns, and alpine grassland occur, indicating flooded soils and/or wetlands conditions, these correspond to the oldest samples (Figure 4b). In contrast, the positive side of Dimension 1 includes variables associated with higher moisture levels and the presence of a well-established water body—such as Zr/Fe (flooding indicator), planktonic–tychoplanktonic and benthic diatoms, testate amebae, invertebrates, and microalgae—together with the youngest samples. Ecologically, this axis reflects a transition from marginal wetland conditions to a more stable and deeper pond system, driven primarily by changes in hydrological balance and water availability through time.
Dimension 2 captures variability within the aquatic ecosystem once the pond was established: negative values represent periods of lower water levels and greater hydrological fluctuations, dominated by benthic and aerophilous taxa and associated with nutrient and organic-matter enrichment; positive values, by contrast, suggest a deeper water column in which the trophic chain—from phytoplankton to zooplankton—is more complex. This axis therefore represents changes in water depth, trophic structure, and ecosystem stability, allowing the identification of phases of environmental instability versus more persistent pond conditions.
Together, the MFA results highlight the strong coherence among biological and geochemical proxies, supporting the interpretation of coupled changes in hydrology, ecosystem development, and vegetation dynamics at Nahualac, and allowing the differentiation of distinct environments divided into: Zone 1 (94–76 cm; 8.4–6.6 ka), Zone 2 (76–57 cm; 6.6–4.7 ka), and Zone 3 (57–2 cm; 4.7–0.15 ka) with two subzones: 3a (57–47 cm; 4.7–3.5 ka), and 3b (47–2 cm; 3.5–0.15 ka).
Geochemical parameters
Total organic carbon (TOC) values ranged from 1.8-7.6%, and total nitrogen (TN) from 0.08% to 0.45% in the 95 cm sequence (Figure 5). In zone 1, the lowest TOC (1.8%) and TN (0.08%) values occurred at 91 cm, followed by an increase to 4% and 0.22% at 79 cm, after which both decreased again. In zone 2, TOC ranged between 3% and 3.5%, and TN between 0.13% and 0.18%, both decreasing toward 56 cm. During zone 3, TOC and TN showed a general increasing trend, both decreased at 10 cm, but reached their highest values toward the top of the sequence (Figure 5).
Elemental ratios (Ti/Ca, Zr/Fe) and phosphorus (P) exhibited similar trends, with low values during zone 1 and 2, reaching minima at the end of zone 2, followed by a marked increase with fluctuations during pond 1. In the zone 3b, Ti/Ca and P decreased, whereas Zr/Fe showed a slight increase before declining again (Figure 5).
Biological proxies
Pollen, non-pollen palynomorphs (NPPs) and charcoal
A total of 40 pollen taxa and 9 NPPs were identified across the sequence, dominated by Pinus, Quercus, and Alnus from temperate forest elements; Poaceae among herbs; and Gymnodinium, Debarya, Spirogyra, and Zygnema amid the NPPs. The charcoal record highlights marked variability in CHAR values and fire intensity (Z-scores), allowing three main zones to be distinguished (Figure 6).
During the zone 1, pollen influx was very low, averaging 37,977 grains cm-2 year-1. Dominant arboreal taxa included Quercus, Alnus, and Pinus, with Poaceae, Artemisia, Polypodium, and Isoetes also present. NPP concentrations were very low, with Gymnodinium, Spirogyra, and Zygnema recorded. CHAR values remain at moderate levels (~425 mm2 cm-2 year-1), and no samples show positive Z-scores, indicating a low recurrence of locally detectable fire events. During this interval, woody particles account for ~75%–80% of the total charcoal assemblage, while herbaceous fragments average ~14%, with local peaks of ~21% around 7.0 ka and ~16% near 8.3 ka (Figure 6).
During the zone 2, pollen influx increased to an average of 81,800 grains cm-2 year-1. Pinus dominated, with Abies, Juniperus, Quercus, and Alnus also present. Poaceae values were lower than in zone 1. NPP concentrations remained low, but increased relative to zone 1, including Gymnodinium, Cosmarium, Debyra, Mougeotia, Spirogyra, and Zygnema. At the end of the zone (58–56; 4.8–4.3%ka), pollen influx decreased to values similar to zone 1, with a decline in Pinus and an increase in Poaceae. There was a slight increase in Carex percentages, along with near absent NPPs. Charcoal data for this interval show intermediate CHAR values between the low-fire conditions of zone 1 and the intense fire activity of zone 3a, suggesting a transitional phase in which fire became more recurrent but still dominated by woody fuels (Figure 6).
In the terrestrial pollen record of zone 3, no marked differences are observed between the two subzones. Pollen influx rose further, reached its maximum (280,200 grains cm-2 year-1; mean 166,880), dominated by Pinus, Abies, Juniperus, Quercus, Alnus, Arceouthobium and Poaceae. Mesic forest taxa and herbs appeared, although at low values. NPP concentrations increased markedly, dominated by Gymnodinium, Debarya, and Zygnema. In the last 12 cm (0.76 ka) Poaceae increased and the NPPs concentrations declined sharply, particularly Gymnodinium, Debarya, Zygnema, and Spirogyra (Figure 3). This interval exhibits the highest fire activity of the entire record: mean CHAR values rise to ~1483 mm2 cm-2 years-1, and approximately 67% of samples display positive Z-scores, reflecting an increase in the frequency of combustion episodes. In the early part of this interval (4.7–3 ka), CHAR values are intermediate (~855 mm2 cm-2years-1), with ~50% of the samples showing Z >0, suggesting a gradual transition into this period of enhanced fire activity. Throughout zone 3, woody charcoal remains the dominant component (73%–82%), while herbaceous fractions are less abundant. During the last 1.3 ka, CHAR values decrease again (~548 mm2 cm-2 year-1), and no samples exhibit positive Z-scores. In the uppermost part of the sequence, (~0.36 ka), the proportion of herbaceous particles shows a modest increase (15%), although woody fragments continue to dominate the assemblage (Figure 6).
Correlation analyses reveal a positive but weak relationship between Z-scores and the woody fraction (r = 0.21; p = 0.34; n = 23), and a negative relationship with the herbaceous fraction (r = 0.27; p = 0.21). Comparable results were obtained using CHAR values. Although not statistically significant, these correlations suggest a general trend of greater fire intensity in association with woody biomass. Detailed correlation coefficients (r, o and n) are presented in Supplemental Table S2.
Diatoms, testate amebae and invertebrates (cladocerans, ostracods and chironomids)
Diatom abundance was very low during the zone 1, ranging from 5.7 to 18.9 × 106 valves g-1 dry sediment (v g-1 ds). The assemblage was dominated by aerophilous species such as Gomphonema acidoclinatum, Stauroneis obtusa, and Pinnularia borealis var. subislandica. Toward the end of this interval, Chamaepinnularia sp.—a species reported exclusively from Nevado de Toluca, Mexico—appeared. At the beginning of this period, planktonic–tychoplanktonic species such as Aulacoseira nivaloides and Fragilaria capucina occurred in very low abundances, below 5% (Figure 7).
Diatom abundance remained low in zone 2, ranging from 7.8 to 19.9 × 106 v g-1 ds, although showing a slight increasing trend. Aerophilous species continued to dominate during the early and middle parts of this interval, while toward its end, planktonic–tychoplanktonic taxa increased in abundance, including A. nivaloides, Cavinula pseudoscutiformis, F. capucina, F. capucina var. gracile, and Punctastriata mimetica.
The invertebrate and testate amebae records begin in zone 3a (4.9 ka). During this period, the highest abundance of cladoceran ephippia (33 resting eggs of Cladocera cm-3) in the entire record was obtained. The cladoceran taxon Alona was the most abundant (198 individuals cm-3), followed by the chironomid Paratanytarsus (18 head capsules cm-3). Ostracods were present in low abundance. Low abundances of the testate ameba Difflugia urceolata “urceolata” were recorded toward the end of this zone.
At the beginning of zone 3a, planktonic–tychoplanktonic species dominated, represented by A. nivaloides, Cavinula pseudoscutiformis, F. capucina, and F. capucina var. gracile. The benthic taxon Psammothidium helveticum was also present and showed a progressive increase in abundance. Also, the highest abundance of the ostracod Heterocypris cf. incongruens (23 carapaces cm-3) was recorded, as well as subfossil carapaces of the two cladocera genera (652 individuals cm-3). At the onset of this zone, Pleuroxus was the dominant cladoceran, but overall, Alona was the most abundant taxon during the early part of zone 3a. Additionally, this section shows the highest testate amebae abundances in the entire record (17 individuals cm−3), dominated by Centropyxis aculeata “aculeata”, C. constricta “aerophila”, and Difflugia urceolata “urceolata” (Figure 2 in Supplemental material).
Zone 3b begins with a marked decrease in diatom abundance, dropping from 133.5 × 106 to 12.5 × 106 valves g-1 dry sediment. Planktonic–tychoplanktonic species declined to less than 3%, while benthic and aerophilous taxa, including S. obtusa, Caloneis aerophila, P. borealis var. subislandica, Pinnularia shimanskii, and P. helveticum, became dominant. Also, the abundance of ostracods (three carapaces cm-3), cladocerans (50 individuals cm-3), and chironomids (four head capsules cm-3) decreased drastically. At 3.1 ka, ephippia reappeared in the record (5 cm-3). Subsequently, between 2.7 and 2.1 ka, benthic invertebrates such as ostracods or chironomids were absent; only cladocerans were recorded, with abundances below 70 carapaces cm-3. As observed in the other proxies, testate amebae remained scarce; in this interval, only Difflugia urceolata “urceolata” was present.
Planktonic–tychoplanktonic taxa increased again (17%–43%) between 2.1 and 1.3 ka, accompanied by a decrease in P. helveticum, while S. obtusa and P. borealis remained present. The three groups of invertebrates reappeared, with the highest abundance of chironomids (more than 30 head capsules cm-3) and the first occurrence of the genus Heterotrissocladius in the record. Around 1.3 ka, the second-highest occurrence of ephippia (22 cm-3) was recorded. Testate amebae abundance remained low, and only taxa of the genus Difflugia were observed.
Toward the top of the record (21.5–5.5 cm, 1.3–0.3 ka), total diatom abundance decreased markedly from 50.8 to 16.1 v × 106 v g-1 ds. Planktonic–tychoplanktonic species declined once more, and aerophilous taxa such as C. aerophila, P. borealis, P. shimanskii, and S. obtusa again became dominant. Regarding the invertebrate groups, Alona increases upward, reflected in the rise of total cladocerans from 227 individuals per cm-3 (1.3 ka) to 498 individuals per cm-3 (0.3 ka), while Pleuroxus remains low throughout. Ephippia occur only at 1.3 ka (22 ephippia per cm-3). Paratanytarsus increases toward the top, matching the rise in chironomids from 17 to 24 head capsules cm-3, whereas Heterotrissocladius is restricted to the lower level. H. cf. incongruens appears at 1.3 ka (eight valves cm-3) and reoccurs at 0.3 ka (four valves cm-3). In this section, testate amebae abundance begins at low values and increases toward 0.3 ka (13 ind. cm−3), accompanied by a rise in diversity marked by the presence of both Centropyxis and Difflugia.
Discussion
The integrated MFA framework shows the strong correspondence between biological and geochemical proxies, indicating that hydrological variability, ecosystem development, and vegetation dynamics at Nahualac evolved in a coupled manner. This coherence provides a robust basis for interpreting the environmental trajectory of the site and supports the identification of four main stages. These include an initial wetland phase (Zone 1: 8.4–6.6 ka), followed by a transitional interval (Zone 2: 6.6–4.7 ka), and a subsequent forested stage (Zone 3: 4.7–0.15 ka). Within this latter stage, two pond substages can be distinguished (Pond 1: 4.7–3.5 ka; Pond 2: 3.5–0.15 ka), corresponding to Zones 3a and 3b, respectively. In the following discussion, these intervals are addressed according to their environmental significance rather than their zonation. It is important to consider that differences in the temporal and spatial resolution of the proxies may influence their comparability. Consequently, proxies with higher sampling resolution, such as pollen, capture a greater number of observations and potentially finer-scale variability. This is particularly relevant for the transitional phase identified in the MFA, which is represented by only two samples in the integrated analysis, but may encompass a more complex sequence of changes when examined through higher-resolution proxies.
Vegetation dynamics: Forest expansion and fires regimes
Plant communities respond to environmental and climatic changes, and palynological records provide valuable information on past vegetation, allowing reconstructions of variations associated with different types of natural and/or anthropogenic disturbances (Chevalier et al., 2020; Stebich et al., 2015). At high-elevation sites, pollen data make it possible to analyze, among other evidence, vegetation responses to climate change, particularly to warming during the Holocene. For the central Mexican mountains, previous studies have also documented changes in vegetation composition during the Holocene (Lozano-García and Vázquez-Selem, 2005; Rodríguez-Pérez et al., 2026; Ruiz-Córdova et al., 2019).
A distinctive feature of the plant communities in these mountains is the dominance of temperate forests composed mainly of Pinus and Quercus species, accompanied by frequent elements such as Alnus and Abies. These taxa are characterized by high pollen productivity and efficient dispersal, which result in their strong representation in both modern and fossil pollen rain (Lozano-García et al., 2014; Lozano-García and Xelhuantzi-López, 1997; Rodríguez-Pérez et al., 2025). Given the high dispersal capacity of these taxa, their presence at high elevations may partly reflect extra-local inputs, which should be considered when interpreting vegetation changes at upper elevations. Studies of modern pollen rain across different vegetation types provide valuable information on the relationship between pollen signals and plant communities in these montane regions. For example, modern pollen rain along an altitudinal transect on Iztaccíhuatl, from 2650 to 4024 m asl, provides a conceptual framework for interpreting the Holocene pollen record of Nahualac (Rodríguez-Pérez et al., 2025).
The Nahualac pollen record, spanning the last 8.4 ka, reveals important fluctuations in pollen assemblages, suggesting vegetation shifts on the slopes of the volcano in response to temperature and moisture fluctuations. During the wetland stage (8.4–6.6 ka), Pinus percentages were relatively low (<60%), while Alnus and Quercus showed intermediate values. This signal is consistent with modern pollen rain on Iztaccíhuatl and Popocatépetl, where these taxa, characterized as long-distance transported, tend to be overrepresented at elevations above 3,900 m asl (Alnus and Quercus are absent nowadays above 3,500 m asl). The pollen composition suggests the presence of alpine grasslands dominated by Poaceae, along with Carex and Isoetes, which were grouped by the MFA, indicating wet conditions, whereas Artemisia represented open meadows and rocky slopes (Lu et al., 2022) or disturbance conditions (Calderón de Rzedowski and Rzedowski, 2001). Polypodium shows high values. In modern alpine communities of Iztaccíhuatl two genera of the Polypodiaceae family are present; however, no Polypodium has been reported (Hernández-Cárdenas et al., 2019); whereas at the alpine area of Cofre de Perote (a mountain 160 km to the east), Polypodium calirhiza is found (Mickel and Smith, 2004). During this period, the pollen influx was low, and the assemblages indicated that Nahualac had swampy conditions and the area around the pond was an alpine grassland. This stage may be considered an environmental prelude to the subsequent forest expansion.
An abrupt change occurred in the transitional stage, beginning at 6.5 ka, when pollen assemblages showed an increase in Pinus values, both in influx and percentages, continuing until 5.2 ka. Simultaneously, Quercus and Alnus decreased, while elements of humid forests such as Carpinus, Fraxinus, Ulmus, and Podocarpus became more prominent. The pollen spectra suggest an altitudinal upward migration of Pinus forests, with percentages reaching 78% at 6.5 ka. Modern pollen rain indicated the highest values (60%–85%) in the zone where pine forests currently occurred. A similar signal has been observed at Valle Agua el Marrano, on the northern slope of the volcano at 3870 m asl, where Pinus expansion was recorded at 6.1 ka (43%), reaching 92% by 6.0 ka (Lozano-García and Vázquez-Selem, 2005). In contrast, at Cofre de Perote (3717 m asl), percentages above 60% are reported earlier, at 7.6 ka (Rodríguez-Pérez et al., 2025). Discrepancies in the timing of upward migration likely reflect local site conditions and/or differences in age models.
As part of this upward migration process, Abies increased to 15% at 6.1 ka. In modern pollen rain, such high values were only reported below 3350 m asl (23%–2%), within the Abies forest, while above that elevation values rarely exceeded 5% (Rodríguez-Pérez et al., 2025). At present, Abies form dense forests between 2800 and 3500 m asl. The increase at Nahualac by 6.1 ka suggests an upslope expansion of this community, probably related to warming. In the Sierra La Negra, another high mountain in central Mexico, Steinmann et al. (2021) documented Abies religiosa nowadays at 4271 m asl, supporting the hypothesis that forest expansion and altitudinal extension reflected rising temperatures during the Holocene. Modeling of future scenarios of suitable climatic habitat changes for Abies religiosa indicates that this taxon can reach up to 4000 m asl in the highest elevations of central Mexico by 2030 (Gómez-Pineda et al., 2020).
Continuing with the hypothesis of upward migration during the transitional stage, the presence of Picea at 5.5 ka is noteworthy. This taxon was reported for the early Holocene (10–9.3 ka) at Cofre de Perote at 3717 m asl (Aguirre-Navarro, 2013; Rodríguez-Pérez et al., 2026), at Malinche at 3100 m asl in the terminal Pleistocene and early Holocene (Heine and Ohngemach, 1976), as well as in several mid-elevation sedimentary sequences (~2200 m asl), such as Chalco, Texcoco, and the Puebla–Tlaxcala basin (Rodríguez-Pérez et al., 2026). It is plausible that some Picea populations migrated upslope, where the high-mountain climate provided refugia, as in Cofre de Perote, and later in Iztaccíhuatl, where sites like Nahualac may have functioned as a microrefugium sensu (Rull, 2009).
A short-lived abrupt change occurred between 4.8 and 4.2 ka at the end of Zone 2 and the beginning of Zone 3, when the pollen spectrum resembled Zone 1, with low pollen influx, a decline in Pinus, and an increase in other temperate forest taxa. However, pollen assemblages recovered after 4.2 ka, during Zone 3. Pinus percentages stabilized between 60% and 80%, also Arceuthobium (a pine hemiparasite) showed an increase, and mesic forest elements along with herbs increased in diversity, suggesting slope stabilization. Nevertheless, a fluctuation in the pollen assemblages was recorded at 0.75 ka, where Pinus declined and Poaceae increased; also, Carex showed higher values suggesting that this change could be related to a dry episode.
Nahualac lies near the lower limit of the current ecotone between Pinus hartwegii forest and alpine grassland. This ecotone is diffuse, characterized by declining tree cover, and ranges from ~3920 to 4111 m asl on Iztaccíhuatl (Beaman, 1962). The ecotone is sensitive to climatic changes, as evidenced by recent warming: the treeline has shifted upward by an average of ~32 m between 1995 and 2017 (22 years), with progress documented predominantly in areas with gentle slopes and NW orientation (Guzmán-Vázquez and Galicia Sarmiento, 2025). Although altitudinal treeline boundaries are inherently difficult to define, particularly due to the influence of long-distance pollen transport, modern pollen data indicate that Pinus values below ~60% occur in alpine grassland areas. Based on this modern analog, the Nahualac record suggests phases of upward expansion of pine forests and probable rise of the upper forest limit on Iztaccíhuatl between ~6.5 to 5.5 ka at ~3900 m asl, followed by a re-expansion at ca. 4 ka with fluctuations. However, these interpretations should be considered as reflecting regional vegetation dynamics rather than precise local timberline positions, thereafter.
In this context, fire activity, as recorded by charcoal, provides additional insight into the ecological consequences of these vegetation changes. The charcoal signal integrates closely with the palynological data and with the proxy associations shown by the MFA analysis. In the intermediate part of the record, increases in Pinus and Quercus, together with declines in Poaceae, coincide with higher CHAR values and the dominance of woody charcoal particles (73%–82%), indicating a tight coupling between forest expansion, fuel availability, and fire activity. This supports the interpretation that woody biomass was the dominant fuel during the most intense burning phases. Given that the strongest charcoal peaks correspond to intervals of higher forest cover rather than increases in herbaceous material, the Nahualac record is consistent with fire regimes primarily governed by the distribution and availability of montane woody fuels. The increase in forest cover during this interval likely enhanced fuel availability, making montane forests more susceptible to fire during severe dry seasons. This relationship also suggests that climate variability indirectly modulated fire regimes through its control on vegetation structure and fuel continuity. Similar feedback between vegetation, temperature, and fire has been reported in tropical Andean records, where fire plays an active role in shaping the vegetation dynamics (Di Pasquale et al., 2008). These relationships should be interpreted cautiously, as charcoal morphometry is intended to provide qualitative information on dominant fuel types within a multiproxy paleoecological framework rather than to establish direct linear correlations with vegetation proxies.
However, the interpretation of the charcoal record must also consider the effects of charcoal taphonomic decay on the interpretation of fire signals. Processes such as oxidation, dissolution, and bioturbation can reduce both the quantity and size of preserved particles, particularly in older sediment sections, leading to an underestimation of early fire frequency. In tropical lacustrine sediments, it has been estimated that a non-negligible fraction of the original charcoal signal may be lost through diagenetic degradation over millennial timescales, due to oxidation, fragmentation, and microbial alteration (Bird, 2013; Feurdean, 2021; Huang et al., 2019; Scott and Damblon, 2010). In Nahualac, the relatively weak charcoal signal observed in the wetland stage could therefore reflect a combination of taphonomic loss and genuinely low fire recurrence associated with wetter conditions and limited ignition windows. Explicitly acknowledging both factors helps avoid over-attributing early low CHAR values solely to diagenesis and maintains consistency with the multiproxy evidence of cool, moist early Holocene environments.
Comparison with other high-elevation tropical sites places Nahualac within a broader regional context. In the northern Andes, records from Guandera (3,650 m asl; 0.6° N) and Quimsacocha (3,780 m asl; 2.9° S) show similar dynamics, with climate-driven fires dominating the mid-Holocene and an increasing human influence during the last two millennia (Di Pasquale et al., 2008; Jantz and Behling, 2012). In Chirripó, Costa Rica (3,820 m asl; 9.5° N), the record shows fire peaks between ~AD 560 and 720 and AD 980–1230, linked to seasonal droughts (Wu and Porinchu, 2020). In Mexico, the record from La Luna (Nevado de Toluca, 4200 m asl; 19.1° N) documents similar climatic and hydrological fluctuations, with reduced fire frequency during late cooling phases of the Little Ice Age (Caballero et al., 2020). By contrast, the Aljojuca maar (2360 m asl; 19.1° N) in the east-central Mexican highlands shows a Late-Holocene fire signal linked primarily to anthropogenic disturbance (Bhattacharya et al., 2015; Bhattacharya and Byrne, 2016). Together, these comparisons indicate that tropical mountain ecosystems respond coherently to moisture and temperature forcing, which modulate the effective ignition window and the dominant type of fuel available to burn, highlighting the sensitivity of high-elevation ecosystems to coupled climate–vegetation–fire interactions.
Aquatic ecosystem establishment and colonization
The aquatic, geochemical, and terrestrial proxies show broadly coherent patterns through time, allowing the reconstruction of coupled changes in hydrology, productivity, vegetation, and fire regime, although some differences among proxies likely reflect their distinct ecological sensitivities and temporal responses. Beyond short-term hydrological variability, the Nahualac record also suggests that longer-term catchment processes may have influenced pond dynamics. Changes in vegetation cover and fire regimes likely affected nutrient delivery, organic matter input, and sediment composition, thereby exerting a delayed control on aquatic communities. (Battarbee et al., 2005; Cole et al., 2007).
The initial stages of the Nahualac wetland (8.4–7.5 ka) are characterized by low Ti/Ca and Zr/Fe ratios the former indicating low detrital input and the latter reflecting flooding events (Davies et al., 2015), together with low percentages of TOC, a proxy for the abundance of organic matter derived from terrestrial vegetation and aquatic organisms. It also shows low TN, which reflects the nutrient status of the water body (Opitz et al., 2012; Wang et al., 2022), as well as low P values. Collectively, these proxies indicate that the initial stages of this pond had a low productivity.
During this time (8.4–5.5 ka) the record also shows the presence of microspores of Isoetes and Carex, both wetland plants, along with a suite of aerophilous diatoms including Gomphonema acidoclinatum, Pinnularia borealis var. subislandica, Caloneis aerophila, and Chamaepinnularia sp. (Figure 7), which typically develop in flooded soils and/or wetlands (Caballero et al., 2020). The occurrence of Chamaepinnularia sp. is restricted to this interval. This species was also recorded in the higher altitude Lake La Luna (Nevado de Toluca, 4200 m asl) during a limited time window around 5 ka. This aerophilous taxa could have been favored by colder-than-present conditions, particularly cooler winters and ice-cover conditions at times of higher seasonality (Caballero et al., 2020). Its presence at Nahualac during an earlier time interval reflects the altitudinal displacement of climatic conditions, forcing this taxon to migrate to the higher altitudes, where colder climates persisted, as warmer conditions were established during the Holocene. The dinoflagellate Gymnodinium, an indicator of cold waters (0°C–6°C) in Lake La Luna, appears only sporadically, as do chlorophyte algae such as Spirogyra and Zygnema. These genera also occur in the high-elevation lakes La Luna and El Sol on the Nevado de Toluca (4200 m asl; Cuna et al., 2022), and have also been reported from modern surface samples of soils, bogs, and lakes in the páramo and superpáramo zones between 3600 and 4150 m asl in the northern Venezuelan Andes (Montoya et al., 2010). The combined evidence from aquatic proxies points to the existence of a cold wetland with oligotrophic, circumneutral to slightly acidic conditions during the first ~3,400 years of the record (Zones 1 and 2). Together with the pollen evidence for alpine grassland and low arboreal cover, these data suggest that both the terrestrial and aquatic systems responded coherently to cool and moist early Northgrippian conditions. During the 8.4–5.5 ka period, the multiproxy data indicate a transitional stage with fluctuating environmental conditions between 5.5 and 4.5 ka. Around 4.9 ka, a rapid flooding event is evidenced by the occurrence of facultative planktonic diatoms such as Aulacoseira nivaloides and Cavinula pseudoscutiformis. This assemblage is very similar to the one present in modern Lake El Sol, suggesting the development of conditions similar to this lake, that is, an oligothrophic, low-conductivity, near-neutral (pH >6), cool permanent water body (Cuna et al., 2015; Ramírez-Nava et al., 2022; Zgrundo et al., 2017). At this time, cladocerans (Alona, Pleuroxus), ostracods (H. cf. incongruens), and chironomids (Paratanytarsus, Heterotrissocladius) appear for the first time in the Nahualac sequence, supporting the interpretation of a permanent, oligotrophic, cool circumneutral water body, that was also rich in vegetation (D’Ambrosio et al., 2020; Verschuren and Eggermont, 2006; Wojewódka et al., 2016). However, the high abundance of cladoceran ephippia (sexual reproduction), which are indicators of unfavorable conditions for these microcrustacean communities—such as a decline in water level—implies that the flooding event was short-lived (Amsinck et al., 2007). This interpretation is consistent with the lowest values of Ti/Ca, Zr/Fe, TN, and TOC, all pointing to a drier phase around 4.7 ka (Figure 5). The short-lived nature of this event is also compatible with the terrestrial record, which shows a brief reversal in pollen assemblages, supporting the interpretation of a basin-wide environmental shift rather than an isolated limnological fluctuation.
Subsequently, between 4.3 and 3.5 ka, the Ti/Ca, Zr/Fe, P, TN, and TOC values increase, indicating enhanced flooding and greater organic matter accumulation, consistent with the establishment of a more stable water body. This period, referred to as Pond-1, corresponds to a distinct cluster in the MFA analysis (Figure 4a and b). During this stage, facultative planktonic diatoms (Aulacoseira nivaloides, Cavinula pseudoscutiformis) and benthic taxa (Psammothidium helveticum) were dominant, indicating moderately humid conditions with colder winters (Caballero et al., 2020). Once again, this assemblage closely resembles that of modern Lake El Sol, indicating an oligotrophic, low-conductivity, near-neutral (pH >6), cool, and permanent water body. At the same time, Oedogonium and Cosmarium reached their maximum abundance in the sequence, the latter being a pioneer and colonizing alga (Huber, 1996). Chironomids such as Paratanytarsus—an indicator of increased aquatic vegetation density (Brooks and Heiri, 2013; Tarrats et al., 2018)—and cladocerans (Pleuroxus) reappeared. By ~3.9 ka, an increase is observed in the abundance of Spirogyra, Alona, and Pleuroxus, along with a higher frequency of Heterocypris cf. incongruens, a littoral ostracod characteristic of shallow and seasonal aquatic habitats (less than 1 m deep; Coviaga et al., 2018; Rossi et al., 2011), as well as testate amebae. Together, the biological and geochemical indicators point to the establishment of a shallower pond.
The period from 3.8 ka to the present (Pond-2) is characterized by a marked increase in TOC, TN, P, and in the detrital (Ti/Ca) and flooding (Zr/Fe) indicators, albeit with fluctuations (Figure 5). Algal taxa such as Gymnodinium, Debarya, Spirogyra, and Zygnema reached their highest abundances, together with the chironomids Paratanytarsus and Heterotrissocladius and testate amebae. These proxies collectively reflect the development of a pond ecosystem with well-established trophic interactions, involving primary producers (diatoms and algae) and herbivorous zooplankton (cladocerans and/or ostracods), which in turn supported higher trophic levels represented by aquatic macroinvertebrates (chironomids). The observed inverse relationship between diatom abundance and chironomids and cladocerans may partly reflect grazing pressure, as both groups are known to consume diatoms. However, this pattern is likely controlled by multiple interacting factors, including water-level fluctuations and nutrient availability, which can simultaneously influence aquatic communities and sedimentary proxies. During this stage, several fluctuations were observed, within which three probable drought intervals were identified. These are supported by the convergence of aquatic biological proxies, geochemical indicators of reduced flooding and detrital input, and changes in the terrestrial record consistent with increased environmental stress.
The first event, between 3.2 and 2.7 ka, is marked by the presence of ephippia and a notable decline in the algal assemblage (Gymnodinium, Debarya, and Zygnema). The diatom assemblage was dominated by aerophilous and benthic species, and the record of microcrustaceans almost disappears (Figure 7). Taken together, these aquatic proxies indicate a marked reduction in water availability and ecosystem instability affecting both aquatic communities and the surrounding vegetation, likely reflecting periods of reduced effective moisture and/or threshold responses of a shallow pond system to climatic variability.
The second event, between 2.0 and 1.9 ka, shows a decrease in the total diatom concentration (Figure 7), dominated by aerophilous species such as Stauroneis obtusa, Pinnularia borealis var. subislandica, and Pinnularia sinistra, suggesting waterlogged soils and cold conditions. At this time, the concentrations of Debarya and Zygnema also decline. Following the second drought event (1.8 ka), the record indicates the occurrence of Heterotrissocladius, a genus of cold stenothermal chironomids associated with low-nutrient water bodies (Hamerlík et al., 2022). This pattern corresponds with its peak abundance observed as phosphorus (P) concentrations began to decline around 1.5 ka and could be related to the Late Antiquity Little Ice Age (LALIA; ca 1.5–1.1 ka; Figure 3). The consistency between low diatom productivity, reduced algal abundance, and the appearance of cold-adapted chironomids supports the interpretation of a shallow, cold, and stressed aquatic environment during this interval, consistent with reduced water availability and/or internal ecosystem responses to climate variability.
The third pulse of dry conditions occurred between ca. 1.2 and 1.1 ka, coinciding with the period of the so-called Maya drought. During this event, the aerophilous species Stauroneis obtusa, Caloneis arophila, Pinnularia borealis var. subislandica and P. shimanskii dominated, indicating low water levels, exposure of the substrate and cold conditions (Hejduková et al., 2024; van Kerckvoorde et al., 2000). The absence of organisms such as ostracods and chironomids also support the interpretation of a renewed phase of reduced water levels in the pond, which may reflect either pronounced drought conditions or the sensitivity of the system to relatively modest hydroclimatic changes. (Figure 7). This signal is consistent with the broader multiproxy pattern of hydrological contraction and environmental stress recognized in the Late-Holocene section of the sequence.
During the Little Ice Age (AD ~1300–1850), diatom abundance remained low; however, facultative planktonic, benthic, and aerophilic species such as Alaucoseria nivaloides, Cavinula pseudoscutiformis, Psammothidium helveticum, Caloneis aerophila, and Pinnularia borealis var. subislandica reappeared, indicating that water levels were higher than during the Maya drought. Nevertheless, these conditions did not allow the establishment of zooplankton and chironomid communities.
Around AD 1550 (5.5 cm), conditions returned to a stable pond environment where algae such as Spirogyra, Pleuroxus, Paratanytarsus, and Heterotrissocladius proliferated, and Heterocypris cf. incongruens reappeared, indicating a brief period during which the pond filled with water. The aerophilous diatom species Chamaepinnularia sp. and Pinnularia borealis var. subislandica indicate low to neutral pH and low organic content. The testate amebae species (Difflugia urceolata urceolata, D. globulosa, and D. oblonga oblonga; Figure 2 in Supplemental material), which are associated with a more disturbed environment, imply human presence during this period. This evidence suggests that the event likely resulted from the construction of a canal by the indigenous people to keep their shrine permanently supplied with water.
Holocene climatic variability in the high-mountain region of central Mexico
High-mountain records are particularly sensitive to climate change (Diaz et al., 2014), and the Nahualac sequence provides valuable information on climatic variability and landscape changes over the last 8.4 ka. The increase in global temperature for the Holocene began around ca. 11.7 ka, and although there was a general temperature rise throughout 6 ka, a cooling trend became evident for the Meghalayan (Liu et al., 2014). Moreover, superimposed on this general warming trend, temperature variations at millennial scales during the Holocene have been described based on compilations of palynological records for North America (Viau et al., 2006).
During the Greenlandian, there was pronounced seasonality with high summer insolation in the Northern Hemisphere. Both the available records and climate models indicate that, from 9 ka onward, a summer-rainfall regime was established (Hardt et al., 2010; Metcalfe et al., 2015), causing increased precipitation across the northern tropics. In central Mexico, precipitation follows a monsoonal pattern related to the latitudinal displacement of the Intertropical Convergence Zone (ITCZ; Mosiño-Aleman and García, 1974). Although the general pattern points to enhanced summer rainfall during the Greenlandian —associated with a northern ITCZ position (Medina et al., 2023) and a strengthened AMOC (Ayache et al., 2018) some lacustrine paleorecords from central Mexico suggest a negative hydrological balance during this period (Figure 8). However, we consider that the combined effect of increased humidity relative to the end of the Pleistocene and higher-than-present temperatures likely enhanced evaporation in the high-elevation basins of central Mexico (~2200 m asl). This has been documented in Chalco, where increased carbonate precipitation, lower lake levels and saline conditions have been reported (Caballero et al., 2019, 2025; Sedov et al., 2010).
Summary of the main variables in terms of paleoecological information of the Nahualac-18-II sequence compared with the following: the AMOC index reconstruction for summer (green), winter (blue) and annual (black; Ayache et al., 2018); the ITCZ variability based on δ18O data from Alfredo Jahn cave (VEAJ), north-central Venezuela (Medina et al., 2023); the precipitation seasonality reconstruction from δ18O data of Buckleye Creek cave, west Virginia, USA (Hardt et al., 2010); and the seasonality insolation data at 19°N (Berger, 1992).
Temperature fluctuations based on the glacial history of Iztaccíhuatl provide a climatic framework for the paleoenvironmental interpretation of the Nahualac record. (Vázquez-Selem and Heine, 2011; Vázquez-Selem and Lachniet, 2017). Milpulco-1 glacial advance took place between 13 and 10.5 ka, with an estimated temperature drop of 5.4°C–4.4°C compared to modern (late 20th century) values. At this time a moraine indicates that the local glacier front was at 3900 m asl, that is, the same elevation of the Nahualac site and only 300 m away from it (Figure 1c). The site was likely a barren area under marked periglacial conditions.
In the early Northgrippian the Milpulco-2 advance (8.5–7.5 ka) is associated with temperatures 3.3°C–4 °C below modern values. The glacier front at the peak of this advance (ca. 8.5 ka) was only 700 m away from Nahualac and 70 m higher, at 3970 m asl. The equilibrium line altitude (ELA) of the glaciers of the mountain (i.e. the elevation of the 0°C isotherm) was at around 4400 m asl (Vázquez-Selem and Heine, 2011). Based on the difference in elevation of ca. 500 m between the ELA and Nahualac, a mean annual temperature of 3°C to 3.5°C can be estimated for the pond at that time, that is, ca. 3°C below modern values.
The first millennium of the sedimentary record of the pond is coeval to the gradual recession of the Milpulco-2 glaciers. Mean temperature anomalies derived from transfer functions from palynological data in mid-altitude regions, such as Lake Chalco, indicate an increase of 0.5°C–1°C above present values following the Milpulco-2 advance (Correa-Metrio et al., 2013). This warming in both high- and mid-altitude areas of the Central Mexico region is also evident in a multiproxy record from the southern margin of Lake Chalco (Tulyehualco), marked by a lake-level drop, the expansion of subaquatic vegetation, and major fire events between 7.3 to 5.5 ka (Caballero et al., 2025). The set of records from Central Mexico shows a temperature increase between 6.5 and 5.5 ka, coinciding with global reconstructions of mean surface temperature, where the warmest millennium of the Holocene centers around 6.5 ka (Kaufman and Broadman, 2023). This is in agreement with an upward shift of the treeline above the Nahualac site around that time. From 6.6 ka onward, after the treeline surpassed the elevation of Nahualac, the mean annual temperature around the site was likely >6.7°C ± 0.8°C, which is the temperature associated with the treeline (Körner and Paulsen, 2004).
At Nahualac, the data indicates climatic control on fire regimes typical of high-elevation tropical environments, where the duration of seasonal droughts and the continuity of fuel determine combustion potential. During the early–middle Holocene (8.4–5.5 ka), the absence of local fire peaks and the slightly higher proportion of herbaceous particles suggest humid conditions with frequent cloud cover, which would have limited ignition and fire spread. This pattern is consistent with a wetter early Northgrippian in the northern tropics, related to a more northerly ITCZ position and intensified summer precipitation (Figure 7; Hardt et al., 2010; Medina et al., 2023).
A transitional period between 5.5 and 4.5 ka is detected in Nahualac, marked by a rise in the pond level followed by an abrupt decline, indicating probable drought conditions around 4.7 ka. These fluctuations may be linked to changes in seasonality, with the gradual reduction in summer insolation followed by significant north–south shifts in the mean position of the ITCZ (Medina et al., 2023), as well as variations in monsoon and in ENSO intensities (Hardt et al., 2010; Lu et al., 2025; Figure 8).
The interval 4–1.3 ka represents the peak phase of local fire activity. The higher proportion of samples with positive Z-scores and elevated CHAR values (~ 1483 mm2/cm2/ year) reflect an increased frequency of burning episodes. Throughout this period, woody charcoal dominates, and the most intense peaks coincide with higher proportions of woody particles. The positive (although weak and statistically non-significant) correlation between Z-scores and woody fractions (r = 0.21, p = 0.34) together with the negative with herbaceous fractions (r = −0.27, p = 0.21) is consistent with the visual associations between major CHAR peaks and increased woody biomass, even if the correlations alone cannot be used as proof of fire-fuel coupling. Taken together, the multiproxy evidence strongly suggests that major fire episodes were primarily fueled by accumulated woody biomass under enhanced seasonality and prolonged dry season. Such a scenario agrees with the southward ITCZ displacements and Atlantic Meridional Overturning Circulation (AMOC) variability, both of which favored drier conditions in tropical Mesoamerica during the Meghalayan (Chiessi et al., 2021).
Although the multiproxy dataset for the last 4.2 ka indicates the establishment of a pond, conditions were far from stable, as evidence suggests the occurrence of three probable drought intervals centered at ~ 3.2–2.7, ~2.0–1.9, and ~1.2–1.1 ka. Comparisons between these droughts and evidence from other lakes in Central Mexico are complex as Late-Holocene sediments records from the mid-elevation sites in the region are generally disturbed by past and present human activity (Caballero et al., 2025). Another source of information on rainfall variability in the region during the last 2.4 ka comes from a speleothem in Juxtlahuaca Cave, in south-central Mexico; variations in δ18O are interpreted as a proxy for summer monsoon intensity and have been linked to Mesoamerican cultural change (Lachniet et al., 2012). The 2.0–1.9 ka drought at Nahualac, which likely falls within the broader Meghalayan pyroclimatic window, is consistent with the visual association between major CHAR peaks and increased woody biomass, even if the correlations alone cannot be used as proof of fire-fuel coupling. Taken together, the multiproxy evidence strongly suggests that major fire episodes were primarily fueled by accumulated woody biomass under enhanced seasonality and prolonged dry seasons correlates with low δ18O values corresponding to the Early Teotihuacan dry period. This early Classic drought has been documented in other sites such as Lake Coatetelco (Caballero et al., 2023). The ~1.2–1.1 ka drought broadly coincides with the first pulse of the Maya drought.
The Nahualac record reveals a maximum in fire activity between 4.0 and 1.3 ka, interpreted as a pyroclimatic window associated with drier conditions, flanked by periods of low fire activity during both the early Northgrippian and the Meghalayan. In the uppermost part of the record, the relative increase in herbaceous particles may indicate the onset of incipient human disturbance.
Fire activity declined markedly during the last millennium (1.3 ka), while the composition of fuel shifted toward higher herbaceous proportions. This reduction in fire frequency coincides with cooler conditions associated with the Little Ice Age, which would have limited ignition and fire spread in mountain ecosystems. Although herbaceous particles increase slightly (~15%), charcoal evidence does not indicate distinct high-intensity fire peaks attributable to anthropogenic burning. Instead, the record suggests a gradual attenuation of fire activity under cooler conditions, accompanied by incipient signs of human disturbance such as vegetation opening. The most recent sample (AD 1600) suggests a slight landscape opening, possibly linked to early human activities such as grazing or small-scale burning, as documented in other highland basins of central Mexico (Caballero et al., 2020).
During the Little Ice Age (LIA), glaciers on Iztaccíhuatl advanced, suggesting a snowline lowering of ca. 250 ± 70 m and a cooling of 1.5°C–1.9°C, with likely impact on the treeline. Above Nahualac, the glacier reached ca. 4,400 m asl. The pollen record from Nahualac shows an increase in alpine grassland during this oscillation, consistent with Charnay’s (1885) sketch depicting sparse tree cover at the site (Figure 1 in Supplemental material), suggesting a lowered treeline elevation. A shallow pond with benthic and aerophilous diatoms and the presence of ephippia indicate a drier environment.
Conclusions
Upper mountain areas contain sedimentary archives that function as sensitive recorders of climate change and, in many cases, are less affected by human activities than those from lowland areas. The multiproxy integration of the Nahualac record (3999 m asl) provides evidence of climate variability and its impacts on this high-altitude landscape, including vegetation dynamics with forest expansion toward higher elevations, fire activity, and pond evolution associated with environmental and climatic changes during the last 8.4 ka.
Using the glacial history of Iztaccíhuatl as an environmental framework, the first millennium of sedimentation in the pond record is contemporaneous with the gradual retreat of the Milpulco-2 glaciers. The pollen results reveal changes in the composition of temperate forests and vegetation structure at upper elevations as responses to temperature and moisture variations during the Holocene. The main changes recorded include the establishment of alpine grassland between 8.4 and 7.8 ka, considered an environmental prelude to later forest expansion. The palynological assemblage suggests that the upward migration of forests occurred around 6.6 ka and continued until 5.2 ka, with significant changes in composition, such as the establishment of Abies religiosa (6.1 ka) at higher elevations than today, consistent with warming, and the presence of other mesic elements. However, this migration was not continuous, the return of alpine grasslands between 4.8 and 4.2 ka interrupted the trend, followed by the re-establishment of Pinus forests and an increase in fire activity. A Late-Holocene decline in Pinus with higher Poaceae around ~0.75 ka suggests a short dry episode affecting the upper elevational ecotone.
Fire regimes at Nahualac were strongly climate-dependent throughout the Holocene. The charcoal record shows that local fire activity was minimal during the early Northgrippian, when cool, moist conditions limited ignition and maintained open alpine landscapes. A sustained accumulation of woody fuels accompanied the expansion of Pinus and Abies, toward higher elevations, enabling a phase of intensified burning between ~4.0 and 1.3 ka. This pyroclimatic window reflects enhanced seasonality and prolonged dry seasons associated with southward ITCZ displacements and AMOC variability conditions favorable for woody fuel drying and high-severity fire episodes. In contrast, the last millennium shows a marked decline in fire frequency and a shift toward herbaceous fuels, consistent with cooler and wetter conditions during the Little Ice Age and only incipient human disturbance near the end of the sequence. Overall, the Nahualac record demonstrates that fire behavior in high-mountain tropical ecosystems is governed primarily by climate-driven changes in fuel availability and type, with woody biomass accumulation acting as the key determinant of high-intensity fires.
The sedimentary record from Nahualac documents a progressive but highly dynamic development of a high-altitude aquatic ecosystem over the last 8.4 ka, closely controlled by hydrological balance and climatic variability. During the early Northgrippian (8.4–5.5 ka), multiproxy evidence indicates the persistence of an oligotrophic wetland under cool, moist conditions, characterized by low productivity, aerophilous diatom assemblages, and sporadic cold-water algae, consistent with shallow, flooded soils rather than a permanent lake.
Between ~5.5 and ~4.5 ka, Nahualac experienced a transitional phase marked by short-lived flooding events followed by drying, reflecting unstable hydrological conditions. A brief establishment of a cool, oligotrophic permanent water body around ~4.9 ka was rapidly interrupted, as indicated by ephippia production and low geochemical inputs, highlighting the sensitivity of the system to short-term climatic shifts.
After ~4.3 ka, a more persistent pond developed, initially under moderately humid and cool conditions (Pond-1), followed by a shallowing trend and, after ~3.8 ka, the establishment of a more complex pond ecosystem (Pond-2) with more complex trophic interactions. Despite this apparent stabilization, the Meghalayan was punctuated by recurrent probable droughts intervals (~3.2–2.7, ~2.0–1.9, and ~1.2–1.1 ka), during which water levels dropped sharply and aquatic communities collapsed or reorganized, reflecting strong sensitivity to regional hydroclimatic forcing.
During the Little Ice Age, cooler conditions limited biological productivity and prevented the re-establishment of stable zooplankton communities, despite episodically higher water levels. The uppermost part of the record reveals brief anthropogenic modification of the hydrological system around AD 1550, likely associated with canal construction to maintain water availability.
The Nahualac record confirms that high-mountain environments in central Mexico responded sensitively and coherently to Holocene climate forcing at both orbital and millennial scales. Early Northgrippian conditions were shaped by cool temperatures, strong seasonality, and enhanced summer precipitation linked to a northerly ITCZ position and a strengthened AMOC, which together limited fire activity and maintained humid conditions despite potentially high evaporation at high elevations.
Northgrippian warming, particularly between ~6.5 and 5.5 ka, promoted glacier retreat, the expansion of forest communities toward higher elevations and the accumulation of woody biomass. This interval corresponds with regional and global temperature maxima and marks a key threshold after which fire regimes became increasingly sensitive to fuel continuity and seasonal drought. The subsequent Meghalayan was characterized by greater hydroclimatic instability, including a transitional phase around 5.5–4.5 ka and the development of a pronounced pyroclimatic window between ~4.0 and 1.3 ka, when fire activity peaked under drier conditions associated with southward ITCZ shifts and AMOC variability and increased ENSO activity.
During the last millennium, cooler Little Ice Age conditions reduced fire frequency, favoring the expansion of alpine grassland and a shallower pond system. Although subtle signals of human disturbance appear late in the record, the dominant controls on fire, vegetation, and hydrology throughout the Holocene were climatic.
Supplemental Material
sj-docx-1-hol-10.1177_09596836261450779 – Supplemental material for Multi-proxy paleoecological analysis of a high-altitude site in central Mexico: Holocene vegetation dynamics, pond development and climate variability
Supplemental material, sj-docx-1-hol-10.1177_09596836261450779 for Multi-proxy paleoecological analysis of a high-altitude site in central Mexico: Holocene vegetation dynamics, pond development and climate variability by Socorro Lozano-García, Diana Avendaño, Montserrat Amezcua-Vargas, Esperanza Torres-Rodríguez, Itzel Sigala-Regalado, Susana Sosa-Nájera, Erandi Rodríguez-Pérez, Margarita Caballero and Lorenzo Vázquez-Selem in The Holocene
Footnotes
ORCID iDs
Author contributions
Funding
The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors received financial support from the Universidad Nacional Autónoma de México for the research, authorship, and/or publication of this article.
Declaration of conflicting interests
The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Supplemental material
Supplemental material for this article is available online.
References
Supplementary Material
Please find the following supplemental material available below.
For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.
For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.
