Fire-human-climate interactions in the Bolivian Amazon rainforest ecotone from the Last Glacial Maximum to late Holocene

1 Introduction

The Amazon rainforest is an expansive and vital biome that displays a diverse distribution of vegetation formations, including tropical rainforests, flooded grasslands, and dry savannahs, creating a mosaic of habitats essential for the region's unparalleled biodiversity. The Amazon Rainforest Ecotone (hereafter the ARF-Ecotone) of the southwestern Amazon Basin is a transitional landscape between the humid evergreen tropical forests (HETF), the seasonally flooded savannahs (SFS) of the Llanos de Moxos (Mayle et al., 2007), the upland (cerrado) savannahs, and the semi-deciduous tropical forests (SDTF) of Chiquitaniá. Fire plays an integral role in maintaining these vegetation transition zones, from fire-limited HETF to fire-prone SFS, SDTF and cerrado vegetation (Figure 1). Modern synergies between climate change (e.g., more frequent, severe droughts (Bowman et al., 2020) and human activity (e.g., increased deforestation for ranching, crop land expansion, and infrastructure development) threaten to accelerate forest die-back through positive-feedback loops by increasing future tree mortality, fire susceptibility, fuel loads, fuel connectivity and fire intensity (Doughty et al., 2015; Feldpausch et al., 2016; Le Page et al., 2017; Brando et al., 2020). These compounding factors have caused the ARF-Ecotone to become one of the most heavily exploited and most threatened regions of the Amazon (Cochrane, 2009; Le Page et al., 2017). This recent increase in fire activity is not just of regional significance, given that changes to the Amazon rainforest have been identified as one of the potential major future tipping points in the Earth's climate state (Lenton et al., 2019).

While fire is an ancient influence on the Earth system human agency has become increasingly important (Bowman et al., 2020). Today, humans are the primary cause of fire ignitions globally (Bowman et al., 2020; Pereira et al., 2022), and in fire-prone ecosystems including the SDTF and cerrado savannahs of the ARF-Ecotone, fires are more frequent when humans are in the landscape (Power et al., 2016; Le Page et al., 2017; Maezumi et al., 2018b; Brando et al., 2019). Indeed, as the keystone species of fire (Pyne, 2021), humans are the only species able to directly and deliberately manipulate fires through a variety of management practices (Levis et al., 2018; Bowman et al., 2020). Humans can shape fire directly, by manipulating the frequency, severity, and intensity of fires, and indirectly, by influencing the abundance and distribution of vegetation and through the introduction, cultivation, and domestication of plant species (Bowman et al., 2016, 2020; Levis et al., 2018; Hoffman et al., 2021; Furquim et al., 2023). For example, ethnographic studies demonstrate hunter-gatherers often influence vegetation structure, composition and diversity near human settlements (Posey and Balée, 1989; Politis, 1996; Denevan, 2006; ter Steege et al., 2013; Clement et al., 2015; Lins et al., 2015; Levis et al., 2018).

While the impact of humans on modern fire is substantial, the human use of fire to modify landscapes was practiced long before the colonization of the Americas (Boivin et al., 2016; Roberts et al., 2017). There is increasing evidence that in the Amazon, Indigenous traditional burning practices (hereafter Indigenous burning) influenced fire activity in mountain cloud forests (Loughlin et al., 2018), humid evergreen tropical forest (Bush et al., 2015, 2016; Carson et al., 2015; Maezumi et al., 2018a,b, 2022), semi-deciduous tropical forests (Bush et al., 2015; Carson et al., 2015), and cerrado savannahs (Carson et al., 2015; Brugger et al., 2016) as early as the mid-Holocene (~6500 cal year BP). The recent archaeological evidence documenting extensive human occupation in the Bolivian ARF-Ecotone region after ~10,500 cal year BP (Watling et al., 2018; Capriles et al., 2019; Lombardo et al., 2020), raises the possibility that human-fire linkages may have influenced fire activity in the region much earlier than previously proposed (Urrego et al., 2013; Carson et al., 2014; Bush et al., 2016; Maezumi et al., 2018b).

In light of the above, it is critical to obtain and compare multidisciplinary, direct proxy records of changes in past fire (hereafter palaeofire) regimes, vegetation change, human activity, and changes in key climate parameters in order to explore changing fire dynamics through time and their implications for future model projections. However, as a result of the complex fire-human-climate interactions and incomplete historical or long-term datasets (O'Connor et al., 2011; Bowman et al., 2020), the role of human-influences on palaeofires in the ARF-Ecotone remains largely unknown. While previous palaeoecological studies in the region have examined the relationship between basin characteristics (e.g., type/size) and fossil pollen archives (Smith et al., 2021), little is known about the relationship between basin characteristics and fossil charcoal records used to reconstruct past fire activity. We hypothesize that there is greater sensitivity of small lake basins to detect local-scale fire activity, thus smaller basins may be more sensitive to picking up Indigenous burning. To address this issue, we compare high-resolution (sub-centennial) Late Quaternary records from sites with varying basin characteristics. In addition to basin size, we explore potential drivers (e.g., climate, human) of the spatial and temporal variability of fire and vegetation change in the ARF-Ecotone. We implement a multi-proxy approach utilizing geochemistry, δ¹³C and δ¹⁵N stable isotopes, and macrocharcoal (>125 μm) analysis from a new ~24,000-year-old sediment core (Cuatro Vientos) that exceeds the earliest evidence of human activity in the region by >~13,500 years. These data are compared with two existing palaeovegetation studies of Laguna Chaplin and Huanchaca Mesetta palm swamp in Mayle et al. (2000), Burbridge et al. (2004), and Maezumi et al. (2015, 2018b), along with regional palaeoclimate (Baker et al., 2001; Whitney et al., 2011) and regional archaeological data (Maezumi et al., 2018b; Capriles et al., 2019; Lombardo et al., 2020; Capriles, 2023; Peripato et al., 2023) (Supplementary Table S4).

2 Materials and methods

2.1 Study area

Noel Kempff Mercado National Park (NKMNP) is a 15,230 km² protected reserve located near the southern margin of the Amazon basin in the department of Santa Cruz, north-eastern Bolivia (14°22′9.336″, 60°51′14.3994″; Figure 1) (Killeen and Schulenberg, 1998). The park has been designated a UNESCO world heritage site due to its exceptionally high beta (habitat) diversity (Killeen et al., 2003; Heyer et al., 2018). The biodiversity of the park is attributed to the high levels of habitat heterogeneity which incorporates five distinct ecosystems: humid upland forest, inundated and riparian forest, semi-deciduous and deciduous forest, upland savannah and savannah wetland (Killeen and Schulenberg, 1998; Killeen et al., 2003; Burbridge et al., 2004). While the park is considered to be largely undisturbed by modern anthropogenic land use, ceramics have been recovered from a soil pit in an interfluve ~25 km northwest of our study sites and abundant ceramics mixed with charcoal (dated to ~400–500 cal year B.P.) were recovered from Amazonian Dark Earths throughout La Chonta (Guarayos province) ~150 km west of NKMNP (Burbridge et al., 2004). To date, no archaeological surveys or excavations have been conducted within NKMNP though there are, increasing numbers of archaeological sites identified in the surrounding region (Figures 1, 9, Supplementary Table S4).

2.2 Climate

The climate of NKMNP is characterized by a tropical seasonal climate with intermediate rainfall between increasingly moist regions to the north and the drier areas to the south (Killeen and Schulenberg, 1998). Mean annual precipitation at NKMNP is ca. 1,400–1,500 mm per year, with a mean annual temperature between 25 and 26°C (Centro de Previsão de Tempo e Estudos Climáticos, Instituto Nacional de Pesquisas da Amazônia [CPTEC/INPA], 2023). There is a 3–5 months dry season during the austral winter, when the mean monthly precipitation is <30 mm (CPTEC/INPA) Precipitation falls mainly during the austral summer, originating from a combination of deep-cell convective activity in the Amazon Basin from the South American summer monsoon (SASM) and the Intertropical Convergence Zone (ITCZ) (Vuille et al., 2012; Novello et al., 2017; Orrison et al., 2022) (Figure 2). The SASM transports Atlantic moisture into the basin and corresponds to the southern extension of the ITCZ driven by seasonal variation in insolation (Bush and Silman, 2004; Vuille et al., 2012). When the ITCZ is located over the equator or further south, the northeast trade winds shift southward and intensify over South America. Diversion of the South American low-level jet (SALLJ) along the eastern foothills of the Andes delivers convective moisture from the Amazon basin to the La Plata Basin. The SALLJ converges with the western edge of the South Atlantic Subtropical High (SASH) and south-south westerly mid-latitude flow to form the South Atlantic convergence zone (SACZ), bringing precipitation to south-eastern Brazil during the austral summer (December to February). During the austral winter, cold, dry polar winds from Patagonia, locally known as surazos, can move northward, penetrating the NKMNP region causing temperatures to drop to 5–10°C for several days at a time (Latrubesse et al., 2012; Mayle and Whitney, 2012).

2.3 Geomorphology and regional vegetation

NKMNP is located on the western reach of the Precambrian Brazilian shield. The eastern half of the park is dominated by the Huanchaca Plateau, a table-mountain ~600–900 m above sea level (a.s.l) comprised of Precambrian sandstone and quartzite (Killeen, 1998). The plateau is dominated by upland cerrado savannah vegetation (Killeen and Schulenberg, 1998). To the west lies a lowland peneplain, where the Precambrian bedrock is overlaid by Tertiary and Quaternary alluvial sediments (Killeen, 1998) which is dominated by HETF vegetation. The Río Iténez and Río Paraguá black-water rivers form the northeastern and western boundaries of NKMNP, respectively. Today, these rivers are lined by riverine evergreen tropical forests with patches of seasonally-inundated savannah in areas where soil drainage is poor (Smith et al., 2021). The southern border of NKMNP defines the modern ARF-Ecotone between the southern Amazonian HETF and the Chiquitano semi-deciduous tropical dry forest (SDTF) of eastern lowland Bolivia (Killeen et al., 2003; Smith et al., 2021).

2.4 Site descriptions

In this study we compare high-resolution (sub-centennial) Late Quaternary charcoal records from three sites in NKMNP with varying basin characteristics (Figure 1). The first site, Cuatro Vientos (CV; 14°31′18.5″S, 61°7′11.3″W; ca. 184 m a.s.l.) is a small-sized (~5 km² area) palm swamp located in the western portion of NKMNP. CV is dominated by HETF to the east and the riparian forest of the Río Paraguá to the west (Smith et al., 2021). Growing within the CV basin itself is a floating mat of sedge/grass swamp vegetation, interspersed with small pools of open water and scattered clumps of Arecaceae (palms) dominated by Mauritiella sp. A recent study of CV (Smith et al., 2021) demonstrates that palm swamps can yield Quaternary pollen archives recording terrestrial vegetation histories beyond the basin margin. To complement the previously published pollen data, here we conducted geochemical, stable carbon (δ¹³C) and nitrogen (δ¹⁵N) isotope, and macrocharcoal analysis from the same CV sediment core.

The second study site is Laguna Chaplin (LCh, 14°28′12″S, 61°2′60″W; elevation ~170 m a.s.l.), a large (~12 km²), shallow (2–2.5 m in the dry season), flat-bottomed lake (within a ~20 km² basin) located ~6.5 km north-east of CV (Mayle et al., 2000; Burbridge et al., 2004). Today, LCh is surrounded by HETF, with a mix of seasonally-inundated riverine forest and terra firme (upland) HETF. Previous pollen and charcoal studies of LCh have documented evidence of climate-driven ecotonal dynamics that replaced flammable SFS with less flammable HETF vegetation, maize crop cultivation, and Indigenous fire use during the Holocene (Mayle et al., 2000; Burbridge et al., 2004; Maezumi et al., 2018b).

The third study site is Huanchaca Mesetta (HM, 14°32′10.66″S, 60°43′55.92″W; elevation: 800 m a.s.l.), a small palm swamp (0.02 km² area) located on the eastern portion of NKMNP atop the Huanchaca table mountain ~40 km east of CV (Figure 1). HM is comprised of a monospecific stand of Mauritia flexuosa palms. While previously published changes in vegetation composition and fuel loads have been documented in phytolith and charcoal data from within the HM palm swamp during the Holocene (Maezumi et al., 2015), the broader table mountain has likely been dominated by edaphically constrained cerrado savannah vegetation for millions of years as a result of the rocky substrate and lack of soil on the table mountain (Killeen and Schulenberg, 1998; Mayle and Whitney, 2012). To date, no archaeological or palaeoecological evidence of human occupation or fire use has been recorded on the table mountain (Maezumi et al., 2015). Based on limited archaeological research on the table mountain, we assume that the 14,500-year-old charcoal record from HM represents a natural, lightning-caused cerrado savannah fire regime (Maezumi et al., 2015). HM is thus used as a general estimate of natural lightning-caused ignitions in the NKMNP region, serving as a “control” site to explore the possible influence of human-caused ignitions at LCh and CV.

2.5 Sediment core

A sediment core from CV palm swamp was collected in August 1995 by FM using a Livingstone modified square-rod piston corer (Wright, 1967) ~300 m from the eastern edge of the palm swamp and described in detail in Smith et al. (2021). The top 20 cm of the core was comprised of a floating mat of sedges and grasses. The 0 cm mark for collecting the core was the surface of the floating mat of grass and sedge which was ~30 cm thick. The water-sediment interface was between 30 cm depth and 155 cm depth. Sediments were solid enough to recover between 155 cm and 309 cm depth. The CV sediment cores were photographed and lithological descriptions were based on Munsell soil color charts and the texture of the sediment core (Figure 3). Visual descriptions, including sediment type, structure, texture, and organic content were undertaken to assist interpretation of the palaeoenvironmental data and comparisons with LCh and HM.

2.6 Chronology

The chronological framework for CV is based on nine Accelerator Mass Spectrometry (AMS) radiocarbon (¹⁴C) dates originally published by Smith et al. (2021) (Figure 3; Supplementary Table 1). Because of the lack of sufficient plant macrofossils, the majority of the dates were obtained from non-calcareous bulk sediment. However, two of the samples (Beta-467884 and Beta-467885) contained enough decayed plant remains following pre-treatment to be dated. All samples selected for dating were treated to remove carbonates, and the plant remains were treated to remove mobile humic acids. The radiocarbon ages were calibrated using the IntCal20 calibration curve (Reimer et al., 2020), and an age-depth chronology was constructed using the Bayesian age modeling software Bacon v2.3.4 (Blaauw and Christen, 2011). As per previous studies (Whitney et al., 2011; Smith et al., 2021), the IntCal20 calibration curve was selected over SHCal20 because of the hydrological connection of the study area to the northern hemisphere, via the SASM (McCormac et al., 2004; Hogg et al., 2013).

2.7 Magnetic susceptibility

Magnetic susceptibility (MS) was measured to identify mineralogical variation in the sediments (Nowaczyk, 2001). The MS of sediments reflects the relative concentration of ferromagnetic (high positive MS), paramagnetic (low positive MS), and diamagnetic (weak negative MS) minerals or materials. Shifts in the magnetic signature of the sediment can be diagnostic of a disturbance event (Gedye et al., 2000). Sediment cores were scanned horizontally, end to end, at 1 cm intervals using a Barrington ring sensor equipped with a 75 mm aperture. Data were calibrated by removing background noise at the beginning and end of each core.

2.8 Loss on ignition

The variability in the organic content of the sediments was used, in conjunction with magnetic susceptibility measurements, to identify periods of variability in sediment composition and organic content throughout the sediment record. Loss-on-ignition (LOI) analysis was carried out at 4 cm intervals through the CV core. After drying at 100°C for 24 h, each 1 cm³ sample was combusted at 550°C for 2 h (LOI550). The relative loss of weight before and after combustion determines the percentage organic carbon content that was present in that sample (Dean Jr, 1974; Heiri et al., 2001).

2.9 X-ray fluorescence

X-ray fluorescence (XRF) analysis of the sediment core was carried out using the ITRAX core scanner at Aberystwyth University at a step size of 2,000 or 5,000 μm. An intense micro-X-ray beam focused through a flat capillary waveguide was used to irradiate samples to enable both X-radiography and X-ray fluorescence (XRF) analysis. Data were acquired incrementally by advancing a split core, via a programmable stepped motor drive, through the flat, rectangular-section X-ray beam (Croudace et al., 2006). Results were normalized using z-scores for comparability (Figure 4).

2.10 Stable isotope analysis

Stable carbon (δ¹³C) isotope analysis was performed as a proxy for changes in vegetation structure and composition. The carbon in soil organic matter (SOM) is mainly derived from plants, resulting in a δ¹³C value that closely reflects that of the standing plant biomass (Ambrose and Sikes, 1991; Bond et al., 1994) and is primarily determined by the photosynthetic pathway of vegetation (Ehleringer et al., 1991; Ehleringer and Cerling, 2002; Sage, 2004). Previous research on δ¹³C values from NKMNP and the surrounding region (Pessenda et al., 1998; Killeen et al., 2003; Mayle et al., 2007; Maezumi et al., 2015; Lombardo et al., 2019) have been used to determine the relative proportions of C₄ savannah grasses vs. C₃ woody and herbaceous vegetation based on differences in the discrimination of the C₄ and C₃ photosynthetic pathways against the heavier ¹³C isotope during CO₂ fixation (Smith and Epstein, 1971). Open vegetation characterized by the predominance of C₄ grasses has a typical δ¹³C signature between −16.0‰ and −10.0‰, whereas dense vegetation characterized by trees is usually dominated by the C₃ photosynthetic pathway with a δ¹³C signature between −26.0‰ and −20.0‰ (Cerling, 1984), although the exact values of these distinctions have been shown to vary for the Amazon Basin (Tejada et al., 2020).

The interpretation of δ¹³C measurements of bulk SOM from lake sediments can be complex, however. Firstly, bulk SOM δ¹³C can be offset from the plant material from which it derives by +1–3‰ as a product of microbial action (Ehleringer et al., 2000). Though this offset cannot mask the non-overlapping differences between C₃ and C₄ taxa (Roberts et al., 2013), the relative changes in δ¹³C values throughout a sequence are often more meaningful than absolute comparisons. Secondly, the terrestrial C₃ and C₄ distinction are potentially complicated by the presence of CAM plants, which can have δ¹³C values overlapping with, or between, both of these photosynthetic pathways (O'Leary, 1981). Thirdly, variation in environmental conditions (e.g., aridity, temperature) has been shown to influence the δ¹³C values of C₃ plants in particular (Ehleringer and Björkman, 1977; Sage and Kubien, 2007; Kohn, 2010). Finally, the bulk organic matter of lake sediments is a product not only of terrestrial vegetation in the catchment but also aquatic plants (Cloern et al., 2002) which can also follow C₃, C₄ or CAM pathways (Keeley and Sandquist, 1992). Additionally, carbon limitation tends to reduce carbon isotope discrimination resulting in higher δ¹³C values in algae and submerged macrophytes (Sharkey et al., 1985; Keeley and Sandquist, 1992; Street-Perrott et al., 2004), making the teasing apart of terrestrial and lacustrine responses to climate and human activity complex. Consequently, comparisons with other, independent, records of vegetation and sedimentological conditions are essential.

Stable nitrogen (δ¹⁵N) isotope analysis provides an index of changing dynamics in the N cycle, particularly denitrification, as well as an indicator of the N₂-fixation processes. Plants, such as aquatic macrophytes (e.g., Sagittaria), can increase available N through N₂-fixation (Arima and Yoshida, 1982). Sediment δ¹⁵N integrates a variety of nutrient cycling processes including the loss of inorganic N to the atmosphere through denitrification which results in the subsequent increases in δ¹⁵N (Robinson, 1991; McLauchlan et al., 2013). Denitrification requires anaerobic conditions coupled with abundant available carbon and nitrate (Seitzinger et al., 2006). Thus, wet, anoxic sediments tend to have higher values of δ¹⁵N. Environmental conditions that alternate between wet (anaerobic) to dry (aerobic) conditions also lead to an increase in δ¹⁵N values (Codron et al., 2005). During dry periods, denitrification is stopped because of an increase in available oxygen in sediments, and thus δ¹⁵N values decrease. If dry soils become hydrated, there is a preferential loss of ¹⁴N, increasing δ¹⁵N values (Jordan et al., 1982; Codron et al., 2005). δ¹⁵N values have also been associated with forest cover in the Amazon (Ometto et al., 2006) and in the anthropogenic enrichment of soils in Southern Africa (Aranibar et al., 2008).

Subsampling for stable isotope analysis was performed at 1.5-cm-resolution throughout the length of the sediment core (Figure 5). One cm³ of bulk sediment was dried, powdered, and treated with 0.5 molar hydrochloric acid to remove carbonates. One to 25 mg of the dried carbonate-free sediment was weighed into tin capsules depending on organic matter content. Following reaction with 100% phosphoric acid, gases evolved from the samples were analyzed for stable carbon and nitrogen isotopic composition using a Finnigan Delta dual inlet elemental analyzer at the SIRFER Lab at the University of Utah. ¹³C/¹²C and ¹⁵N/¹⁴N ratios were presented in delta (δ) notation, in per mil (⁰/₀₀) relative to the international scale standard used for carbon, Vienna Pee Dee Belemnite (VPDB; USGS40 = −26.24⁰/₀₀; USGS41 = +37.76⁰/₀₀) and the standard for nitrogen, atmospheric N_2AIR(air) (USGS40 = −4.5⁰/₀₀ and USGS40 = −4.5⁰/₀₀) (Codron et al., 2005). δ¹³C and δ¹⁵N values were compared against International Standards [PLRM-1 (δ¹³C_VPDB = +23.96‰; δ¹⁵N_AIR = +49.63‰); PLRM-2 (δ¹³C_VPDB = +28.18‰; δ¹⁵N_AIR = −4.56‰); and SLRM MT soil (δ¹³C_VPDB = +7.5‰; δ¹⁵N_AIR = +16.96‰)]. Replicate analysis of SLRM standards suggests that machine measurement error is c. ± 0.2‰ for δ¹³C_VPDB and ± 0.2‰ for δ¹⁵N_AIR. Overall measurement precision was studied through the measurement of repeat extracts from the SLRM standard.

2.11 Macrocharcoal

Sediment samples were analyzed for macrocharcoal pieces (>125 μm in minimum diameter) using a modified macroscopic sieving method (Whitlock and Larsen, 2001; Leys et al., 2015) to reconstruct the history of local and extra-local fires. Macrocharcoal was analyzed using 1 cm³ volume of sediment in contiguous 0.5 cm intervals for the length of the sediment core. Samples were treated with 10% potassium hydroxide in a hot water bath for 15 min. The residue was sieved through a 125 μm sieve. Macroscopic macrocharcoal was counted in a gridded petri dish at 40 × on a dissecting microscope. Macrocharcoal counts were converted to macrocharcoal concentrations (number of macrocharcoal particles cm⁻³) and macrocharcoal influx rates by dividing by the deposition time (year cm⁻¹). CHAR statistical software was used to decompose charcoal data to identify distinct charcoal peaks using standard methodology (Higuera et al., 2007, 2009, 2010). Charcoal peaks are interpreted as a fire episode (a period of increased burning). The time difference between peaks is reflected in the fire return interval (fire frequency) for every 1,000 years.

2.12 Existing pollen data

The CV core was previously analyzed for pollen by Smith et al. (2021). Here, the Species Richness Index (SRI) was calculated from the existing pollen data which has previously been used to estimate floristic richness and biodiversity (Heck et al., 1975). Species richness was calculated from the total terrestrial fossil pollen data from CV (Smith et al., 2021) and LCh (Burbridge et al., 2004) using rarefaction analysis (Hurlbert, 1971; Heck et al., 1975) in the Vegan Community Ecology Package in R (Oksanen et al., 2017; R Core Development Team, 2020). Rarefaction enables the estimation of pollen richness using standardized pollen counts (Birks and Line, 1992). Standardization is important as palynological richness depends in part on the pollen sum because, as the sum increases, more taxa can potentially be detected. The units of rarefaction analysis represent the estimated number of pollen taxa reduced to a certain pollen sum, usually the lowest sum of the group of samples being compared.

2.13 Laguna Chaplin and Huanchaca Mesetta sediment cores

For complete palaeoecological analyses of LCh and HM see Mayle et al. (2000), Burbridge et al. (2004), and Maezumi et al. (2015, 2018b). Details on the recalibrated age-models follow. LCh and HM were cored in 1998 and 1995 by FM, with the results of the palaeoecological analyses presented in subsequent publications (Mayle et al., 2000; Burbridge et al., 2004; Maezumi et al., 2015, 2018b). Here, we replot the LCh and HM macrocharcoal and stable carbon isotope data with updated Bayesian age-depth models (Blaauw et al., 2018). The chronological framework for LCh presented here is based on 14 Accelerator Mass Spectrometry (AMS) radiocarbon (¹⁴C) dates (Supplementary Table 2, Supplementary Figure 1) and six AMS radiocarbon dates for HM (Supplementary Table 3, Supplementary Figure 2). Updated age models implement the Bacon Bayesian age modeling software package using the IntCal20 calibration curve (Reimer et al., 2020). CV, LCh, and HM are plotted from 24,000 to 3,500 to allow direct comparison with the 24,000-year CV macrocharcoal and isotope record.

2.14 Regional archaeological activity

The absence of archaeological sites in NKMNP, and the broader region of eastern Santa Cruz, the largest department in Bolivia, is attributed to the lack of systematic archaeological research in the region (Capriles, 2023). To the eastern edge of the park, across the Brazilian boarder, there are over 50 archaeological sites documented <50 km from the KNMNP (Figure 1). Additionally, recent maximum entropy models estimate 85% likelihood of pre-Columbian earthworks in the surrounding area (Peripato et al., 2023). Thus, in the absence of archaeological excavations and surveys within NKMNP, we relied on existing, regional archaeological data (Figure 1). We infer that these data are representative of NKMNP as these dated archaeological sites come from analogous settings across the ARC-Ecotone including sites that are coupled with palaeoecological studies documenting crop cultivation and agroforestry (including maize cultivation, similar to LCh) and Indigenous burning (Carson, 2014; Brugger et al., 2016; Maezumi et al., 2022). Additionally, the archaeological ceramics present within an interfluve ~25 km northwest of our study sites are similar to ceramics documented outside of the park (Riester, 1981; Carson et al., 2014; Barreto et al., 2016; Maezumi et al., 2022).

Identifying change in regional human activity levels and population size over time is frequently done with a “dates-as-data” approach (Rick, 1987). Under this paradigm, the abundance of dated archaeological deposits is thought to correlate over time with regional human activity in the sense that more dated deposits (more dates) in a given time and place indicates more human activity. This logic holds because a change in activity levels and/or an increase in population corresponds (though imperfectly correlated), to a change in the abundance of datable deposits entering the archaeological record at a given place and time. Following this logic, we compiled a radiocarbon date database (Supplementary Table S1) comprising archaeological deposits across the Bolivian ARF-ecotone into the Llanos de Moxos savannah ranging from 19,500 BP to 3,500 BP. We compared the abundance of those dates to our palaeoecological proxies for CV, LCh, and HM. The aim was to evaluate if regional human activity corresponded to changes in the palaeoecological record (i.e., fire activity).

To accomplish this, we used two visualizations to summarize the information contained in the radiocarbon database. The first visualization is Summed Probability Density (SPDs). SPDs have frequently been used as an indirect proxy for identifying change over time in human activity and demographic patterns (Shennan et al., 2013; Bevan et al., 2017; Ramsey, 2017). They are produced by calibrating each date in a given sample and then summing the probability distributions to produce a single density curve. That curve represents the temporal distribution of events in the relevant database along with their combined chronological uncertainty (Shennan et al., 2013; Timpson et al., 2014; Downey et al., 2016; Goldberg et al., 2016; Zahid et al., 2016). While popular, SPDs have significant limitations. Crucially, they conflate temporal event density or count with dating uncertainties and should largely serve only as a visual aid with some significant caveats regarding biases, calibration artifacts, and their point-wise interpretation (Carleton and Groucutt, 2021; Price et al., 2021; Crema, 2022). Nevertheless, they have become a staple of regional archaeological syntheses and can be deployed to help visualize the overall distribution of anthropogenic deposits and their relevant dating uncertainties in relation to other proxies.

The second visualization used here is called a Radiocarbon-dated Event Count Ensemble (RECE). A RECE provides visual information along three dimensions: one is time, another is count, and the other is probability (Carleton, 2021; Carleton and Groucutt, 2021). The ensemble refers to a collection of randomly sampled, probable event count time series that can be produced given a sample of radiocarbon dates. The dates are first calibrated, and then the calibrated dates are sampled at random in proportion to their densities producing a probable set of event dates. Each time the set of calibrated radiocarbon densities is sampled, the sampled event dates are then binned at a given resolution into a pre-determined temporal grid in order to count the number of events (dates, and by extension dated deposits) that fall into each bin. This produces an event count time series at the given resolution. The process is then repeated (i.e., 10,000 times) to produce the ensemble. A 2d grid histogram, is created where one axis represents time corresponding to the selected bin resolution and the other axis represents event count. The ensemble is then coded on the 2d grid based on the probable count sequence for each bin. These relative frequencies are then visualized with a heat-map spectrum so that the final image (i.e., the plotted RECE) provides a visual of the most likely count sequences along with uncertainty about both count and time in a single image. In this way, RECEs can help to separate event abundance from chronological uncertainty, suffer less from calibration artifacts than SPDs, have a clearer interpretation for comparison with other proxies, and indicate both overall sample size (maximum counts in one dimension) and total uncertainties in two dimensions.

To produce these visualizations, a total of 418 AMS dates were compiled from 98 archaeological sites from the Bolivian lowlands (ca. 69°-58°W, 9°-16°S; Supplementary Table S3). A binning procedure was applied to account for sites that have multiple dates within a phase (Shennan et al., 2013; Timpson et al., 2014; Goldberg et al., 2016). Dates within sites were ordered and those occurring within 100 years of each other were grouped into bins and merged. Each bin had a maximum width of 200 years. SPDs were then built in the package “rCarbon” (Crema and Bevan, 2021) and a RECE plot was produced using the package “chronup” (Carleton, 2021) using the IntCal20 calibration curve (Reimer, 2020; Reimer et al., 2020). In response to recent critiques of the SPD method such as minimum sample size (n = 500) (Williams, 2012) and the conflation of chronological uncertainty (Carleton and Groucutt, 2021), we refrain from interpretations about past population levels, and rather use the SPD curve and RECE plots as indicators of human presence in the region at a coarse temporal resolution with no reference to rates of change or palaeodemographics.

3 Results

The CV Bacon age-depth model was derived from seven of the nine AMS-radiocarbon dates and has previously been published (Smith et al., 2021), although here we update the chronology using the IntCal20 calibration curve (Reimer et al., 2020). The dates at 240 cm and 276 cm were rejected based on the outlier identification in Bacon. Ages for the top 160–155 cm and bottom 292 to 309 cm of the core were based on interpolation and extrapolation. The lithology of CV, as described in Smith et al. (2021), was divided into three main stratigraphic sections (Figure 3): (1) the Last Glacial Period (309–255 cm; 24,000–11,000 cal year BP), (2) the early to middle Holocene (255–230 cm; 11,000–7,000 cal BP), and (3) the middle to late Holocene (230–155 cm; 7,000–3,700 cal year BP). These three stratigraphic sections provide the chronological framework for discussing the new MS, LOI, XRF, stable carbon and nitrogen isotope, and macrocharcoal results from CV.

3.1 Last glacial period (24,000–11,000 cal year BP; sediment hiatus 18,000–11,000 cal year BP)

Magnetic susceptibility (MS) values from CV ranged between 4 and 8 and <10% organic material was present, as indicated by the LOI values. The XRF data indicated high frequency, low amplitude variability in Si, Rb, Ca, AL, Fe, Ti, and K (Figure 4). These data suggest a higher energy, nutrient-poor environment that may have been associated with substantial riverine influence. The pollen during this period suggests abundant Poaceae (grass 40%−60%) and Cyperaceae (sedge, 5%−10%) presence, along with cold adapted taxa including Podocarpaceae (Podocarpus), Betulaceae (Alnus), and Aquifoliaceae (Ilex). The savannah indicator Curatella americana (Dilleniaceae) was present (<1%). Only a few grains of the aquatic/semi-aquatic taxa from Alismataceae (Sagittaria sp.) and Isoetaceae (Isoetes sp.) were recovered (Smith et al., 2021). Pollen data are summarized in Figure 6.

Total %N remains low (<0.1%) and δ¹⁵N ranges between −0.6 and 3.4‰. Total %C was low (<1%) and the δ¹³C values range between −16.0 to −8.0‰, indicating the highest contribution of C₄ vegetation in the CV record during this period (Figure 5). The δ¹³C values at LCh (−18.0) and HM (−18.0 to −16.0‰) also indicate a contribution of C₄ vegetation at this time. The C/N values at CV fluctuate between 14 to 22 (Figures 6, 7), suggesting a mixed contribution of aquatic (~10) and terrestrial (~20) organic sources The existing palynological data from CV and LCh (Figure 8) indicate open, seasonally flooded savannah (Burbridge et al., 2004; Smith et al., 2021). The lack of phytolith preservation at HM prior to 15,000 cal year BP limits vegetation reconstruction on the table mountain (Maezumi et al., 2015). The absence of phytolith preservation at HM was attributed to the limited vegetation growing within the small palm swamp basin. Nevertheless, the δ¹³C values (−19 to −15‰, Figures 7, 8) suggest C₄ grasses were growing on the table mountain. The sediment hiatus between ~18,000 and 11,000 cal year BP at CV is synchronous with low sedimentation rates (~0.02 mm year⁻¹) at LCh (Burbridge et al., 2004) and a sediment hiatus at HM from ~18,000 to 16,000 cal year BP (Maezumi et al., 2015). Despite the presence of C₄ fire-adapted vegetation, such as grasses, macrocharcoal values were low at all three sites during this period. Human activity in the region was minimal between 19,500 and 11,000 cal year BP based on the regional radiocarbon dated archaeological sites (Figure 9).

3.2 Early to middle Holocene (11,000–7,000 cal year BP)

The end of the sediment hiatus at CV ~11,000 cal year BP (Figure 3) was associated with decreased MS values from ca. 8 to 4 ~10,000 cal year BP (Figure 4), suggesting a shift to a lower energy depositional environment. LOI values remained <10% indicating a continuation of organic-poor sediments. The XRF data, particularly Si, exhibits greater amplitude variability after 10,000 cal year BP, reaching record high values (~3) between 10,000 and 8,000 cal year BP that correspond with a sand layer (Figure 4). These lithological changes were likely associated with changing fluvial dynamics and greater riverine influence that increased Si-rich sand deposition at CV. The sand layer corresponds to poor pollen preservation at CV between 8,750 and 7,000 cal year BP (Smith et al., 2021; Figure 8).

Total % N remained low (<0.1) and δ¹⁵N values increased slightly from the previous period (1.0 and 5.0‰). Total % C was low (<1%) and the δ¹³C values indicate a progressive replacement of C₄ vegetation with C₃ vegetation at all three sites toward the late glacial-Holocene transition (Figures 6, 7). Coupled with an increase in C/N values (~16‰), these data indicate increased terrestrial organic matter input at CV (Figures 5, 6). Fire activity NKMNP increases at HM ~ 12,000 to 9,000 cal year BP and at CV ~10,500 to 9,600 cal year BP, indicated by charcoal influx and fire frequency values (Figures 6, 9). Macrocharcoal values remain low at LCh (<1) throughout this period. Fire activity was low between ~9,000 and 6,500 cal year BP at CV and LCh, (<2) coupled with a sediment hiatus at HM from ~9,000 to 8,300 cal BP prior to the onset of the mid-Holocene Dry Period (MHDP; 6,000–4,000 cal year BP). The SPD and RECE values indicate increased human activity at this time (Figure 9, Supplementary Figure S3).

3.3 Middle Holocene (7,000–3,700 cal year BP)

Around 6,500 cal year BP there were abrupt declines in MS (6 to −2) and an increase in organic content (>90%) at CV (Figure 4). This was synchronous with declines in Si, Rb, Fe, Ti and K values. Ca and Al exhibit higher frequency and greater amplitude variability after ~6,000 cal year BP. Pollen data show the presence of a dry forest-savannah mosaic during this time period (Smith et al., 2021). The aquatic/semi-aquatic taxa Sagittaria and Isoetes increase, with Sagittaria becoming a large proportion of the total pollen sum (~20%) between ca. 6,000–5,500 cal year BP (Figure 8). The increase in the aquatic C₄ plant Sagittaria corresponds to a shift in δ¹³C values at CV (from −24.0 to −12.0). These data correspond with low C/N values (<12; Figures 5, 7), indicating a greater contribution of aquatic vs. terrestrial carbon input. The increase in sediment organic matter at CV during this period suggests a shift toward more eutrophic conditions indicating lower water levels (Smith et al., 2021).

Between 6,000 and 4,000 cal year BP the δ¹³C values at LCh also indicate an increase in C₄ vegetation (Figure 6) that was not associated with a significant increase in grass pollen (Figure 8), that did correspond to a 10% increase in the aquatic CAM plant, Isoetes which has higher δ¹³C values (−20 to −10‰) compared to woody C₃ vegetation (<-24‰).

The MHDP was associated with increased variability of fire between CV, LCh, and HM. There was an increase in fire activity at CV ~6,500–6,000 cal year BP indicated by higher charcoal influx and fire frequency values that were not documented at LCh despite similar SDTF vegetation (Figures 8, 9) and close proximity (~5 km) of the lakes. At both CV and LCh there was a negative correlation between pollen species richness and increased macrocharcoal (Figure 8). This corresponds with an increase in human activity ca. 6,500 and 5,000 cal year BP indicated by a rise in SPD and RECE values (Figure 9). Fire activity at HM was low to absent during this period, despite abundant available fuel indicated by the presence of ca. 80% woody vegetation documented in the phytolith data (Figure 8). Fire increased between 5,000 to 4,000 cal year BP at CV and LCh, followed by decreased fire at both sites after ~4,000 cal year BP.The CV record terminates at 3,750 cal year BP, prior to the establishment of modern palm swamp vegetation.

4 Discussion

4.1 Climate and ARF-ecotone vegetation dynamics during the last glacial period (24,000–11,000 cal year BP)

During the Last Glacial Maximum (LGM; ~21,000 cal year BP) atmospheric CO₂ levels were ~76 ppm lower (Monnin et al., 2001), temperatures were ~5°C cooler than pre-industrial values (Seltzer et al., 2002; Whitney et al., 2011; Metcalfe et al., 2014), and decreased precipitation, in response to lower evaporative moisture, resulted in drier than present conditions across tropical lowland South America (Burbridge et al., 2004; Whitney et al., 2011; Woillez et al., 2011; Novello et al., 2019; Orrison et al., 2022). C/N values (>12) and δ¹³C (−16 to −12‰; Figures 7, 8), combined with the low abundance of C₄ aquatic plants documented in the pollen records, suggest that C₄ terrestrial vegetation dominated the NKMNP landscape between the LGM and onset of the Holocene (~24,000 to ~11,500 cal year BP). This interpretation corroborates existing pollen data indicating glacial-age savannah vegetation at CV (Smith et al., 2021), savannah and semi-deciduous dry forest communities at LCh (Burbridge et al., 2004), and C₄ dominated landscapes at HM (Maezumi et al., 2015; Figure 6). Similar glacial vegetation patterns of C₄ dominated landscapes have been documented at sites south of NKMNP in the Pantanal at Laguna La Gaiba (LLG; ~500 km southeast) (Whitney et al., 2011; Metcalfe et al., 2014) and Jaraguá cave (~850 km southeast of NKMNP) (Novello et al., 2019), suggesting the dominance of C₄ savannah grasslands during the last glacial period in the region. These data are consistent with previous modeling studies that simulate a replacement of rainforest with savannah/seasonally dry forest in ecotonal areas during the LGM (Beerling and Mayle, 2006). During the late glacial-Holocene transition, stable carbon isotope data indicate a progressive increase in C₃ vegetation (Figure 6) that corresponds with a period of poor pollen preservation at CV, an increase of semi-deciduous tropical dry forest at the expense of savannah vegetation at LCh after ~13,000 cal year BP, and a decrease in C₄ grass phytoliths at HM palm swamp after ~12,000 cal year BP (Figure 8), though it is likely that C₄ cerrado savannah persisted on the broader table mountain as a result of its rocky substrate and lack of soil (Mayle and Whitney, 2012; Maezumi et al., 2015). Model simulations demonstrate late glacial climate change (i.e., precipitation), drove the competitive replacement of drought-adapted vegetation (e.g., savannah or deciduous/semideciduous dry forest) by rainforest in drought-sensitive ecotonal areas (Beerling and Mayle, 2006). Precipitation increased in tropical lowland South America after ~12,800 cal year BP, as indicated by more frequent flooding and rising lake levels at LLG (Whitney et al., 2011), and in the Andes, water levels increased at Salar de Uyuni (Servant-Vildary et al., 2001) and Lake Titicaca (Baker et al., 2001). Reconstructions from LLG indicate that deglacial temperatures in the lowland Bolivian Amazon rose by 4°C at 19,500 cal year BP and have likely remained relatively consistent until present (Whitney et al., 2011), supporting the interpretation that temperature was not the dominant driver of late glacial-Holocene vegetation change at NKMNP. It is important to note however, that temperature may not have been uniform across lowland Amazonia. For example, the presence of Podocarpaceae (Podocarpus sp.) at Lagoa da Curuça (0°46′S, 47°51′W) in the eastern Brazilian Amazon has been attributed to late glacial-Early Holocene cooling (Behling, 1996).

The decline in drought-adapted C₄ vegetation associated with regionally wetter climates suggests that increased precipitation was likely the dominant driver shaping C₃-C₄ ARF-Ecotone vegetation dynamics in NKMNP during the glacial-interglacial transition. This interpretation is supported by palaeoecological studies from other drought-sensitive ecotonal areas around the Amazon Basin (Absy, 1991; Behling and Hooghiemstra, 1999; Sifeddine et al., 2001; Mayle et al., 2004). Interestingly, this interpretation contrasts with a recent record from Jaraguá cave (ca. 320 km east of NKMNP) that identified CO₂, opposed to temperature or precipitation, as the dominant driver of glacial-interglacial vegetation dynamics in western Brazil (Novello et al., 2019), demonstrating considerable spatial variability and the potential influence of different proxy data to determine drivers of late glacial vegetation responses in the tropical lowlands.

4.2 Documenting fire in C₄ Savannah grasslands during the last glacial period

During the LGP, minimal charcoal was present at NKMNP, however, determining the cause this absence of macrocharcoal is complicated by deposition and preservation issues associated with grass charcoal. Charcoal morphometrics based on aspect ratio can successfully differentiate woody and grass fuels (Vachula et al., 2021), however, previous studies document the overrepresentation of woody charcoal in grass dominated systems (Leys et al., 2015). To complicate this issue, recent charcoal dispersal models indicate that the dispersal of wood versus grass charcoal is not uniform, with fine grass charcoal particles being transported further than woody charcoal particles of the same size class (i.e., 125 μm) (Vachula and Rehn, 2023). Thus, the grass-charcoal component present in charcoal records from grass-dominated and mixed wood-grass ecosystems, such as NKMNP, likely represent broader regional fire histories than woody charcoal records from wood-dominated ecosystems. The pollen data at CV and LCh indicate the continued presence of woody fuels on the landscape during the LGP, thus if fire was present, we would anticipate woody charcoal to be present, if not overrepresented. In the absence of significant charcoal at any of the sites, we infer minimal LGP fire activity, across NKMNP.

The cooler glacial climate conditions and decreased evaporative moisture (Baker et al., 2001; Novello et al., 2019) in turn decreased the incidence of lightning-caused ignitions (Ramos-Neto and Pivello, 2000). The dominance of C₄ grasses would have reduced fuel loads and decreased fuel connectivity (Cochrane and Ryan, 2009). Coupled with the lower incidence of lightning-caused fires, these factors likely contributed to decreased fire activity in NKMNP during the LGP. Previous research suggests that ecological disturbance by fire is a fundamental component of C₄ grassland ecosystems (Keeley and Rundel, 2005; Beerling and Osborn, 2006; Edwards et al., 2010) and demonstrates that fire was an essential feedback mechanism leading to C₄ grassland expansion during the late Miocene (8–3 MYA) (Edwards et al., 2010; Strömberg, 2011; Bouchenak-Khelladi et al., 2014; Kirschner and Hoorn, 2019). Additionally, studies from the southwestern African grasslands of Namibia demonstrate that intensified fire activity caused by aridification triggered the initial establishment and later dominance of C₄ plants (Hoetzel et al., 2013). Although the presence of fire activity may have been essential in the initial evolution and establishment of C₄ grassland ecosystems, the absence of fire activity in our study sites during the LGM suggests that fire was not necessary for the presence and maintenance of C₄ savannah grasslands in NKMNP for over 24,000 years. These data challenge recent model simulations that indicate fire and low CO₂ were responsible for the expansion of the savannah vegetation in the southern ecotone margins of the Amazon during the LGM (Sato et al., 2021). Alternatively, our data support previous interpretations that low LGM CO₂ values and drier conditions likely helped to maintain the C₄ dominated system in NKMNP (Burbridge et al., 2004).

4.3 The influence of basin size in recording fire and vegetation dynamics in the ARF-ecotone

The glacial-Holocene transition (ca. 13,500–10,000 cal year BP) was characterized by a shift toward wetter conditions after 12,800 cal year BP in tropical lowland South America, coupled with oscillating periods of dry and wet conditions (Baker et al., 2001; Whitney et al., 2011; Novello et al., 2017). Despite increased precipitation variability, stable carbon isotope data indicate a progressive replacement of C₄ savannah grasslands with C₃ woody vegetation at our study sites (Figure 8). The stable isotope data corroborate existing pollen records indicating increased woody vegetation at CV and LCh (Burbridge et al., 2004; Maezumi et al., 2015) and decreasing C₄ phytoliths within the HM palm swamp (Maezumi et al., 2015). There was an increase in macrocharcoal values at HM ~12,000 to 9,000 cal year BP and at CV between 10,500 to 9,600 cal year BP.

The transition from cool-dry glacial climates to warm-wet climates in the early Holocene likely increased the incidence of convective thunder storms and lightning-induced fires at HM. Empirical studies of fire activity in the Emas National Park, Brazil (~900 km southeast of NKMNP) indicate that the majority of modern fires in the cerrado savannahs occur during the wet season (October–April), with a significantly higher occurrence of lightning-caused fires in open vegetation (wet field or grassy savannahs) on flat plateau areas, similar to the table mountain in NKMNP (Ramos-Neto and Pivello, 2000). Over a four-year period (1995–1999), wet season lightning fires were frequent (n = 31 fires burning ~30,000 ha), patchy, and extinguished by rain, likely representing the natural fire regime in the region (Ramos-Neto and Pivello, 2000). These data support the interpretation that HM likely represents the natural lightning-induced fire regime since the LGM.

Below the table mountain at CV, wetter climate conditions after 12,800 cal year BP (Whitney et al., 2011) likely increased fuel accumulation, fuel connectivity, and fuel moisture. The later may have made it more difficult for lightning-ignited fires to spread. LCh does not record an increase in fire activity during this period (Figure 9), despite similar declines in C₄ to vegetation indicated by the decrease in δ¹³C (from −16 to −24‰) and C/N values (>12; Figures 6, 7). The synchronicity of fire activity at the smaller lake basins of CV and HM, but not the larger lake basin of LCh, suggests basin size may have influenced the detection of fire activity between the records. In most forested ecosystems, macrocharcoal deposition typically occurs within 1 km of the lake (Clark et al., 1998; Lynch et al., 2004; Oris et al., 2014). However, in open grasslands in the South African Savannah, Afrotropics and North American tall-grass prairies, macrocharcoal deposition typically occurs within 5 km of the lake (Duffin et al., 2008; Aleman et al., 2013; Leys et al., 2015), though recent charcoal dispersal models suggest fine grass charcoal particles may be transported further (Vachula and Rehn, 2023). The close proximity of CV to LCh (<6 km) suggests these study sites should have similar fire histories. However, additional palaeoecological studies indicate small lake basins (<5 km²) are more sensitive to detecting local-scale (1–10⁶ m²) fire activity, while large sedimentary basins (>50 km²) capture regional (>10¹⁰ m²) macrocharcoal signals (Power et al., 2010). The presence of fire activity recorded in the small basins of CV and HM, but not the large basin of LCh, suggests these small basins are likely more sensitive to detecting local-scale fire activity (i.e. slash and burn cultivation fires) in the ARF-Ecotone. This interpretation corroborates existing studies documenting the influence of basin size on fire detection sensitivity in lake sediment basins ~300 km north of NKMNP (Carson et al., 2014).

4.4 Fire and human linkages during the early to middle Holocene (11,000–3,700 cal year BP)

In addition to basin size, the increase in the spatial variability of fire after 10,500 cal year BP at CV and LCh raises the possibility that climate (i.e., lightning) was not the dominant cause of ignitions. Previous studies have shown that increased spatial variability of fire, despite similar vegetation and climate, may be attributed to increased human influence (Bush et al., 2007a; Mayle and Power, 2008; Marlon et al., 2013; Maezumi et al., 2018b). A large spike in macrocharcoal has been documented at numerous sites in the Amazonian lowlands when areas were first occupied by humans (Bush et al., 2007a, 2016). Recent studies demonstrate that the ARF-Ecotone was extensively occupied since the early Holocene (after 10,500 cal year BP) with thousands of human-managed artificial forest islands documented across the Llanos de Moxos plains (~350 km west of NKMP) (Lombardo et al., 2013, 2020; Capriles et al., 2019). Additionally, recent distribution models estimate between 10,000 and 24,000 sites remain to be discovered in the region (Peripato et al., 2023).

The ARF-Ecotone region was also a center for early Holocene plant domestication (Oliver, 2008; Clement et al., 2010, 2015; de Cristo-Araújo et al., 2013; Brugger et al., 2016; Hilbert et al., 2017; Kistler et al., 2018; Watling et al., 2018; Lombardo et al., 2020). Cultivating crops (i.e., maize, manioc, sweet potato, squash, arrow root) was often coupled with Indigenous burning (Brugger et al., 2016; Bush et al., 2016; Maezumi et al., 2018b, 2022; Iriarte et al., 2020). Maize and Indigenous burning has been well-documented in the region in the anthropic forest islands [~6,850 cal year BP; (Lombardo et al., 2020)], Lake Rogaguado [~6,500 cal year BP; (Brugger et al., 2016)], Lake Versalles [~5,700 cal year BP; (Maezumi et al., 2022), and the Monte Castelo shell-mound [~5,300 cal year BP (Hilbert et al., 2017; Furquim et al., 2021)]. Previous palaeoecological studies demonstrate that crop cultivation was not uniform across the landscape; while some Amazon lake records had limited evidence of cultivation and fire activity, neighboring lakes with similar basin characteristics (<50 km away) had histories of crop cultivation and fire activity for millennia (Burbridge et al., 2004; Bush et al., 2007a,b; Carson et al., 2014; Brugger et al., 2016; McMichael and Bush, 2019). This may account for the spatial variability of fire and presence of maize cultivation at LCh (Maezumi et al., 2018b) but not at CV.

Another compelling line of evidence documenting the prevalence of past Indigenous burning comes from Amazonian Dark Earth Soils (ADEs). ADEs are human made soils that have high charcoal concentrations and organic matter, higher pH values and greater concentrations of P, Ca, and Mg which enable them to maintain nutrient levels over hundreds of years (Lehmann et al., 2003; Glaser and Woods, 2004; Woods et al., 2009). The earliest ADEs are documented in the Upper Madeira ~6,000 years ago (Mora et al., 1991; Watling et al., 2018), at the Teotônio site (ca. 600 km NW of NKMNP) where crops such as Cucurbita (squash) and Phaseolus (beans) were grown (Watling et al., 2018). While the widespread formation of ADEs across the Amazon did not begin until after the termination of the CV record ~2,500 cal year BP (Iriarte et al., 2020), the incipient land use practices trace their origins to the early to middle Holocene.

Archaeological studies have also identified the importance of riverine bluff habitats for Indigenous settlement, ADE soil formation, crop cultivation, and migration corridors (Denevan, 1996, 2005; Heckenberger et al., 2008; McMichael et al., 2014). The close proximity of CV to the Rio Paraguá (Figure 1) would have made it an optimal location for Indigenous utilization of the naturally open, mosaic landscape in the early Holocene (Figure 10). At CV there was a significant increase in macrocharcoal values at CV ~6,500 to 6,000 cal year BP that was not documented at LCh despite similar savannah-SDTF mosaic vegetation at both sites. The spatial variability of fire coupled with the synchronous increase in human activity in the region indicated by the SPD and RECE values (Figure 9), may be attributed to increased human-caused ignitions at this time. Climate likely also promoted fire by decreasing fuel moisture and increasing flammability during the MHDP, that would have in turn promoted natural fires (Bowman et al., 2009) and made Indigenous burning a more effective tool for land management (Mayle and Power, 2008; Bowman et al., 2020). However, while it is possible that humans have been influencing fire activity in NKMNP since the early Holocene (~10,500 cal year BP), additional archaeological evidence is needed to support this interpretation.

4.5 Climate and vegetation linkages during the mid-Holocene dry period (6,000–4,000 cal year BP)

The MHDP was characterized by significantly drier-than-present climate conditions with longer, more extreme dry seasons across much of tropical South America (Baker et al., 2001; Wang et al., 2007; Mayle and Whitney, 2012; Cheng et al., 2013; Kanner et al., 2013; Bernal et al., 2016). The MHDP has been attributed to lower southern-hemispheric summer insolation levels driven by the precessional cycle of Earth's orbit (Berger and Loutre, 1991). This restricted the southerly migration of the ITCZ (Haug et al., 2001) and decreased the strength of the SASM (Cruz et al., 2009; Baker and Fritz, 2015; Novello et al., 2017). The CV and LCh pollen records (Burbridge et al., 2004; Smith et al., 2021) and the phytolith record from HM (Maezumi et al., 2015) indicate a savannah-SDTF mosaic dominated the NKMNP landscape during the MHDP. Around 7,000 cal year BP higher δ¹³C values at CV correspond to an increase in aquatic/semi-aquatic macrophytes Sagittaria (~30%−40%) and Isoetes (~10%), indicating a shift from a clear open lake to a shallower more eutrophic environment (Smith et al., 2021). The dry conditions of the MHDP likely lowered lake levels creating optimal conditions for Sagittaria which grows in shallow waters (15–30 cm deep) along lake shores and requires exposed mud for successful seedling establishment (Jepson, 1993). Higher levels of %N (3–5; Figure 5) indicates increased nitrogen fixation, likely from the increase in Sagittaria which has high rates of N₂ fixation (Arima and Yoshida, 1982) and lower δ¹⁵N values (−2‰). At LCh the mid-Holocene expansion of Isoetes was also attributed to shallower conditions that allowed this aquatic plant to spread from the lake shore-line across the lake (Burbridge et al., 2004). The increased abundance of these aquatic macrophytes at CV and LCh was coupled with consistent levels of Poaceae (~50%) at both sites, suggesting C₄ grass was not the dominant driver of the δ¹³C values during this period. Some aquatic plants, such as Sagittaria (Sage, 2004), predominantly utilize the C₄ photosynthetic pathway (δ¹³C −12 to −13‰) (Ehleringer et al., 1991), while Isoetes uses the CAM photosynthetic pathway which can have a similar δ¹³C value (range −20 to −10‰) to C₄ plants (O'Leary, 1988). Together, the increase in Sagittaria and Isoetes, as a result of lower lake levels, likely drove the increased δ¹³C values at CV and LCh during the MHDP (Figure 10).

Additional evidence of vegetation response to regional drying has been documented in the Llanos de Moxos, where seasonally flooded savannahs were replaced by dry cerrado savannah vegetation (Mayle et al., 2007; Lombardo et al., 2019). The longer and more severe dry seasons were associated with a decrease in seasonal flooding during the MHDP (Mayle and Power, 2008), increased available land for cultivation and enabled the expansion of C₃ tree populations, intolerant of waterlogging, into the low-lying plains of the Llanos de Moxos (Mayle et al., 2007; Lombardo et al., 2019).

4.6 Fire activity and pollen species richness in the middle Holocene (7,000–3,700 cal year BP)

After ca. 5,000 cal year BP the spatial variability of fire in NKMNP begins to decrease, as indicated by similar macrocharcoal signatures at CV and LCh. The similarity in the macrocharcoal records suggest climate was likely the dominant driver of the regional increase in fire at CV, however, as humans are present in the region surrounding NKMNP (Figure 9), human-caused ignitions should not be ruled out and future archaeological research within the park may help to resolve this issue.

Interestingly, there is a slight negative correlation between fire activity and pollen species richness (Figure 8) at CV and LCh during this time that is likely attributed to declines in non-fire tolerant rainforest taxa (e.g., Moraceae/Urticaceae) immediately following a period of increased fire activity. However, at both CV and LCh, pollen diversity levels return to, or exceed pre-fire levels, suggesting that fire did not have a long-lasting negative influence on pollen species richness. Ethnographic studies have shown that through agroforestry practices (Denevan, 2006; Levis et al., 2018), even small groups of hunter-gatherers can change the configuration of useful plant species and create more productive environments surrounding human settlements (Posey and Balée, 1989; Politis, 1996; Denevan, 2006; Levis et al., 2018). For example agroforestry and Indigenous burning can increase alpha (local) and beta (regional) biodiversity in anthropogenic forests (Balée, 2000; Junqueira et al., 2010; Neves, 2012; ter Steege et al., 2013; Clement et al., 2015; Lins et al., 2015; Arroyo-Kalin, 2018). These studies, combined with the potential influence of Indigenous burning throughout the Holocene, raise the intriguing possibility that the composition of modern flora at NKMNP developed as part of a co-evolutionary process between people and plants that started at the beginning of its occupation as proposed in other sites in the Amazon (Furquim et al., 2023).

Fire activity decreases at CV and LCh after 4,500 cal year BP, likely as a result of the progressive increase in precipitation associated with the end of the MHDP (Mayle et al., 2000; Baker et al., 2001; Wang et al., 2007; Mayle and Whitney, 2012; Cheng et al., 2013; Kanner et al., 2013; Bernal et al., 2016; Lombardo et al., 2018) in response to the strengthening of the SASM driven by gradually increasing insolation levels (Berger and Loutre, 1991; Cruz et al., 2005; Baker and Fritz, 2015). This resulted in local seasonal flooding at CV and the establishment of M. flexuosa consistent with regional records indicating the expansion of both riverine and terra firme tropical rainforest (Mayle et al., 2000; Carson et al., 2014; Brugger et al., 2016). Unfortunately, the CV record terminates at 3750 cal year BP prior to the establishment of the modern palm swamp and the apex of regional human activity in the region (Heckenberger and Neves, 2009; Maezumi et al., 2018b; de Souza et al., 2019; Arroyo-Kalin and Riris, 2021).

5 Conclusions

This study reconstructs vegetation and palaeofire activity in the Bolivian ARF-Ecotone for the past 24,000 years (i.e., from the LGM to the late Holocene). Our data indicate that C₄ terrestrial savannah vegetation dominated the CV, LCh, and HM catchments during the glacial-Holocene transition. Additionally, fire was likely not a requisite for the maintenance of C₄ savannah-woodland ecosystems during the late glacial Period. The influence of basin size on the detection of palaeofire activity compliments existing pollen studies indicating the greater sensitivity of small lake basins to detect local-scale fire activity, supporting our hypothesis. Increased spatial variability of fire coupled with continuous presence of human occupation and regional maize cultivation suggests that humans may have contributed to the increased heterogeneity of fire at NKMNP during the Holocene. The close proximity of CV to the Rio Paragua would have made it an ideal settlement location for Indigenous occupation (likely less so once it became a swamp after ~3,750 cal year BP); however, archaeological surveys and excavations in the NKMNP area are still needed. While pollen richness declined with peak fire activity, richness values recovered to equivalent or higher levels following fire events, suggesting that the increase in fire at CV and LCh during the Holocene did not have a lasting negative influence on savannah-SDTF vegetation at NKMNP. While NKMNP has long been celebrated for its “pristine” ecosystems, the data presented here add to the increasing body of evidence that Indigenous peoples may have had a long, multi-millennial history of managing Bolivian Amazon ecosystems through fire and agroforestry. Hence, it is increasingly likely that, through millennia of land use, humans may have helped shape patterns of biodiversity in NKMNP, that potentially has a strong legacy in today's forests. Incorporating insights from Indigenous traditional ecological knowledge may contribute useful insights in combining biodiversity conservation and maintenance of ecosystem services with sustainable land use.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

Conceptualization: SM, MP, KM, and FM. Methodology, resources, project administration, and funding acquisition: SM, MP, and FM. Fieldwork and sample collection: FM. Formal analysis, data curation, writing—original draft preparation, and visualization: SM. Writing—review and editing: SM, RS, MP, AK, PR, CC, KM, and FM. Supervision: MP and FM. All authors contributed to the article and approved the submitted version.

Funding

Funding for this research was provided to SM by Global Change and Sustainability Center, the Graduate Research Fellowship, the Don Currey Graduate Research Fellowship, and the PAGES Graduate Research Fellowship. The University of Leicester provided a start-up grant to FM (1995) to support the fieldwork. SM and PR would like to thank the Max Planck Society for funding.

Acknowledgments

We thank Mary McIntyre and Daniel Harris for their help in sample preparation and analysis. We thank Dr. Tim Killeen and the Museo de Historia Natural Noel Kempff Mercado, Santa Cruz, Bolivia for providing logistical support, and in particular Rene Guillen and local guides from the village of Florida (e.g., Juan Surubi) for assistance with accessing the site, coring the site, and plant identification. We also thank the two reviewers for their critical reading and suggested improvements to help clarify our manuscript.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling editor GS declared a past collaboration with the author(s).

The author(s) AB, PR declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Publisher's note

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

Supplementary material

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

References