Introduction

After European contact and the Spanish conquest of the territory representing modern Guatemala (1524–1541 ce) the anthropogenic use of land was drastically altered to facilitate the rearing of livestock, agriculture, large scale timber extraction, mining and new settlements. These were heterogenous across the uplands and lowlands varying in extent, duration and intensity (Thompson 1958). These practices initially attained to Spanish colonial aspirations following three important ideologies: (i) the church; (ii) the state; and (iii) the ambition of individuals to attain and generate wealth (Lovell 1985). The Cuchumatanes highlands (Fig. 1; ESM 1 Table 1) were subject to drastic land use changes under colonial rule documented in various historical accounts (e.g. Ingersoll 1826; Thompson 1958) and official records (e.g. Cortez y Larraz 1958). In the absence of chronicled accounts or direct measures of land use change, palaeoenvironmental reconstructions can be used to infer impacts and the mechanistic responses of anthropogenic practices across space and through time (e.g. Dull 2004a, b, 2007; Velez et al. 2011).

Fig. 1
figure 1

Regional Map of the Maya area depicting: (i) Columbian palaeo-climatic and palynological records; (ii) the current departments for the Cuchumatanes Highlands; and (iii) colonial ecclesiastical administrative areas and churches: (I) San Andres Cuilco; (II) San Marcos; (III) San Pedro Soloma; (IV) Chiantla; (V) Huehuetenango; (VI) Santa Ana Malacatan; (VII) Santa Maria Nebaj; (VIII) San Miguel Uspantan; and (IX) Totonicapan (Lovell 1985). Ecclesiastical administrative areas are depicted and linked to the local parishes through the connected dots

In this study we apply a multi-proxy approach to independently reconstruct the land use history of the area surrounding the Cuchumatanes highlands after European contact, with the aim of identifying how and when anthropogenic practices involving different types of land use change impacted vegetation composition during the past 500 years, and to independently validate the catchment area for burning represented by macroscopic and microscopic charcoal using satellite earth observation data.

Integration of remotely sensed satellite data with palaeoecological records is often challenging and rarely attempted due to the resolution of palaeolimnological records and their associated age-depth uncertainties (Metcalfe et al. 2015). Lakes with high sedimentation rates may be precisely dated using 210Pb measurements to give approximately annually resolved sedimentary sequences going back up to 200 years (Oldfield and Appleby 1984). Precision dating of lake sediments allows for palaeoecological data and satellite data to be directly compared enabling contemporary validation of the palaeoecological proxy data (e.g. charcoal representing burning) for as long as the satellite instrument has been active (e.g. MODIS since 2000 ce). We present the first sub-decadal dated palaeolimnological record from this region enabling a high-resolution integrated reconstruction of land-use history since the Spanish Conquest (1524 ce).

Climate

In the uplands of Guatemala, climate (precipitation and temperature) is highly variable and is primarily determined by the altitude rather than latitude (Ingersoll 1826). Mean annual temperatures range between 14 and 25 °C, while average rainfall across the region lies between 900 and 3,700 mm depending on the elevation (Kappelle 2006). The region is defined as being within the Central American Dry Corridor (CADC), a region of Central America defined by its high rainfall variability between the wet and dry seasons (Bouroncle et al. 2017). There is currently only one palaeoclimatic reconstruction from the Cuchumatanes Highlands, a high resolution isotope record from Cenote Kail, which suggests that conditions became increasingly wetter during the Little Ice Age (c. 1400–1850 ce) (Stansell et al. 2020). These findings are also corroborated by a nearby palaeoclimatic reconstruction using diatoms from Lake Amatitlán (Velez et al. 2011). The only other palaeoenvironmental precipitation record to have been produced from within the CADC is a δ18O lake record from lowland Central Pacific Nicaragua, 11.906° N, 85.918° W (Stansell et al. 2013). This record presents increasingly drier conditions for the Little Ice Age juxtaposed to the currently available records from the highlands. Little else is known about the spatial and temporal dynamics of precipitation for this upland region and therefore further research is warranted to determine if precipitation patterns are homogenous with the northern lowland regions (Fig. 1) or respond variably (Metcalfe et al. 2015).

Vegetation

Vegetation for this region of Guatemala is primarily comprised of mixed deciduous and coniferous forests. The most outstanding characteristic of this biome is the diversity of pine (> 115 Pinus spp.) and oak (> 150 Quercus spp.) species (Muller 1942; Kappelle 2006), which are well adapted to variable climatic conditions and natural fires (Corrales et al. 2015). These pine-oak forest formations often form intricate mosaics and complex successional interactions extending up into broad leaved cloud forests at higher altitudes (Kappelle 2006; Rzedowski 2006). This biome is currently threatened by agricultural expansion, logging, firewood extraction, forest fires and pests (Veblen 1978; Harvey et al. 2019). The vegetation response to landscape disturbance over the past c. 500 years is particularly important for understanding the response and resilience of the mixed hard wood forests (MHWF) and coniferous forests (CF), e.g. pine forests (PF) and mixed pine-oak forests (POF), to past and current management practices, particularly those involving natural regeneration and municipal-communal management (Ewel 1979; Miles 1987; Holder and Chase 2012). To assess the status of these forests, it is essential to distinguish the stages of succession that have led to its current state (Figueroa-Rangel et al. 2012). Analysis of the lake sediments from Cenote Kail prior to European arrival discern the transformation of the POF by pre-Columbian populations, by identifying agrarian practices and burning for lime production as the dominant drivers of vegetation change between 4000 and 1500 ce (Harvey et al. 2019).

Forest cover across the whole of Guatemala is reported to have declined at an annual rate of 1.2% between 1990 and 2005 ce (Holder and Chase 2012). There are very few remaining primary forests outside of the cloud forests in the Sierra de las Minas and several regions in the Petén district (Veblen 1976; Holder 2006; Holder and Chase 2012). Secondary forests account for more than 40% of the total forest area in the tropics (Brown and Lugo 1990; Corlett 1995; Willis et al. 2004) forming complex successional groups as they have expanded to recolonise former agricultural landscapes (Tucker 2008; Brady 2009; Redo et al. 2009; Bass 2010). In Guatemala, most of these secondary forests occur on municipal and communal lands, which the local human populations depend on for timber and fuel (Urrea 1995; Reyes 1998; Pira et al. 1999; Gibson 2001; van Kempen et al. 2009; Holder and Chase 2012). While the high expense of transporting timber out of the Cuchumatanes Highlands has previously precluded this area for commercial logging (Veblen 1974), a history of deforestation for agrarian expansion has led to many hectares of forests being felled.

Today, the vegetation surrounding the Cuchumatanes highlands is comprised of PF, POF, cloud forests, and agricultural/pastoral grasslands (Breedlove 1981; Ramı́rez-Marcial et al. 2001; Rzedowski 2006; Chiabai 2015; Franco-Gaviria et al. 2018). Of these three vegetation types the PF are the most widespread, dominated by Pinus oocarpa at the canopy level and Acacia pennatula, Myrica cerifera, and Vernonia canescens at the understory level (Franco-Gaviria et al. 2018). The formation of these PF is mostly associated with sustained anthropogenic disturbance relating to anthropogenic use of fire and extensive grazing (Ramı́rez-Marcial et al. 2001).

The impact of European contact

The legacy left by the Maya caused widespread changes to the forests within the Maya Area through burning, e.g. for lime production (Hansen 2001, 2012; Anderson and Wahl 2016), and agrarian practices involving Milpa cycling (Wilken 1971, 1987; Nigh 2008; Ford and Nigh 2009), particularly during the Late Pre-Classic (350 bce–250 ce) and Classic Periods (250–950 ce), which is when the Maya civilization expanded towards its climax (Ford and Nigh 2009). Accompanying the Spanish conquest (1524 ce), European diseases, slavery and warfare drastically reduced the native population, resulting in reduced indigenous anthropogenic pressure on the upland forests (Dull 2004a, 2007; Velez et al. 2011); however, this was counterbalanced by new pressures of colonial settlements (Jáuregui 1894; De Remesal 1932; Lovell 1985; van Oss 1986), grazing (Carranza 1897; Vázquez and Lamadrid 1944; Veblen 1976) and resource extraction (Recinos 1954; Sherman 1979; Lovell 1983) as indicated by historical records and accounts. In the Cuchumatanes Highlands the native population was reduced from 260,000 to 73,000 between 1520 and 1550 ce (Lovell 1985).

Previous palynological reconstructions linking European contact to the environmental changes in the Maya area are limited; most available records are of low temporal resolution, greater than 200 years between samples (e.g. Deevey et al. 1979; Curtis et al. 1998; Mueller et al. 2009) or are missing sediment from the last c. 500 years (e.g. Neff et al. 2006; Caffrey et al. 2011). One palynological study has been previously conducted south of the Cuchumatanes highlands (Lago Amatitlán), which reports evidence for deforestation and increased agriculture from c. 1600 ce, and greater tree diversity after 1875 ce, including increasing abundance of Alnus, Eugenia and Morella (Velez et al. 2011); and two studies to the north (San Lorenzo, and Esmeralda), which report declines in Pinus and Alchornea attributed to a reduction in environmental and anthropogenic pressure on the vegetation (Franco-Gaviria et al. 2018).

Across the wider Maya area (Yucatan Peninsula, Belize, Guatemala, El Salvador, and Northern Honduras), 29 sites have been previously investigated presenting palynological reconstructions spanning from c. 1500 ce through to the present (Fig. 1). Evidence for widespread forest recovery is supported after the initial arrival of the Spanish (c. 1500–1600 ce) across lowland Guatemala (Brenner et al. 1990; Anderson and Wahl 2016; Wahl et al. 2016), Mexico (Leyden 1987) and El Salvador (Dull 2004a, 2007). However, greater land clearance and deforestation are reported from Belize at this time (Walsh et al. 2014). This widely reflects the reductions of indigenous populations and corresponds to reduced burning across the larger Maya area (Dull 2004a, b; 2007; Walsh et al. 2014; Anderson and Wahl 2016). After initial forest recovery, few palynological studies report upon vegetation changes during the past 500 years (e.g. Velez et al. 2011; Franco-Gaviria et al. 2018); however, Dull (2007) comments upon the reductions of Quercus during the nineteenth and twentieth centuries, which is supported by palynological evidence from Laguna del Espino (Dull 2004a), Laguna Cuzachapa (Dull 2007), and Lago Amatitlán (Velez et al. 2011). By integrating high resolution proxy data for vegetation, burning and herbivory with the historical record and modern remote sensing data, this study aimed to determine the long-term ecological implications of anthropogenic forcing upon the vegetation of the areas surrounding the Cuchumatanes Highlands since the Spanish Conquest (1524 ce).

Methods

Field and sampling techniques

Cenote Kail is situated in the uplands of Guatemala (16° 00′ 00.0″ N, 91° 33′ 14.4″ W; 1,534 m a.s.l.), surrounded by a mixed coniferous forest mosaic, ranches, agricultural fields and mines (Fig. 2). In 2015 a 105 cm composite core, with overlapping sections, was extracted from Cenote Kail using a combination of a Livingstone piston corer (Livingstone 1955) and a piston corer with polycarbonate tubing. The uppermost 20 cm of surface sediment were sampled in the field at 0.25 cm intervals into plastic Whirlpak® bags in order to preserve the sediment-to-water interface.

Fig. 2
figure 2

Aerial view of Cenote Kail depicting land cover: a forests, agriculture and pastures (Bing Maps 2019a, b). b Photograph looking over Cenote Kail exemplifying contemporary vegetation assemblage and anthropogenically modified land cover (Google Maps 2019a, b)

Thirty subsamples along the sedimentary sequence (1 g wet weight) were extracted for biological proxy analysis of fossil pollen, Sporormiella, macroscopic charcoal (> 150 μm) and microscopic charcoal (> 10 μm, counted on pollen slides) based upon information inferred from the age depth model. The upper 36 cm were sampled at a resolution of 2 cm representing 2–13 years between samples. This enabled a high-resolution comparison of the palaeoecological data to satellite derived data covering the past c. 20 years. The lower half of the core (41–105 cm) was sampled at 5–10 cm resolutions representing 25–50 years between samples.

Chronology

A high-resolution age depth model was constructed using ten 210Pb dates and four calibrated 14C dates (Table 1). Charcoal fragments were carefully selected for radiocarbon dating. Samples were pre-treated using standard acid–base-acid protocols (Abbott and Stafford 1996). Radiocarbon dates were generated at the W.M. Keck Carbon Cycle Accelerator Mass Spectrometry Laboratory at the University of California, Irvine. The IntCal13 calibration curve (Reimer et al. 2013) was used to calibrate the measured radiocarbon dates and OxCal (v.4.3) was utilised to construct an age-depth model applying a Bayesian approach (Ramsey 2009). Outliers were identified using the general outlier model implementing an outlier probability of 0.05 (Ramsey 2008). Sedimentation rate was calculated using this age-depth model. Years ce were assigned to the sediment subsamples based upon this age depth model.

Table 1 Measured, calibrated and modelled 210Pb and 14C ages for Cenote Kail

Fossil pollen and Sporormiella analysis

Fossil pollen was used to reconstruct the abundance and composition of past vegetation dynamics. Fossil pollen extraction followed standard palynological procedures applying the Oxford Long-Term Ecology Laboratories protocol (OxLEL 2016). Silicone oil was used to mount pollen grains on slides to allow for rotation, easing identification. Samples were spiked with known concentrations of an exotic marker, Lycopodium spores (batch No. 20848 or 9666), to calculate pollen influx and accumulation rates (Stockmarr 1971). Counting and identification of pollen grains were conducted at 400 × and 1,000 × magnification. For each level a minimum of 300 terrestrial pollen grains were counted. Morphological identification was achieved using: (i) pollen databases (APSA Members 2007; Bush and Weng 2007; Martin and Harvey 2017); (ii) published plates (Roubik and Moreno 1991; Willard et al. 2004); and (iii) botanical reference materials from the Oxford Long-Term Ecology Laboratory reference collection (ESM 1 Table 2; ESM 2). To interpret the relative composition of the forest (deciduous/coniferous), coniferous and mixed-hard-wood forest canopy taxa were compared as a ratio.

The abundance of coprophilous fungal spores (Sporormiella) was used to indicate herbivorous animal presence and abundance. Sporormiella spores were counted and morphologically identified on the same slides (Davis and Shafer 2006; Baker et al. 2016). The abundance of spores in sediments was calculated as accumulation rates (spore cm−2 year−1) using the spore concentration (spore cm−2) and the sedimentation rate (cm year−1) (Maher 1981; Bennett 1994; Bennett and Willis 2001; Baker et al. 2016).

Pollen counts were converted to percentages based on the terrestrial pollen sum. The annual influxes of Sporormiella, macroscopic and microscopic charcoal were calculated by using the area in mm2 and number of microscopic and macroscopic charcoal fragments, respectively, the number of Sporormiella spores, the exotic marker spores, sample volume, and the sedimentation rate (Maher 1981; Bennett 1994; Bennett and Willis 2001; Whitlock and Larsen 2002; Baker et al. 2016). To identify discrete zones in the resulting palynological diagrams, constrained hierarchical clustering upon the palynological assemblage was applied and statistical significance of zones was assessed using the broken stick model (Bennett 1996). Presentation of the palynological data was preformed using the packages rioja (Juggins 2009) in base R (R Core Team 2012).

Ordination analysis of the palaeoecological data

Before performing all ordinations analysis, we square-root transformed the palynological percentage data to normalise the distribution. A square root transformation was chosen because it can be applied to data sets containing zero values (Malik and Piepho 2016). To check if it was appropriate to apply a linear ordination method (e.g. Principal Component Analysis) or unimodal ordination method (e.g. Correspondence Analysis), detrended correspondence analysis (DCA) was conducted upon the palynological assemblage data (ter Braak and Prentice 1988). To calculate the species turnover we extracted the scores for the first axis of the DCA. Principal component analysis (PCA) was used to infer similarities between samples and the change in trajectories of composition of taxa through time applying a singular value decomposition of the centered but not scaled data matrix. This method of ordination was deemed appropriate to infer associations between taxa; however, limitations for the application of this method upon data sets containing zero values should be noted (ter Braak and Prentice 1988). Finally, we performed a canonical correspondence analysis (CCA) to quantify the relationship between environmental variables (macroscopic charcoal, microscopic charcoal and Sporormiella) and the palynological assemblage data. Ellipses representing the discrete zones were calculated using standard parameterization (cos(θ + d/2), cos(θ − d/2)), where cos(d) is the correlation of the parameters (Murdoch and Chow 1996). Statistical analysis and presentation of data was performed using the packages vegan (Oksanen et al. 2015) and rioja (Juggins 2009) in base R (R Core Team 2012).

Fire analysis: integrating fossil charcoal with satellite images

Past fire was quantified by counting and analysing fossil charcoal fragments. The area (in mm2) of microscopic charcoal (> 10 μm and < 150 μm) fragments was quantified on the same slides applying the point counting method at 400 × magnification and counts were recorded until a minimum of 50 Lycopodium spores and 200 fields of view were encountered for each level (Clark 1982). Macroscopic fossil charcoal fragments were counted from the > 150 μm fraction of each 1 g sample at 10 × magnification. Charcoal volume was calculated from loss on ignition density and weight of each sample. These data were presented as charcoal influx (mm2 or particles cm−2 year−1) as well as concentration (mm2 or particles per cm−3). Macroscopic charcoal is typically used to infer local burning (Gavin et al. 2003; Lynch et al. 2004; Petersn and Higuera 2007; Higuera et al. 2007, 2011; Anderson and Wahl 2016), while microscopic charcoal is typically used to infer more regional burning (Clark 1988). To independently validate this assumption at Cenote Kail, we followed the approach of Adolf et al. (2018a), where micro- and macroscopic charcoal influx were compared to the amount of fire pixels detected by Moderate Resolution Imaging Spectroradiometer (MODIS) sensors in increasing radii around the study sites. We increased the radii encompassing the MODIS active fire pixels in 1 km steps from 1 to 25 km radius and in 10 km steps from 30 to 200 km radius around Cenote Kail. As opposed to the approach in Adolf et al. (2018a), we do this validation for a single site which enables us to focus on the particular effects basin properties of Cenote Kail may have on charcoal source area. MODIS data was downloaded for the years 2000–2015 ce (MCD14ML version 6.1 Global monthly fire product, Justice et al. 2002, 2006). A MODIS active fire location represents the centre of a 1 km pixel that is flagged by the algorithm as containing one or more fires within the pixel size of fire detected (Giglio et al. 2003). MODIS data were taken at 3 year averaged moving windows to account for age-depth uncertainty in the sedimentary record (Giglio et al. 2003). Filtering of MODIS fire time series data was conducted using R to remove false fires (e.g. hot factory roofs) using an urban fire filter from Globcover 2009 (Arino et al. 2009). Spatial and temporal autocorrelations in the MCD14ML product were removed following Giglio (2013), Oliveras et al. (2014) and Adolf et al. (2018a).

Results

Chronology and resolution

The age-depth model indicates that the recovered composite core (105 cm) covers the past 500 years in continuous time (Fig. 3). The general outlier model did not identify any dates to be removed from the sequence. Overall model agreement was high (Index = 143.6) and represented an average age range of 206 years for the radiocarbon dates and 15.1 years for the 210Pb dates surrounding the median modelled age.

Fig. 3
figure 3

Age depth model of Cenote Kail post European contact

Palaeoecological trends

Three statistically significant zones were identified from the broken stick model. Zone 1 (1550–1725 ce) represents an average temporal resolution of 50 years between samples spanning 225 years (with a range of 49–50 years), Zone 2 (1725–1940 ce) has an average resolution of 17 years spanning 215 years (with a range of 3–50 years) and the temporal resolution of Zone 3 (1940–2015 ce) represents an average of 5 years spanning 75 years (with a range of 2–11 years). Seventy-six taxa were recognised in the palynological sequence, which were grouped based upon their lowest identifiable common taxonomic denomination. Pinus, Liquidambar, Quercus, Morella cerifera, Anacardiaceae and Fabaceae dominate the tree component while Asteraceae, Poaceae and Cyperaceae are the most abundant herbaceous taxa (Fig. 4).

Fig. 4
figure 4

Palynological percentage diagram of taxa appearing in an abundance greater than 2% showing: Forest structure; Coniferous to hardwood ratio; Pollen influx; Sporormiella influx; macroscopic and microscopic influx; DCA axis 1

The PCA displays several gradients and associations between samples and taxa (Fig. 5). The first axis represents 16.3% of the variation, while axis 2 represents 14.1% of the variation. To help visualise the most distinctive gradient we superimposed the independently calculated palynological zones upon the quadrants of the PCA. Zone 1 is most associated with the bottom left quadrant along with canopy taxa Quercus, Liquidambar and Bursera simaruba and understory taxa Anacardiaceae, Fabaceae and Roupala. Zone 2 spans the right top and right bottom quadrants represented by the canopy taxon Alnus; understory taxa Myrica and Pourouma; and herbaceous taxa Poaceae, Brassicaceae, Cyperaceae, Apiaceae and Primulaceae. Zone 3 is best represented by the top left quadrant associated with the canopy taxon Abies; understory taxa Terminalia, Rubiaceae, Aphelandra and Juniperus; and the herbaceous taxon Asteraceae. Results from the CCA show that microscopic charcoal and Sporormiella are significant environmental variables most associated with Zone 2 (Fig. 5). Macroscopic charcoal is most associated with Zone 3 but is not a statistically significant.

Fig. 5
figure 5

Canonical correspondence analysis (a). Principal component analysis (b). Zones derived from the broken stick model represented by ellipses at a confidence of 95%. Zone 1 = green triangles; Zone 2 = yellow squares; Zone 3 = blue circles

Zone 1 (105–74 cm, 4 samples, 1550–1725 ce) spans 175 years and is defined by canopy taxa Liquidambar (1.3–30.7%) and Quercus (0–13%). Total tree cover is very high (69.3–86.3%) with a dominance of canopy taxa from 1550 ce (82.7%); however, this dramatically decreases towards 1700 ce (28%) with understory (44.7%) and herbaceous (23.3%) taxa becoming more abundant from 1600 ce (Fig. 4). The reduction of Liquidambar (30.7–10.3%) from 1550 to 1600 ce followed by Quercus (12–0%) and a further reduction in Liquidambar (10.3–1.3%) from 1600 to 1700 ce marks the end of MHWF dominance in this sequence (Fig. 4). Fabaceae (35.7%) becomes prevalent around 1600 ce after which Anacardiaceae steadily increases through to 1700 ce (8.3–18.7%). From 1650 ce Morella cerifera rises to fill the tree component (5–17%) along with Pinus (14–20.3%). Alnus establishes in low abundance from 1700 ce (6.3%). Asteraceae (7.6–12.3%) and Poaceae (3.3–7.3%) increase from c. 1600 to 1700 ce. Capsicum is present in low abundances at 1550 ce (1.3%) and 1600 ce (2.7%). Sporormiella influx is in low and stable abundance (0.8–1.5 spores cm−2 year−1) until 1650 ce when it disappears. Macroscopic charcoal influx remains low (0.6–2.1 particles cm−2 year−1) and microscopic charcoal influx decreases from 1550–1700 ce (481.8–59.8 mm2 cm−2 year−1). Microscopic charcoal particularly decreases between 1550 and 1600 ce (481.8–154 mm2 cm−2 year−1). Pollen influx is highest between 1550 and 1650 ce, (201,170–173,958 grains cm−2 year−1) decreasing towards 1700 ce (80,513 grains cm−2 year−1).

Zone 2 (65–26 cm, 12 samples, 1725–1942 ce) covers 193 years and is characterised by the increasing abundances of understory taxa Morella cerifera (18.6–32%) and Anacardiaceae (0.6–8%) as well as the sustained abundance of herbaceous taxa Asteraceae (3.3–11.3%), Poaceae (3.3–16.3%) and Cyperaceae (0.7–7.7%). Total tree cover continues to slowly decline between 1750 and 1938 ce (73.3–68.7%). Increases in Euphorbiaceae (0–3.3%) and Myrica from 1860 to 1938 ce (0–3.7%) are joined by increases in Cyperaceae between 1860 and 1910 ce (0.7–7.7%). Zea mays punctuates the sequence at 1900 ce (0.3%), while Capsicum appears at 1910 ce (0.7%) and 1938 ce (0.3%). Sporormiella influx is very high between 1750 and 1800 ce (7.5–8.6 spores cm−2 year−1), absent between 1800 and 1900 ce (0%), and then high again between 1900 and 1938 ce (0–8.5 spores cm−2 year−1). Macroscopic charcoal (1.9–4.6 particles cm−2 year−1) and microscopic charcoal (115.2–160.5 mm2 cm−2 year−1) increase from 1700 to 1800 ce and both decrease between 1800 and 1900 ce (macroscopic 4.6–1 particles cm−2 year−1; microscopic 160.5–75.6 mm2 cm−2 year−1). Pollen influx remains stable (c. 174,604 grains cm−2 year−1).

Zone 3 (24–0 cm, 14 samples, 1942–2014 ce) encompasses 67 years and is defined by the increased abundance of the herbaceous taxon Asteraceae (17–26.3%). Total tree cover reduces from 1948 to 2014 ce (65.3–53%). Canopy (28.9%), understory (31%) and herbaceous (33.2%) taxa are, on average, in equal abundance. Morella cerifera decreases in abundance after 1997 ce (20–12.3%), while Liquidambar (0.6.3%) and Anacardiaceae (2–9.3%) increase from 1948–2003 ce. Macroscopic charcoal remains in low abundance; however, it increases from 1988 to 2006 ce (0.2–3.2 particles cm−2 year−1) and also peaks in 2010 ce (2.1 particles cm−2 year−1). Microscopic charcoal decreases from 1948 to 1988 ce (201.5–2.4 mm2 cm−2 year−1) and is almost absent from 1988 to 2014 ce (2.4–0.01 mm2 cm−2 year−1). Sporormiella remains in high abundance from 1948 to 1968 ce (1.8–6.5 spores cm−2 year−1) after which influx decreases and is absent after 2008 ce. Pollen influx also decreases after 1968–2014 ce (196,916–5,855 grains cm−2 year−1).

Charcoal area validation

Linking the charcoal influx values to MODIS-inferred modern day burning regimes, an approach used by Adolf et al. (2018b), showed that the source area for macroscopic and microscopic charcoal is estimated to be 16 km for macroscopic charcoal and 180 km for microscopic charcoal, assuming that distance is indicated by the strongest correlations between charcoal and satellite data (Fig. 6). The strongest correlations between charcoal influx and the MODIS product are slightly stronger for microscopic charcoal (r: 0.76) compared to macroscopic charcoal (r: 0.72). Microscopic and macroscopic charcoal are not significantly correlated (p-value: 0.3) between charcoal data sets supporting the two different source areas.

Fig. 6
figure 6

Calibration of macroscopic charcoal influx in particles cm−2 year−1, microscopic charcoal influx in mm2 cm−2 year−1. rFN stands for the correlations coefficient r between MODIS-inferred fire number and macroscopic and microscopic charcoal influx, respectively. The spatial intervals chosen for comparison are 1–25 km radius in 1 km steps around Cenote Kail and from 30 to 200 km radius, in 10 km steps including the radii 75 and 125 km

Discussion

How and when did anthropogenic practices involving different types of land use change impact floral assemblage changes around the Cuchumatanes Highlands?

Three stages of compositional change can be discerned from the palynological assemblage since European contact in 1524 ce. These changes represent responses to anthropogenic activities in the Cuchumatanes highlands corresponding to pastoral activities, agriculture, timber extraction, mining, and new settlements. These three phases encompass: (i) the decline of MHWF, associated with the building of new settlements, agriculture and timber extraction for fuel (1550–1675 ce); (ii) pastoral expansions involving the rearing of livestock (1700–1800 ce); and (iii) the expansions of urban settlements and the increasing anthropogenic management of the surrounding land (1821–2015 ce).

Colonial settlement construction

Evidence from Cenote Kail indicates that land clearance for the establishment of new settlements around the Cuchumatanes Highlands began from at least 1550 ce and by 1650–1675 ce most hard-wood taxa (Liquidambar and Quercus) had been removed from the landscape (Figs. 4 and 5). The palynological evidence suggests that hard-wood taxa were removed due to land clearing for new settlements under a regime called ‘congregación’ or ‘reducción’ from 1537 ce (van Oss 1986), which involved forced removal of indigenous peoples from their villages into colonial settlements (Fig. 1). Historical records report that these settlements were mandated to have a church and administrative buildings, for which substantial amounts of timber were required (Jáuregui 1894; De Remesal 1932). Administrative towns, such as Purification Jacaltenango (meaning “wall of cabins” in Nahuatl) located 40 km from the study site would have required even larger quantities of timber for both construction and for fuel (Lovell 1985).

By the late 16th/early 17th century, historical (Thompson 1958) and palynological accounts (Velez et al. 2011) report rapid deforestation across Guatemala to meet the needs of these towns. Timber was extracted for buildings, reportedly a combination of pine (Thompson 1958) and oak (Dull 2004a; Velez et al. 2011); however, it would appear from the historical records that timber extraction was non-discriminatory (Thompson 1958). For example, the forests surrounding Lake Amatitlán (Fig. 1) were rapidly being felled to meet the needs of Antigua for which evidence is reported from historical chronicles (Thompson 1958) and from a palynological record from Lake Amatitlán (Velez et al. 2011).

Pastoral activities: Ranching

It is likely that as well as being preferentially extracted for wood, the increasing pressures of pastoral activities damaged the reproductive success of taxa such as Quercus through the consumption of acorns by cattle and sheep (van Wieren 1996; Kappelle 2006). Bovids, ovine and porcine (cattle, sheep and pigs) can digest fibrous content and plant lignin (van Wieren 1996), meaning that predation of seeds such as acorns would have dramatically increased with pastoral expansions. This would prevent further reproduction and recovery of the MHWF in this area (Kappelle 2006; Wallace et al. 2015). During the Pre-Columbian era Sporormiella influx into Cenote Kail was found to be higher than for most of the Colonial and Post-Colonial Periods (Harvey et al. 2019; Stansell et al. 2020). This is thought to be due to the higher abundance of native herbivores, such as white tailed deer (Odocoileus virginianus), the Baird tapir (Tapirus bairdii) and peccary (Tayassuidae), which were reduced through habitat loss and hunting after European contact.

During the early 17th century Spanish settlers began to expand into the countryside (MacLeod 1973), and by the end of the century the anthropogenic impact caused to the forests and soils was evident though the practices of excessive forest clearing and overgrazing (Veblen 1976). In 1699 ce and 1701 ce, the villages surrounding Antigua were forbidden from grazing or cultivating the surrounding hillsides by royal decree to attempt to preserve the depleted soils (MacLeod 1973). These edicts are suggested to be the first conservation laws to be implemented anywhere in Latin America (Veblen 1976). Substantial degradation to the underlying soil is first reported in the 19th century for the north-eastern half of the Department of Totonicapan (Carranza 1897), which was the leading sheep-rearing district of Guatemala over the previous two centuries (Vázquez and Lamadrid 1944).

Sporormiella evidence from Cenote Kail suggests that pastoral activities expanded from 1700 to 1800 ce around the Cuchumatanes Highlands, and while no clear change in the overall vegetation composition can be discerned, sustained absence of Liquidambar and Quercus is evident (Fig. 4). This could indicate the translocation of cattle and sheep from a neighbouring province such as the Department of Totonicapan over to the Department of Huehuetenango (Fig. 1). Between 1800 and 1870 ce the influx of Sporormiella diminished, which could indicate a similar edict of land conservation as that recorded for the north-eastern half of the Department of Totonicapan (Carranza 1897). One other such explanation could revolve around civil unrest pertaining to the Act of Independence of Central America in 1821 ce, which led to land use change and abandonment due to war.

Land-use management

Burning, locally and regionally, is limited and in line with the modern regime, which suggests that fire has been managed and controlled at least since European contact (Fig. 4). During the Pre-Columbian era (before 1524 ce) fire was identified as a significant driver of forest turnover in this region (Harvey et al. 2019); however, land-use change to incorporate extensive ranching may have reduced natural burning regimes and/or previous anthropogenic uses of fire, such as burning for Milpa or lime production.

The conservation of the MHWF taxa around the Cuchumatanes Highlands may therefore be achieved by continuing to limit the presence of fire in these forests, and by controlling or removing domestic bovids, ovine and porcine. This would allow for MHWF taxa such as Quercus and Liquidambar to re-establish over the course of 80–200 years (Kappelle 2006).

There is no evidence of deforestation from mining activities within the palaeolimnological record from Cenote Kail. The possibility that the Cuchumatanes Highlands might contain significant mineral wealth led to open-cast mines opening in Pichiquil, San Francisco Motozintla and Chiantla shortly after the Spanish conquest (Recinos 1954; Sherman 1979; Lovell 1983). The gold mines of Pichiquil and San Francisco Motozintla did not prove successful (Recinos 1954); however, success at Chiantla led to silver and lead extraction from 1537 ce (Sherman 1979). The mines of Chiantla produced modest quantities of silver and lead throughout the colonial period (Lovell 1983).

What is Cenote Kail’s charcoal source area?

The size of the source area for macroscopic charcoal and microscopic charcoal into a depositional environment is contentious. Macroscopic charcoal is widely suggested to reflect burning from a local to sub-regional source area between 0 and 10 km radius (Gavin et al. 2003; Lynch et al. 2004; Petersn and Higuera 2007; Higuera et al. 2007, 2011; Anderson and Wahl 2016), while microscopic charcoal is suggested to primarily reflect regional burning within a radius of up to several hundred kilometres (Clark 1988); however, recent validation of charcoal dispersal across Europe linked both macroscopic and microscopic charcoal influx to a regional source area of around 40 km radius (Adolf et al. 2018b). Our results from Cenote Kail indicate that the catchment area for macroscopic charcoal and microscopic charcoal varied by a degree of magnitude. Validation through the comparative integration of the charcoal records with MODIS-derived fire time series suggests two distinct biomass burning histories, the first reflecting extra-local to sub-regional burning at a radius of 16 km, and secondly a clearly regional burning signal from up to 180 km. The direct distance of Cenote Kail to the Pacific coast is around 160 km suggesting that the regional charcoal signal is indicative of burning inclusive of this wider region, encompassing anthropogenic settlements built in the Cuchumatan Highlands as well as lowland activities. The macroscopic charcoal signal is likely to reflect a more local signal, encompassing the modern pastoral, agricultural and mining towns of Aguacate, Yalambojoch and Yulaurel. The subtle differences between the correlation scores of macroscopic and miscroscopic charcoal are likely not enough to truly discriminate between local vs regional fire; however, they do provide a best-case scenario estimate.

Conclusions

The establishment of colonial settlements is evident from the palynological record from Cenote Kail as well as nearby Lake Amatitlán, indicating that around 150 years of intensive logging led to the permanent decline of the MTHFs. The recovery of these forests shows evidence of suppression since through the expansion of pastoral practices. Seed predation is suggested to be the dominant factor preventing MHWFs from re-establishing in the Cuchumatanes Highlands during the past 500 years. Evidence for edicts of conservation can be inferred from the Sporormiella record for different parts of the Cuchumatanes Highlands during the eighteenth and nineteenth centuries. No direct evidence for the impacts of mining have been inferred from our data sets.

The high sedimentation rate and precise 210Pb dating of Cenote Kail has enabled a rare validation of the palaeolimnological catchment area through cross comparisons of satellite data with the charcoal proxy data. This approach has enabled a calibration of the charcoal data sets to allow for a quantitative analysis of burning, comparable to modern day observations. Radial burning was identified at 16 km for macroscopic charcoal and 180 km for microscopic charcoal. To properly interpret the signals represented in the palaeolimnological record it is recommended that all future palaeolimnological work should attempt to fully understand the catchment area for each proxy represented using this sort of quantitative approach. The integrated approach of palaeoenvironmental techniques with historical records and remotely sensed satellite data has enabled a unique and high-resolution reconstruction of the flora, fire, and animal abundance of the areas surrounding the Cuchumatanes Highlands over the past 500 years, detailing anthropogenic impacts since European contact.