# Eccentricity-driven glacial climate variability and its influence on speciation in the tropical Andes

Manuel Chevalier, William D Gosling, Henry Hooghiemstra, Olivier Cartapanis, Brian M Chase, Stefanie Kaboth-Bahr

► To cite this version:

Manuel Chevalier, William D Gosling, Henry Hooghiemstra, Olivier Cartapanis, Brian M Chase, et al. Eccentricity-driven glacial climate variability and its influence on speciation in the tropical Andes. Quaternary Science Advances, 2025, 18, pp.100278. 10.1016/j.qsa.2025.100278. hal-04991097

HAL Id: hal-04991097

https://hal.science/hal-04991097v1

Submitted on 14 Mar 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Eccentricity-driven glacial climate variability and its influence on speciation in the tropical Andes

Manuel Chevalier<sup>a, *, ID</sup>, William D. Gosling<sup>b, ID</sup>, Henry Hooghiemstra<sup>b, ID</sup>, Olivier Cartapanis<sup>c, ID</sup>, Brian M. Chase<sup>d, e, ID</sup>, Stefanie Kaboth-Bahr<sup>f, ID</sup>

## ABSTRACT

This study presents a detailed pollen-based climate reconstruction from Lake Fúquene, Colombia, spanning the past 284,000 years and revealing significant changes in temperature and hydroclimate modes of variability, particularly between glacial periods. Our findings highlight the central role of eccentricity in modulating the influence of high-latitude glacial climate variability on the northern Andes. During periods of high eccentricity (~70-260 ka), we conclude that the region was shielded from northern high-latitude glacial ice-sheet dynamics by an intensified Walker circulation. After ~70 ka, a weakening of this zonal circulation allowed high-latitude, millennial-scale glacial climate variability to penetrate deeper into the tropics, influencing the position of the ITCZ and the rainbelt. This shift in climate dynamics led to increased regional millennial-scale climate variability. A detrended correspondence analysis (DCA) on fossil pollen data highlights a contemporaneous compositional shift previously masked by glacial-interglacial vegetation changes. The observed ecological shift coincides temporally with predictions of elevated speciation rates in mechanistic simulations. Combined, these results suggest a long-lasting environmental impact from glacial millennial-scale climate dynamics, supporting the hypothesis that increased environmental dynamism during glacial periods may trigger the speciation pump in the Northern Andes more than glacial inceptions or terminations. These findings underscore the evolving nature of regional climate drivers and their complex interplay with ecology, offering new insights into the role of climate variability in shaping biodiversity in this hotspot.

## 1. Introduction

Tropical mountains are home to the most biodiverse ecosystems on the planet (Antonelli et al., 2018; Myers et al., 2000; Rahbek et al., 2019b), a distribution largely explained by diverse modern climates (e. g., temperature and humidity (Antonelli et al., 2018)) and environmental factors (e.g., soil properties and topographic ruggedness (Antonelli et al., 2018; Flantua et al., 2019; Hoorn et al., 2013)). Recently, Coelho et al. (2023) highlighted the equally critical role of the spatial distribution of these variables—whether continuous or scattered—in determining global biodiversity patterns. However, the processes by which these evolving factors may have triggered regional “speciation pumps” remain poorly understood. One central debate revolves around whether biodiversity in tropical mountains arises from long-term climatic stability—allowing low extinction rates and gradual species accumulation—or from dynamic processes, where spatiotemporal

climatic variability and environmental heterogeneity fuel speciation (Flantua et al., 2019; Rahbek et al., 2019b). Given the escalating threats to biodiversity hotspots from human activities worldwide, identifying the mechanisms and timescales that generate and/or sustain biodiversity is crucial for effective conservation strategies.

The Northern Andean biodiversity hotspot, where ancient and recent lineages coexist (Rahbek et al., 2019a), has been extensively studied for its ecological history. A modelling experiment reconstructing South American biodiversity patterns over the past 800,000 years could reproduce the emergence of new lineages during the Quaternary glacial cycles, suggesting that modern biodiversity patterns reflect past interactions between climatic and topographic factors (Rangel et al., 2018). Intriguingly, these simulations revealed a marked increase in speciation rates in the Northern Andes during the last glacial period compared to earlier ones. However, empirical palaeoecological records have not been able to confirm this, and no explicit mechanisms were

<footnote><sup>a</sup> Institute of Geosciences, Sect. Meteorology, Rheinische Friedrich-Wilhelms-Universität Bonn, Germany</footnote>

<footnote><sup>b</sup> Institute for Biodiversity & Ecosystem Dynamics, University of Amsterdam, Amsterdam, the Netherlands</footnote>

<footnote><sup>c</sup> CEREGE, Aix Marseille Université, CNRS, IRD, INRAE, Coll. France, Technopole Arbois, 13545, Aix en Provence Cedex 4, France</footnote>

<footnote><sup>d</sup> Institut des Sciences de l'Evolution-Montpellier (ISEM), University of Montpellier, Centre National de la Recherche Scientifique (CNRS), EPHE, IRD, 4: Department of Environmental, Montpellier, France</footnote>

<footnote><sup>e</sup> Geographical Science, University of Cape Town, South Lane, Upper Campus, 7701, Rondebosch, South Africa</footnote>

<footnote><sup>f</sup> Institute for Geological Sciences, Freie Universität Berlin, Berlin, Germany</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: mchevali@uni-bonn.de (M. Chevalier).</footnote>
identified to explain this timing. Additionally, contrasting dynamics
across the two other South American tropical biodiversity hotspots point
to region-specific drivers or varying sensitivities to shared factors
(Rangel et al., 2018). Understanding these heterogeneous dynamics
requires robust, empirical reconstructions of past vegetation and envi-
ronmental change, although few existing records provide the temporal
resolution and duration needed to detect abrupt transitions across
multiple glacial-interglacial cycles (Hooghiemstra et al., 2022b).

Situated in the Eastern Cordillera of the Colombian Andes (5°27'N, 73°46'W, 2540 masl), the Lake Fúquene composite pollen record stands out for its extraordinary length (~284,000 years, nearly three complete glacial-interglacial cycles) and a high temporal resolution (~60-year intervals on average) (Bogotá-A et al., 2011; Groot et al., 2011), making it comparable in resolution to the best South American speleothem records that span similar periods (Deininger et al., 2019). To empirically characterise regional climatic variability, we analysed 4688 pollen samples from the Fúquene record to reconstruct the climatic conditions under which these botanical assemblages occurred. By integrating pollen-derived palaeoclimate reconstructions with independent regional and global records of past climate change, we investigate how Northern Andean environments responded to internal and external climate drivers during the past three glacial-interglacial cycles and explore possible mechanisms driving species diversification during the last glacial period.

## 2. Data and methods

### 2.1. *The Fúquene pollen datasets: Chronologies and data curation*

The composite Fúquene pollen dataset comprises two independent records: Fuq-2 (Van Geel and Van der Hammen, 1973) and Fuq-9C (Bogotá-A et al., 2011; Groot et al., 2011). Fuq-9C is a composite record derived from parallel sediment cores Fúquene-9 and Fúquene-10, which were retrieved from the lake's centre. It spans 27,000 to 284,000 years BP, with a mean temporal resolution of ~18 samples per millen-nium (4586 samples in total). In contrast, Fuq-2 was retrieved from the lake shore and covers the past 45,000 years (102 samples). Despite differences in sampling locations, vegetation trends in the overlapping sections (27,000 to 45,000 years BP) are consistent (Bogotá-A et al., 2011), suggesting a similar spatial representativity. We combined these two datasets to create a continuous record, anchoring the long, high-resolution Fuq-9C sequence to the present.

Fuq-2 was radiocarbon-dated (Bogotá-A et al., 2011), providing a robust chronology. In contrast, establishing a chronology for Fuq-9C, which extends beyond the radiocarbon horizon, was a complex task that required assumptions (Groot et al., 2014). The current chronology was developed by tuning arboreal pollen percentages to orbital obliquity, reflecting a strong response of vegetation to glacial-interglacial cycles (Groot et al., 2011). The pollen-based climate estimates presented in this study align with regional and global palaeoclimate patterns, thus supporting this chronology. Nonetheless, some discrepancies persist, such as a possible ~5000-year delay in the MIS 4 to MIS 3 transition (see Results). Similarly, while Dansgaard-Oeschger (D-O) oscillations are a major feature of the pollen record (and the reconstructions, see Results), the temporal uncertainties prevent us from precisely connecting each D-O event to the Greenland ice core benchmark. Despite such mismatches, we consider the chronology sufficiently precise for the analyses and discussions presented in this study.

### 2.2. *The CREST climate reconstruction method*

Traditional pollen-based climate reconstruction methods, such as the
Modern Analogue Technique (Anderson et al., 1991; Guiot, 1990) and
Weighted-Average Partial-Least-Square (ter Braak and Juggins, 1993),
face two limitations that restrict their application in tropical regions
(Chevalier et al., 2020). First, they rely on extensive collections of

modern pollen samples as input data to infer pollen-climate relation-
ships, collections that are currently insufficiently developed in the tro-
pics to reconstruct long-term climate change (Chevalier et al., 2020).
Second, their accuracy is constrained by the common identification of
pollen grains to genus or family levels, as these methods estimate only
the optimal climate conditions for taxa without fully accounting for the
broader range of climate tolerances under which the associated plants
can persist. This limitation becomes particularly problematic when these
tolerances are broad or when taxa exhibit multimodal climate responses,
such as in the tropics (Birks et al., 2010; Chevalier et al., 2020). To
overcome these challenges, we used the CREST (Climate REconstruction
SofTware) method. This "indicator species" climate reconstruction
technique estimates past climate conditions by combining the modern
climate responses of individual pollen taxa observed in fossil records
(Chevalier, 2022; Chevalier et al., 2014). Because CREST does not
require that the past vegetation communities have modern analogues
(Chevalier et al., 2020), it can be used to produce reliable climate esti-
mates from periods with distinct boundary conditions (e.g. Chevalier
and Chase, 2016, 2015; Hui et al., 2023).

CREST uses presence-only data, i.e. modern observations of where plant species associated with specific pollen morphotypes grow, instead of modern pollen rain samples to model the pollen-climate relationships (Fig. 1). This approach is particularly well-suited for regions like the northern Andes, where presence-only data are abundant and provide a robust alternative to modern pollen rain databases (Chevalier, 2019). A key innovation of CREST lies in its focus on estimating the climate responses of all species contributing to a pollen taxon within the study region, thereby integrating biodiversity and plant climate tolerances into its statistical framework (Fig. 1). This methodological approach is particularly advantageous for studying tropical records, where high plant diversity and less resolved pollen identification keys present significant modelling challenges.

At the species level, climate responses are modelled as parametric
probability density functions (PDFs), which are estimated from the
climate values associated with each species' modern distribution
(Fig. 1a,b). Species PDFs are modelled as normal distributions for
temperature-like variables, while for moisture-like variables (e.g., pre-
cipitation or aridity), log-normal distributions are preferred to prevent
possible negative values (Fig. 1c). The species PDFs are then combined
into pollen taxon PDFs to meet the taxonomic resolution of the studied
fossil record. Although some species are more abundant today than
others and are more likely to inform about local climate in modern
settings, we cannot determine how such relative contributions changed
in the past in response to changing and environmental interaction dy-
namics. As a result, we chose to assign the same weight to all the species
to estimate the pollen PDFs. Ultimately, the shape of the pollen PDFs is
unconstrained, allowing them to reflect the diversity of their contrib-
uting species (Fig. 1d).

Once fitted, pollen PDFs become the foundational units for recon-
structing climate parameters from fossil assemblages. For each fossil
pollen sample, CREST generates a likelihood distribution by weighting
and multiplying the pollen PDFs, with weights derived from normalised
pollen percentages. The weights are calculated as follows: each pollen
taxon time series is normalised by its mean percentage when the taxon is
present in the record. An average presence corresponds to a weight of 1,
irrespective of the absolute percentage observed, ensuring that even
low-abundance taxa with significant climatic relevance contribute to the
reconstructions (Chevalier, 2022; Chevalier et al., 2014). This normal-
isation of the percentages is central to accounting for the signals ported
by rare taxa, which often have more restricted climate requirements and
are thus key climate indicators. CREST derives climate estimates (most
likely value) and their associated tolerance (e.g., 95 % of the likelihood
distribution) from these likelihood distributions.

A distinguishing strength of CREST is its incorporation of taxonomic
uncertainties directly into the reconstruction framework. Since these
uncertainties are propagated into the final reconstructions, the results
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Conceptual representation of CREST illustrated with artificial pollen data. (a) Modern distribution of a climate variable to reconstruct (e.g. temperature). (b) Occurrences of four species part of the same pollen morphotype exhibiting marked preferences for the lowest values of that climate (e.g., dark/cold values) across the study area. (c) Combination of the four species PDFs (in colour) and the resulting pollen taxon PDF (in black). The histogram in the inset panel represents the proportion of the modern climate space (white) occupied by at least one of the four species (black), highlighting the higher chances of observing the taxon at the lower end of the climate gradient. (d) Taxon PDF produced by combining the four species PDFs and the type of synthetic statistics that can be derived from it (e.g., optimum, mean, uncertainty range). While the shape of the species PDFs is constrained (normal distribution on (c)), taxon PDFs can be irregular, skewed (such as on (d)), and even multi-modal (See Figs. S2–S4).</figcaption></figure>

exhibit broader uncertainty ranges than those produced by traditional statistical methods. However, this broader uncertainty offers a more accurate representation of the epistemic uncertainty of pollen-based climate reconstructions, better reflecting some of the “known unknowns” defined by Jackson (2012)— i.e. factors that influence results but are complex to quantify directly. By addressing these complexities, CREST provides a robust statistical framework for reconstructing past climates with improved reliability, even in tropical biodiversity hotspots.

## 2.3. Calibration data and variable selection

As with any transfer function, defining an appropriate calibration dataset is crucial and involves balancing two key considerations: 1) ensuring that all relevant conditions to be reconstructed are included and 2) avoiding the inclusion of improbable climate conditions and vegetation types in large quantities (Chevalier et al., 2020). To achieve this balance, we used the global calibration dataset *gbif4crest*, which integrates gridded vegetation and climate data, to estimate the species PDFs (Chevalier, 2019). To ensure the most extensive observation distributions and state-of-the-art climate data, we updated the *gbif4crest* dataset (v3, Chevalier (2025)) using the same code as used by Chevalier (2019) (v1) and Chevalier et al. (2021a) (v2). We also improved the grid resolution from 15' to 5' to better capture the topography of the Andes. This updated dataset is available through the crestr package (v1.4.0 and newer; Chevalier (2022)).

As the *gbif4crest* dataset is global, we restricted it geographically to focus on the vegetation of the northern Andes (defined here between 82°W and 68°W and 5°S to 12.5°N) and altitudinally to elevations above 1200 m asl (Fig. 2A). This altitude threshold corresponds to the altitudinal range of the vegetation documented by the Lake Fúquene pollen record, which spans the transition from sub-Andean forests during peak interglacials to páramo vegetation during peak glacials (Bogotá-A et al., 2016; Groot et al., 2013). Some fine-scale regional topographic details are inevitably lost due to the dataset's gridding. Because mountains are vertically compressed, the selected threshold of 1200 m asl corresponds to all observations with a “true” elevation ranging from 1000 to 1500 m asl, which matches the region's lower modern limit of sub-Andean forests. Consequently, warm tropical lowland vegetation, which was never present at the elevation of Lake Fúquene, was excluded from the calibration data.

We selected Mean Annual Temperature (MAT), Mean Annual Precipitation (MAP), and Aridity Index (AI) as key climate variables to

reconstruct (Fig. 2B). AI, which is calculated as the ratio of precipitation to potential evapotranspiration (Trabucco et al., 2008), was square root transformed to account for its non-linearity (Turner et al., 2020). These three variables are crucial in defining the vegetation of the Northern Andes. MAT and MAP were obtained from the WorldClim2 dataset (Fick and Hijmans, 2017), while AI came from the Global Aridity Index and Potential Evapotranspiration (ET0) Climate Database v3 (Zomer et al., 2022). We opted for annual rather than seasonal variables because the modern climate at Lake Fúquene is not strongly seasonal, which limits our ability to reliability estimate pollen-climate relationships for seasonal indicators.

The vegetation changes inferred from the pollen record document altitudinal migrations of vegetation zones, supporting the possibility of reconstructing MAT. MAP and AI also show significant variation across the region covered by the calibration data. The relatively low correlation between MAT and MAP/AI (0.47/0.17) suggests that MAP/AI can be independently reconstructed from MAT using the pollen data (Fig. S1). However, AI and MAP are strongly correlated (0.91), indicating a strong relationship between these variables in modern settings. Still, studies have shown that the rainfall-AI relationship may change over time and that large temperature changes, especially glacial-interglacial shifts, can overprint rainfall changes, thus breaking this relationship (Baxter et al., 2023; Chevalier and Chase, 2016). Overall, the climate breadth of the selected calibration dataset is substantial, with over 20 °C for MAT, 3000 mm for MAP, and 1.2 for AI (Fig. 2C–E, and S1). This range is more extensive than the climate changes needed to convert sub-Andean to páramo vegetation and, thus, sufficient to analyse the Fúquene pollen record.

## 2.4. Taxa classification and selection

Pollen-based climate reconstructions with high signal-over-noise ratios are obtained when only the terrestrial pollen taxa responsive to the climate variables being reconstructed are considered. To meet this requirement, we first excluded aquatic plants, ferns, and local indicators (e.g., plants growing on lake shores) from both datasets. We then recalculated the percentages of the remaining regional terrestrial taxa accordingly. The Fuq-2 dataset contains 73 distinct terrestrial taxa. However, 14 of these are observed in only one sample, and about half never exceed 1 % of the pollen sum.

In contrast, the Fuq-9C dataset includes 50 taxa due to its targeted counting strategy, which prioritised identifying key regional taxa to allow for more samples (Bogotá-A et al., 2011). Most Fuq-9C taxa are
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Calibration data and climate reconstructions. (A) Map describing the spatial and density distributions of the modern vegetation observations used to calibrate the pollen-climate relationships, with the climate space corresponding to these observations in (B). (C) MAT distribution in the selected climate space (orange). The orange circles indicate the climate preferences of the pollen taxa observed in the fossil record, and the diamond represents the modern climate at the site. (D) Same as (C) for MAP. (E) Same as (C) for AI.</figcaption></figure>

present in more than 100 samples, and they all exhibit maximum percentages above one. These differences in rare taxa between the two records likely reflect variations in counting strategies (exhaustive for Fuq-2 versus targeted for Fuq-9C) and potentially distinct taphonomic processes in their depositional settings. To account for these differences, we applied the pollen percentage normalisation procedure of CREST to each dataset individually (see CREST section). After harmonising taxonomic names across the two datasets, the combined Fúquene pollen record comprises 4688 samples and 81 unique terrestrial pollen taxa (See details in Table S1).

As an indicator species method, CREST evaluates each species independently, allowing different sets of climate-sensitive taxa to be selected for each variable (Chevalier et al., 2021b). This flexibility, however, introduces a challenge: rare taxa, represented by only one or

two pollen grains, which may have been transported over long distances, can disproportionately influence reconstructions. Since such taxa are not generally representative of the regional vegetation, we excluded all that do not constitute more than 1 % of the terrestrial pollen in any of the samples. This criterion excluded 28 taxa, impacted approximately 1.5 % of the samples (mostly from Fuq-2), and corresponded on average to about 1.7 % of the total terrestrial vegetation. Excluding such rare taxa reduced the sample-to-sample variability caused by taphonomy and counting uncertainties, thus increasing the reconstructions' reliability.

The *crestr* R package requires creating a proxy-species equivalency (PSE) table to assign species from the calibration dataset to pollen taxa (Chevalier, 2022). We adopted a conservative classification approach for most taxa, assigning all species within the calibration dataset corresponding to the pollen name. However, for common, hyperdiverse
taxa that occur under a wide range of environments, such as Asteraceae and Poaceae (Chevalier et al., 2021b), we refined their taxonomic resolution to the most common páramo genera (24 and 9 genera, respectively), as their occurrences in the record strongly correlated with other identifiable páramo vegetation taxa (see Table S1). To ensure robust species PDFs were estimated, we only included species with at least 20 distinct presence-only occurrences across the calibration data (Fig. 2). This conservative threshold is required to estimate the two parameters that define the species PDFs.

As a result, we could not estimate the climate response of ten taxa, two that were already excluded at the previous step (Didymopanax-type and Panopsis), three that were close to being excluded with a maximum percentage of 1.2 (Lysipomia, Daphnopsis, and Styloceras) and five potentially more important indicators (Draba, Eugenia, Thalictrum, Symplocos, and Lycopodium foveolate). However, the limited data available for the latter can have two different causes: 1) the taxa are rare, which would, in any case, prevent estimating reliable PDFs for them, or 2) the taxa mainly exist beyond the defined calibration dataset (i.e. further south of the Andes, or in the tropical lowlands), which means they are unlikely to have significantly contributed to the vegetation surrounding Lake Fúquene. With all these parameters and taxa classification, we estimated individual species PDFs for 493 species from 20,364 unique plant occurrence data across the study area, covering the three selected climate variables (Fig. 2). All these data were retrieved from the gbif4crest v3 database (Chevalier, 2018, 2019, 2025) that compiled vegetation observation datasets (GBIF.Org, 2024a, 2024b; 2024c; 2024d; 2024e; 2024f; 2024g; 2024h; 2024i), all sourced from the Global Biodiversity Information Facility (www.gbif.org; last access: August 25, 2024).

Ultimately, 51 taxa were retained for climate reconstructions, each evaluated for their suitability based on combined autecological and statistical assessments. Initially, all taxa were considered suitable climate indicators, and individual suitability evaluations followed. We used the breadth of the taxon's climate range as a key criterion for exclusion. Taxa with broad ecological ranges (e.g. mimosoid clade) or weak sensitivity to specific climate variables were excluded to avoid generating uninformative PDFs (Chevalier et al., 2021b). In contrast, the taxa for which our exclusion criteria were undecisive were kept as climate indicators. Exclusion criteria include:

* Amaranthaceae and Cactaceae were retained only for reconstructing moisture-related variables (MAP and AI) due to their strong association with seasonally dry soils and weak temperature sensitivity.

* Taxa with broad climatic distributions, such as Brassicaceae, Ericaceae, Hedyosmum, Melastomataceae, Miconia, and Morella, were excluded from MAT reconstructions and sometimes MAP/AI, depending on their observed responses.

* Quercus (Q. humboldtii) was excluded from all variables due to its broad climatic tolerance (Hooghiemstra et al., 2022a).

* Alnus (A. acuminata) was excluded due to its high pollen production, transportability, and preservation rates, often resulting in false observations in many Andean regions. In Lake Fúquene, changing percentages of Alnus mainly relate to lake-level changes (i.e. the potential surface at lake shores to bear Alnus-dominated swamp forest) and are, therefore, only indirectly related to climate change.

Importantly, all the exclusions of regional terrestrial taxa described in this section did not require recalculating percentages, as all excluded taxa remain representative of regional vegetation, even if they are unreliable climate indicators or cannot be modelled due to insufficient calibration data. Table S1 provides an exhaustive summary of all these decisions.

# 3. Results

## 3.1. Assessing reconstruction reliability

After designing a calibration dataset encompassing all plausible past climate and vegetation conditions at Lake Fúquene during the past 284,000 years, we obtained three reconstructions that confirm the strong influence of glacial-interglacial cycles in the northern Andes (Bogotá-A et al., 2011; Flantua and Hooghiemstra, 2018; Groot et al., 2011; Hooghiemstra et al., 2022b) (Fig. 3). The modern MAT values fall within the reconstruction uncertainties of the top samples, indicating that our modelling decisions were coherent. In contrast, the reconstructed MAP and AI values appear too high relative to modern values. This is not, however, surprising, considering that our calibration dataset “lacked” low rainfall/dry conditions relative to the modern Fúquene climate (Fig. 2B), thus shifting the entire reconstructions towards rainier or more humid conditions. The average amplitude of the three glacial-interglacial cycles is ~6.4 °C in MAT, ~642 mm in MAP, and ~0.11 AI units between full glacial and peak interglacial intervals (Fig. 3, Table 1). The most prominent climate change occurred at the transition between MIS 6 and MIS 5e, during which MAT increased by ~9 °C, MAP by > 800 mm, and AI by 0.17. However, the absolute amplitude of the hydroclimate reconstructions may be an underestimation of the true range due to the missing low rainfall/dry conditions in the calibration dataset.

Previous temperature reconstructions from Fuq-9C suggested a much larger amplitude (9–10 °C) between glacial and interglacial periods (Groot et al., 2011), which we cannot explain with the same argument of a truncated calibration dataset (Fig. 2B). However, several factors may explain this difference: 1) the percentage of tree pollen used by Groot et al. (2011) may not be strictly linearly correlated with temperature across the extensive range of tree pollen variability observed at Fúquene (from <10 to >90 % of the total pollen sum); 2) their method converted tree percentage changes into temperature, conflating multiple climate aspects into a single variable; and 3) the spatial resolution of the gbif4crest calibration data compresses the temperature elevation gradient, with colder values near the mountaintops being underrepresented. To estimate how this might influence our results, we calculated the average MAT values of our calibration data for the modern vegetation zones classified as Andean Forest (~16.1 °C on average, with extreme values over 21.3 °C) and páramo (~7.6 °C on average, with extreme values as low as 4.6 °C; Table 1). These two vegetation types represent the boundaries for our pollen data, suggesting our potential reconstructed range is empirically restricted between these values. CREST being based on likelihood, estimating the extreme values is nearly impossible, and this is why our calibration data does not allow reconstructing a thermal amplitude larger than ~9–10 °C, except for the transition from MIS 6 to MIS 5e during which the pollen transition was the most significant (Fig. 3A). We assume that hydroclimate reconstructions are less affected by these limitations, as these variables are more heterogeneously distributed along the mountain slopes (e.g. some low or mid-elevations can be dry due to a rain shadow effect).

Even if Fuq-2 and Fuq-9 show coherent vegetation change patterns in their overlapping sections, we observe differences in the amplitude of the reconstructions obtained from Fuq-2 (the most recent, lower-resolution samples) and Fuq-9. The amplitude of MAT and MAP reconstructions seems lower in Fuq-2 samples, while it appears slightly larger for AI. While it may be possible that the last 25,000 years were less extreme than the previous ~260,000 years for reasons that remain to be determined, these differences are more likely to reflect the distinct environments from which the sediments were recovered. Fuq-2 was retrieved near the modern shore of the lake, and the pollen record may be more sensitive to local conditions, such as lake-level changes, which strongly influence water access, a parameter ultimately tied to AI (Bogotá-A et al., 2011). This could explain why AI appears to have a stronger imprint on Fuq-2 samples than MAT and MAP. Fuq-9C is less
<figure><img src="image_4.png" /><figcaption>Fig. 3. Climate reconstructions. (A) MAT reconstructions (black curve) with their 50 % uncertainty range (coloured envelope). The triangle on the left y-axis represents the modern climate at the site. (B) Same as (A) for MAP. (C) Same as (A) for AI.</figcaption></figure>

impacted by such local variations.

The robustness of the three reconstructions was further investigated using leave-one-out (LOO) analyses (see Figs. S8–S10). In LOO analyses, the reconstructions are repeated multiple times, each time excluding one taxon (Chevalier and Chase, 2015). The reconstructed values are then compared to the original reconstructions. The differences highlight the influence of the excluded taxon on the results. High leave-one-out values were individually investigated to detect taxa that might have an undue influence on the results, using autecological insights to exclude specific taxa from the reconstruction of some climate parameters. The reconstructions were deemed robust after several iterations to identify and exclude taxa that were not sensitive to the reconstructed climate variable. Notably, the LOO results highlight the non-linearity of CREST reconstructions relative to the input taxa. They show that, despite the normalisation of the pollen percentages, some taxa provide significantly more information for reconstructing MAT (e.g., the cold *Plantago rigida*, and Poaceae ( páramo ), and the warm *Acalypha* and *Alchornea*), MAP (e.

g., the rain-loving *Acalypha*, *Alchornea*, and *Spermacoce*, and the drier Proteaceae and *Dodonaea*), or AI (e.g., the dry indicators such as Proteaceae, *Dodonaea*, and Amaranthaceae and the humid indicators *Sapium*, *Spermacoce*) at Fúquene, while many others are more “passive”. The relative importance of each taxon depends on the target variable (see Figs. S8–S10).

3.2. *Dominating trends and variability changes*

The strong correlations between the three climate recon-
structions—MAT, MAP, and AI (cor(MAT-MAP) = 0.84, cor(MAT-AI) =
0.66, cor(MAP-AI) = 0.85)—suggest either 1) an accurate climatic
signal, where the three reconstructions are statistically independent but
the climate parameters covaried over time in response to global and
local climate drivers, or 2) a statistical bias, where the reconstructions
are not fully independent due to relationships between the modern
variables used for calibration, representing the same primary climate
<caption>**Table 1**</caption>

Summary statistics of key intervals. Mean climate values and the 5–95 % uncertainty intervals of the reconstructions were estimated for each period by generating 5000 climate histories with a Monte-Carlo approach that sampled the climate uncertainties of each sample. Each history was averaged to estimate the interval’s mean climate. Here, ‘modern’ was defined as the past 2000 years, full glacial conditions as the periods between 19 and 22, 134–138, and 244–245 ka, and full interglacial conditions as the intervals between 3 and 7, 129–131, and 202–210 ka. The observed modern values correspond to the climate value of the Fúquene grid cell in the gbif4crest v3 dataset. Similarly, the values for the vegetation types were estimated by mapping the global biome dataset (Olson et al., 2001) on our grid. Here, we defined Páramo vegetation as a combination of the ‘Santa Marta páramo’, ‘Cordillera Central páramo’, ‘Cordillera de Merida páramo’, and ‘Northern Andean páramo’ ecoregions, and the Andean Forest is a combination of the ‘Magdalena Valley montane forests’, ‘Cordillera Oriental montane forests’, ‘Eastern Cordillera real montane forests’, and ‘Santa Marta montane forests’ ecoregions.

<table><thead><tr><th></th><th>MAT (°C)</th><th>MAP (mm)</th><th>AI (no units)</th></tr></thead><tbody><tr><td>Observed modern (from calibration dataset)</td><td>14.20</td><td>1213</td><td>0.956</td></tr><tr><td>Reconstructed modern (8 samples &lt; 2ka)</td><td>12.92<br>[11.14–14.69]</td><td>1561<br>[1317–1841]</td><td>1.085<br>[1.000–1.177]</td></tr><tr><td>Average full glacial conditions</td><td>8.87<br>[8.47–9.28]</td><td>1262<br>[1217–1309]</td><td>1.064<br>[1.048–1.082]</td></tr><tr><td>Average full interglacial conditions</td><td>15.19<br>[14.21–15.63]</td><td>1902<br>[1825–1980]</td><td>1.180<br>[1.156–1.204]</td></tr><tr><td>Average difference glacial - interglacial</td><td>6.41<br>[5.81–7.02]</td><td>642<br>[553–733]</td><td>0.114<br>[0.085–0.144]</td></tr><tr><td>Páramo vegetation</td><td>7.60<br>[4.61–11.44]</td><td>1156<br>[750–1789]</td><td>1.070<br>[0.883–1.288]</td></tr><tr><td>Andean forest vegetation</td><td>16.12<br>[9.67–21.34]</td><td>1705<br>[834–2808]</td><td>1.157<br>[0.793–1.504]</td></tr></tbody></table>

variability. However, the absence of correlation between MAT and the hydrological variables in modern settings (Figs. S5–S6), coupled with the demonstrated ability of CREST to produce independent reconstructions of these three parameters in other tropical environments (Chevalier and Chase, 2015, 2016), supports the former interpretation that the reconstructions capture accurate climate signals that covary at these timescales.

To identify the primary modes of climate variability governing the region, we performed a Principal Component Analysis (PCA) using normalised (i.e. centred and scaled) reconstructions (Fig. 4). A PCA is particularly useful in identifying the dominant patterns of variation across the three climate parameters. Since we have three input variables, only the first two components were interpretable (only two degrees of freedom). The third component corresponds to residuals (i.e., statistical noise). To improve the signal-to-noise ratio and reduce the sample-to-sample variability, we first applied a Gaussian smoothing with a kernel of 60 years, which corresponds to the average temporal resolution of the record (Rehfeld et al., 2011), before running the PCA. This transformation had marginal effects on the overall results (Fig. S11), and it smoothed out some of the high-frequency noise.

The first principal component (PC1) explains 85.7 % of the variance and strongly correlates all three variables—MAT, MAP, and AI (Fig. 4A and B). This suggests that the primary climatic variability driving the region is a combined signal from temperature, precipitation, and aridity, with generally warmer, rainier, and wetter interglacial periods contrasted with cooler and drier glacial periods. Typical of most tropical systems, this result aligns with the upward migration of vegetation bands along the elevation gradient of the Northern Andes, as inferred from the pollen data (Bogotá-A et al., 2016; Flantua et al., 2019; Groot et al., 2011). Orbital forcing factors such as precession, obliquity, and eccentricity contribute to the signal (Fig. 4D and F). An interesting observation is the attenuated expression of the MIS 7e interglacial period, in contrast to MIS 7c which shows more similarities to MIS 1 and MIS 5e, although it doesn’t reach their levels (Fig. 4C). This attenuated expression of MIS 7e is also observed in pollen records at lakes Junín

(Schiferl et al., 2023) and Titicaca (Hanselman et al., 2011) in the Peruvian and Bolivian Andes. This suggests that a mild MIS 7e interglacial may be a regional feature of the tropical Andes, potentially influencing vegetation belts and the extent of tropical glaciers (Rodbell et al., 2022). Another key feature of PC1 is the much higher millennial-scale variability reconstructed for the last glacial period than the previous ones (Fig. 4C and D). This stark difference highlights the distinctiveness of the last glacial period, with large and rapid fluctuations in temperature, precipitation, and aridity.

The second principal component (PC2), explaining 11.8 % of the variance, contrasts MAT (negative loadings) and AI (positive loadings), with MAP contributing marginally (Fig. 4A). This component highlights the inverse relationship between temperature and aridity, where higher temperatures lead to increased evapotranspiration (Van Boxel et al., 2014), thus reducing moisture availability to plants. PC2 also shows a loadings reversal between the last glacial period, marked by negative loadings, and the two preceding glacial periods, which have positive loadings. This reversal is independent of the transition between Fuq-9 and Fuq-2, further emphasising the uniqueness of the last glacial period relative to the earlier ones. The correlation between PC2 and the lake-level record derived from aquatic pollen and XRF data from Lake Fúquene (Fig. S12) further supports the interpretation of PC2 as representing the temperature-driven modulation of aridity during this period (Bogotá-A et al., 2011; Vriend et al., 2012).

# 4. Discussion

## 4.1. Drivers of glacial tropical climate variability

The individual reconstructions and PC1 highlight significant climate changes in temperature and hydroclimate during the last glacial cycles (Figs. 3 and 4). While the imprint of glacial cycles at Lake Fúquene had already been established (Bogotá-A et al., 2011; Groot et al., 2011), we present the first quantified analysis of three key parameters governing tropical environments. Temperature changes in the tropics are often considered less important than changes in the hydrological cycle, partly due to the limited availability of temperature-sensitive proxies. However, we show that the impact of temperature variations on vegetation at c. 2500 m elevation is comparable to changes in precipitation and humidity, contributing similarly to PC1. PC2 shows that temperature also interacts with the aridity index (AI), further extending findings from African records that highlighted the importance of temperature when interpreting hydrological data during periods of large temperature changes (Baxter et al., 2023; Chevalier and Chase, 2016). Although these African and Fúquene records come from mid-to-high elevation sites, where temperature changes are more important relative to the lowlands, accounting for temperature remains crucial when studying tropical mountains.

One prominent feature of the three reconstructions captured by PC1 is the millennial-scale climate variability during the last glacial period. Despite uncertainties in the chronology, the reconstructed changes’ relative succession, duration, and amplitude are consistent with the strong imprint of Dansgaard-Oeschger (D-O) events. D-O events, originating in the northern high latitudes, are known to affect tropical rainfall systems by latitudinally shifting the Intertropical Convergence Zone (ITCZ) and the tropical rainbelt in response to Northern Hemisphere warming and cooling phases (Cheng et al., 2016). These events have been widely documented in numerous South American hydrological records (Cruz et al., 2005; Deplazes et al., 2013; Fritz et al., 2010; Kanner et al., 2012), and here we show that they are also associated with substantial temperature changes (~3 °C on average), equivalent to an elevational shift of ~500–600 m in vegetation zones at Fúquene. In addition, these changes occur rapidly within a few pollen samples, which, thanks to the high temporal resolution of the Fúquene pollen record, indicates that vegetation can rapidly respond to abrupt climate forcing over decades.
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Principal Component Analysis (PCA) and wavelet analysis. (A) Distribution of the climate reconstructions in the PCA space, where the first (x-axis) axis represents 85.7 % of the total variance and the second (y-axis) 11.8 %. (B) Eigenvalues of the three components. Most of the variance is explained by PC1. (C) PC1 and (E) PC2 time series highlight our reconstructions' main trends. The curves are colour-coded by the value of the global $\partial^{18}$O marine stack (LR04 (Lisiecki and Raymo, 2005)). Wavelet analysis for (D) PC1 and (F), where the white dashed lines indicate the precession (23kyr), obliquity (41kyr), and eccentricity (100kyr) periodicities. The 1.5kyr period outlines the millennial-scale variability discussed in this paper (average duration of D-O cycles). The wavelet cone of influence (COI), delineated by the light grey shading, indicates the periodicities not impacted by edge effects.</figcaption></figure>

In contrast, earlier glacial periods (MIS 6 and MIS 8) show no similar high-amplitude variability, even though northern high-latitude variability persisted (Barker et al., 2011; Martrat et al., 2007; Rousseau et al., 2020) (Fig. 5a). Despite the decadal resolution of our record, we find no detectable changes in the pollen data or sediments that would suggest a loss of sensitivity to D-O events during these periods. This lack of rapid, high-amplitude variability in earlier glacial periods, coupled

with the millennial-scale changes observed during the last glacial period, leads us to propose two potential explanations. First, the intensity of D-O events during MIS 6 and MIS 8 may not have been strong enough to propagate to the latitudes of Lake Fúquene (~5°N). However, long-term D-O records indicate that, while the amplitude of these events was possibly smaller, it was unlikely small enough to mask their impact entirely. Second, the Northern Andes may have been shielded from these
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Regional and global comparisons. (f) PC1 from this study compared to (a) the $\delta^{18}$O record from NGRIP (North Greenland Ice Core Project members, 2004) (black) and the emulated Greenland $\delta^{18}$O record (Barker et al., 2011) (grey), the eccentricity parameter (b) and June-July-August insolation at 5N (c) (Berger and Loutre, 1991; Laskar et al., 2004), (d) the Chinese speleothem $\delta^{18}$O stack (Cheng et al., 2016), and (e) the Cariaco Basin Mo record from (Deplazes et al., 2013). The background is colour-coded by eccentricity changes to highlight low (green) and high (blue) eccentricity periods.</figcaption></figure>

events by an intensified Walker circulation during periods of higher eccentricity (~70–260 ka), possibly as a response to increased SST seasonality at high latitudes (Beaufort and Sarr, 2024). If correct, this means that tropical climates were more sensitive to low-latitude and local factors before ~70 ka, and this enhanced sensitivity blocked influences of high-latitude ice-sheet dynamics. We hypothesise that this shift was facilitated by a decrease in eccentricity, primarily linked to the low phase of the 400-kyr cycle (Fig. 5b). The decrease in eccentricity led to a weakening of the zonal circulation, ultimately allowing high-latitude climate signals to penetrate deeper into the tropics, influencing and modulating the position of the ITCZ and the tropical

rainbelt.

This absence of D-O-like variability during MIS 6 and MIS 8 is also reflected in other high-resolution tropical climate records, with the nature of the replacement varying depending on regional settings. In the Cariaco Basin, for example, D-O events from the last glacial period are replaced by lower-amplitude, longer-duration events (Deplazes et al., 2013) (Fig. 5e). Similarly, speleothem records from southern South America exhibit pure precessional signals before 70 ka (Cheng et al., 2013; Cruz et al., 2005). Similar shifts are seen in other regions, such as China (Cheng et al., 2016)(Fig. 5d), Ethiopia (Trauth, 2021), and South Africa (Chase et al., 2021; Partridge et al., 1997), and the global $\delta^{18}$O$_{atm}$
record tracking monsoon strength (Landais et al., 2010). These records suggest a global shift in tropical climate dynamics around 70,000 years ago, from locally-driven variability during high-eccentricity phases to increased influence from high-latitude forcing.

The loadings reversal of the “warmer-drier” tendency captured in the second principal component (PC2), which modulates the “warmer-rainier” signal of PC1, also shifts around 70 ka (Fig. 4E). Before this transition, colder and more humid conditions characterise glacial periods, while relatively warmer and drier conditions characterise the last glacial period (after ~70 ka). At the multi-millennial timescale, we observe a relationship between PC2 and the sea surface temperature (SST) gradient across the equatorial Pacific (Fig. S12), which we use here as a proxy for the strength of the Walker Circulation. This relationship evolves from a positive correlation during MIS 6 and MIS 8 to a negative correlation during MIS 2 to MIS 4. A steeper Pacific SST gradient is associated with increased convection over northwestern South America, where Martínez et al. (2003) showed that the high SST gradient during the LGM increased precipitation over Peru. However, this moisture was blocked by the Eastern Cordillera, creating dry conditions at Lake Fúquene. We hypothesise that deeper convective processes in the Neotropical atmosphere during higher eccentricity phases allowed Pacific moisture to cross the Western Cordillera's altitudinal barrier, leading to wetter conditions at Lake Fúquene during MIS 6 and MIS 8 and a positive correlation with the strength of the Pacific Walker cell.

Interestingly, the recently published speleothem $\delta^{18}$O record from the nearby Carracos Cave deviates significantly from our reconstructions (and most Northern Hemisphere tropical records) when interpreted through the standard framework for the tropics, where more depleted $\delta^{18}$O values are typically linked to increased rainfall (Ramirez et al., 2023). This leads, for instance, to infer elevated rainfall during the LGM, while our vegetation record suggests that conditions were much drier. However, this record is consistent with our findings if interpreted as being source- rather than amount-driven, and thus dominated by the interplay between the region's three primary moisture sources—the Atlantic and Pacific Oceans and recycled water from the Amazon Basin (Ramirez et al., 2023). Specifically, the $\delta^{18}$O record is positively correlated with PC2 before 60–70 ka and negatively correlated since then, mirroring the pattern observed with the Pacific Ocean SST gradient (Fig. S12). This suggests that the speleothem record did not directly record rainfall amount but instead reflected the interplay between zonal and meridional atmospheric factors, which modulated the relative influence of these moisture sources, ultimately determining rainfall at Lake Fúquene.

In conclusion, the two main features of PC1—change in millennial-scale variability—and PC2—inversion of glacial loadings—coupled with a pantropical compilation of records with sufficient temporal coverage and resolution, converge to indicate a significant shift in climate dynamics and sensitivities around 70 ka. This shift marks a transition from a state dominated by low-latitude forcing to one where high-latitude events have a more substantial impact, determining the most significant fraction of the variability observed. However, more long and high-resolution records, like the Fúquene pollen record, and climate model simulations with sufficient spatial resolution to capture the Andean topography are necessary across the Neotropics to better understand these shifts and their interactions with global climate patterns.

## 4.2. Possible implications for regional eco-evo processes

The climate reconstructions presented above reveal substantial climate variability during the past 284,000 years, especially during the last glacial period, and suggest that climatic dynamics played a central role in shaping the environment of Lake Fúquene. In particular, the rapid and pronounced temperature and hydroclimate fluctuations during the last glacial period likely had significant implications for local vegetation

and ecosystem structure. To investigate how these climatic fluctuations influenced vegetation dynamics, we performed a detrended correspondence analysis (DCA) of the pollen data, a standard technique used to extract dominant environmental gradients and track shifts in vegetation composition over time (Legendre and Legendre, 2012). The first axis of the Detrended Correspondence Analysis (DCA1) captures the glacial-interglacial variability in vegetation at Lake Fúquene that is obvious from the pollen diagram and the climate reconstructions (Fig. 6A). DCA1 demonstrates that temperature and moisture changes, as expected from glacial and interglacial phases, are strongly correlated and govern shifts in vegetation types. This aligns with previous studies of tropical mountain ecosystems, where temperature and precipitation are the primary drivers of vegetation dynamics during glacial-interglacial cycles. DCA1 thus confirms the general climatic framework observed in this region.

The second axis of the DCA (DCA2), which explains almost as much variance as DCA1, reveals a more intriguing trend (Fig. 6B). Unlike DCA1, which follows predictable glacial-interglacial patterns, DCA2 uncovers a striking shift in the vegetation state around 60,000 years ago at the transition from MIS 4 to MIS 3 (the trend begins precisely at the reconstructed transition, but the exact timing of this transition is delayed by a few millennia due to chronological uncertainties), i.e., when the eccentricity-driven influence of the Northern Hemisphere high latitude driven millennial-scale variability on many tropical records becomes most substantial (Fig. 5). Although happening in the time interval where the two Fúquene records overlap (see Methods), the DCA performed on the Fuq-9 data only shows the same trends (Fig. S16). The clear separation between earlier and later periods in the DCA2 suggests that significant ecological processes occurred during this time, triggering a distinctive change in the regional vegetation.

This notable finding did not occur during a glacial termination, which corresponds to massive shifts in climate states, but rather in the middle of the last glacial, when variability was high, but the amplitude of the change was only about half what is observed during a termination. This suggests that the cumulative effect of mid-amplitude events could have had more lasting impacts than rare high-amplitude events. Interestingly, this shift observed in DCA2 is contemporaneous with the pulse of diversification simulated by Rangel et al. (2018) (Fig. 6C-D-E). While the exact mechanism behind this shift remains speculative, this finding highlights the complex and dynamic relationship between climate variability and ecological change in the Northern Andes. Still, the shift observed in DCA2 supports the hypothesis that the millennial-scale variability of the last glacial period may have brought additional dynamism to the region's environments and triggered the speciation pump. It is also important to note that the intermediate peak in cumulated novelty simulated between 230 and 260ka (Fig. 6E) is not reconstructed in DCA2. We posit that it is because the data suggests an amplitude of climate change associated with the MIS 8 – MIS 7e transition much lower than that used in the model.

The dynamic interplay between temperature, precipitation, and aridity is central to understanding vegetation changes. The repeating climate events would have created environments with fluctuating ecological pressures, driving the response observed in DCA2. The timing of this shift, as identified in the DCA, suggests that the changing sensitivity of this region to high amplitude millennial-scale climate change events (e.g. D-O events) during periods of low eccentricity created ideal conditions to boost the speciation pump, thereby promoting a new vegetation composition during MIS 3 and MIS 2. The rapidly evolving climate of the topographically complex Northern Andes created intricate habitats where populations would have been regularly connected and isolated following a “flickering connectivity”, a mechanism crucial for the páramo vegetation near the mountaintops (Flantua et al., 2019).

In addition, the orientation of the Andes across the Equator means that D-O events manifest differently north and south of the Equator, with increased rainfall in the northern Andes and reduced rainfall southward (e.g. Cheng et al., 2013; Kanner et al., 2012). Long-distance
<figure><img src="image_7.png" /><figcaption>**Fig. 6. Biodiversity dynamics.** (A) First axis of the Detrended Correspondence Analysis (DCA1), (B) Second axis of the DCA (DCA2), (C) and (D) Simulated cradle (blue) and grave (red) richness in the Northern Andean Biodiversity hotspot from Rangel et al. (2018). (E) Cumulated difference between the cradle and grave richness over time, i.e. the cumulated sum of cradles – graves (the difference between the (C) and (D)) from 284,000 years ago to the present. Values close to 0 mean an overall compensation of new species and extinctions, and the more positive values indicate an accumulation of new species over time, which is largest after 70 ka. The data for (C), (D), and (E) were digitised from Rangel et al. (2018).</figcaption></figure>

habitat-tracking along the Andes would have also favoured interspecific
hybridisation (Flantua and Hooghiemstra, 2018), at least near the nexus
of the transition between dry and wet D-O events around the Equator
(DO events are dry in the Andes at ~6°S, Cheng et al. (2013)). Spatially
distinct responses across the Andes to the same climate events of the last
glacial cycle may thus have contributed to creating new vegetation as-
semblages and possibly enriching one of the most diverse biodiversity
hotspots of the world.

The increased environmental dynamism during MIS3, as revealed by
the climate reconstructions and DCA analysis, likely played a crucial role
in increasing the region's biodiversity, potentially driving speciation
and contributing to the exceptional biodiversity observed in the
Northern Andes. Although further mechanistic studies are needed, and
possibly longer records documenting previous periods of intense
millennial-scale variability in the region, it seems that the variability in
temperature and moisture played a crucial role in creating ecological

conditions conducive to this shift. While the exact nature of the speci-
ation or ecosystem change remains elusive, we propose that the
increased environmental dynamism during MIS 3, and possibly other
periods with similar conditions, contributed to increasing the region's
biodiversity.

5. Conclusion

In this study, we present a detailed climate reconstruction for Lake
Fúquene over the past 284,000 years, revealing significant variability in
temperature and hydroclimate, especially during the last glacial period.
Our findings highlight the central role of eccentricity in modulating the
influence of high-latitude forcing on tropical climates. During periods of
high eccentricity (~70–260 ka), the tropics were shielded from the
northern high-latitude ice-sheet dynamics by an intensified Walker cir-
culation, which reduced the impact of millennial-scale variability
originating from the Northern Hemisphere high latitudes. However, after ~70 ka, a weakening of this zonal circulation allowed such climate signals to penetrate deeper into the tropics, particularly influencing the position of the ITCZ and the rainbelt. This shift in climate dynamics led to significant changes in the regional environment, with increased climate variability likely driving ecological change. While the primary mode of vegetation variability follows glacial-interglacial cycles, our analysis reveals a striking shift in vegetation composition around 60,000 years ago, suggesting that the DO events of the last glacial had lasting ecological consequences. This shift coincides with the pulse of diversification previously simulated by Rangel et al. (2018), which supports the hypothesis that the increased environmental dynamism during this period may have triggered the speciation pump in the Northern Andes. While the exact mechanisms remain speculative, our findings underscore 1) the evolving nature of regional climate drivers and 2) the complex relationship between climate variability and ecological change in this biodiversity hotspot. Finally, this study also demonstrates the need for long, continuous, and high-resolution palaeoecological records, such as the Fúquene record analysed here, for improving our understanding of past climate and environmental changes and for addressing key questions regarding the interactions between climate, ecology, and biodiversity in the tropics.

### CRediT authorship contribution statement

**Manuel Chevalier:** Writing – review & editing, Writing – original draft, Visualization, Validation, Project administration, Methodology, Investigation, Formal analysis, Data curation. **William D. Gosling:** Writing – review & editing, Validation, Investigation, Conceptualization. **Henry Hooghiemstra:** Writing – review & editing, Validation, Methodology, Data curation. **Olivier Cartapanis:** Writing – review & editing, Visualization, Methodology, Investigation, Formal analysis. **Brian M. Chase:** Writing – review & editing, Methodology, Investigation, Formal analysis. **Stefanie Kaboth-Bahr:** Writing – review & editing, Validation, Investigation.

### Data and code availability

All the scripts and data produced by this study are available on Zenodo (https://doi.org/10.5281/zenodo.14959694). The reconstructions have also been submitted to www.pangaea.de The CREST analyses were performed with the crestr R package (v1.4.0)(Chevalier, 2022), available from http://doi.org/10.5281/zenodo.14718209 (or https://github.com/mchevalier2/crestr for the latest version), calibrated using the gbif4crest v3 dataset (https://doi.org/10.6084/m9.figshare.6743207.v12).

### Funding

MC is supported by the German Federal Ministry of Education and Research (BMBF) with the Research for Sustainability initiative (FONA) through the PalMod Phase II (grant no. FKZ: 01LP1926D) and the Phase III (grant no. FKZ: 01LP2308B) projects. HH produced the data set of Fúquene-9C over 7 years with financial support from the Netherlands Foundation for the Advancement of Tropical Research (WOTRO), grants W75-168, WB84-552 and WB84-636; The Netherlands Organization for Scientific Research (NWO), grants H75-284/CCH, 854.00.007/ALW, and 2012/13248/ALW); the University of Amsterdam grants, European Union Programme of High-Level Scholarships for Latin America (AlBan) grant E04D033907CO; and the Hugo de Vries Foundation (Amsterdam). OC was supported by ANR NEHOS (PI: T. Devièse), ANR MARCARA (PI: E. Bard) and CYCLOCARB MSCA (PI: O. Cartapanis).

### Declaration of competing interest

The authors declare that they have no known competing financial

interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgements

HH thanks Gavesa Co. (Bogotá) for drilling the cores and Bartriz Chaparro of the Corporacion Auto'noma Regional (CAR) for permission to work in Lake Fúquene. Annemarie Philip is thanked for preparing almost 5000 pollen samples. Pollen analysis of the Fúquene-9C record was carried out by a team of trained palynologists consisting of Mirella Groot, Giovanni Bogota', Juan Carlos Berrio, Alejandra Betancourt, Catalina Giraldo, Natalia González, Lineth Contreras, Diana Ortega, and Yahir Valderrama.

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.qsa.2025.100278.

### References

Anderson, P.M., Bartlein, P.J., Brubaker, L.B., Gajewski, K., Ritchie, J.C., 1991. Vegetation-pollen-climate relationships for the arcto-boreal region of north America and Greenland. J. Biogeogr. 18, 565-582.

Antonelli, A., Kissling, W.D., Flantua, S.G.A., Bermúdez, M.A., Mulch, A., Muellner-Riehl, A.N., Kreft, H., Linder, H.P., Badgley, C., Fjeldså, J., Fritz, S.A., Rahbek, C., Herman, F., Hooghiemstra, H., Hoorn, C., 2018. Geological and climatic influences on mountain biodiversity. Nat. Geosci. 11, 718-725. https://doi.org/10.1038/s41561-018-0236-z.

Barker, S., Knorr, G., Edwards, R.L., Parrenin, F., Putnam, A.E., Skinner, L.C., Wolff, E., Ziegler, M., 2011. 800,000 Years of abrupt climate variability. Science 334.

Baxter, A.J., Verschuren, D., Peterse, F., Miralles, D.G., Martin-Jones, C.M., Maitiuerdi, A., Van Der Meeren, T., Van Daele, M., Lane, C.S., Haug, G.H., Olago, D. O., Sinninghe Damsté, J.S., 2023. Reversed holocene temperature-moisture relationship in the horn of Africa. Nature 620, 336-343. https://doi.org/10.1038/s41586-023-06272-5.

Beaufort, L., Sarr, A.-C., 2024. Eccentricity forcing on tropical ocean seasonality. Clim. Past 20, 1283-1301. https://doi.org/10.5194/cp-20-1283-2024.

Berger, A.L., Loutre, M.-F., 1991. Insolation values for the climate of the last 10 million years. Quat. Sci. Rev. 10, 297-317.

Birks, H.J.B., Heiri, O., Seppä, H., Bjune, A.E., 2010. Strengths and weaknesses of quantitative climate reconstructions based on Late-Quaternary biological proxies. Open Ecol. J. 3, 68-110. https://doi.org/10.2174/1874213001003020068.

Bogotá-A, R.G., Groot, M.H.M., Hooghiemstra, H., Lourens, L.J., Van der Linden, M., Berrio, J.C., 2011. Rapid climate change from north Andean Lake Fúquene pollen records driven by obliquity: implications for a basin-wide biostratigraphic zonation for the last 284 ka. Quat. Sci. Rev. 30, 3321-3337. https://doi.org/10.1016/j.quacirev.2011.08.003.

Bogotá-A, R.G., Hooghiemstra, H., Berrio, J.C., 2016. North Andean environmental and climatic change at orbital to submillennial time-scales: vegetation, water-levels and sedimentary regimes from Lake Fúquene between 284 and 130ka. Rev. Palaeobot. Palynol. 226, 91-107. https://doi.org/10.1016/j.revpalbo.2015.09.007.

Chase, B., Harris, C., De Wit, M.J., Kramers, J., Doel, S., Stankiewicz, J., 2021. South African speleothems reveal influence of high- and low-latitude forcing over the past 113.5 k.y. Geology 49, 1353-1357. https://doi.org/10.1130/G49323.1.

Cheng, H., Edwards, R.L., Sinha, A., Spötl, C., Yi, L., Chen, S., Kelly, M., Kathayat, G., Wang, X., Li, X., Kong, X., Wang, Y., Ning, Y., Zhang, H., 2016. The Asian monsoon over the past 640,000 years and ice age terminations. Nature 534, 640-646. https://doi.org/10.1038/nature18591.

Cheng, H., Sinha, A., Cruz, F.W., Wang, X., Edwards, R.L., d'Horta, F.M., Ribas, C.C., Vuille, M., Stott, L.D., Auler, A.S., 2013. Climate change patterns in Amazonia and biodiversity. Nat. Commun. 4, 1411. https://doi.org/10.1038/ncomms2415.

Chevalier, M., 2025. GBIF for CREST database for crestr - version 3. https://doi.org/10.6084/m9.figshare.6743207.v12.

Chevalier, M., 2022. crestr: an R package to perform probabilistic climate reconstructions from palaeoecological datasets. Clim. Past 18, 821-844. https://doi.org/10.5194/cp-18-821-2022.

Chevalier, M., 2019. Enabling possibilities to quantify past climate from fossil assemblages at a global scale. Global Planet. Change 175, 27-35. https://doi.org/10.1016/j.gloplacha.2019.01.016.

Chevalier, M., 2018. GBIF database for CREST [WWW Document]. https://doi.org/10.6084/m9.figshare.6743207.

Chevalier, M., Chase, B.M., 2016. Determining the drivers of long-term aridity variability: a southern African case study. J. Quat. Sci. 31, 143-151. https://doi.org/10.1002/jqs.2850.

Chevalier, M., Chase, B.M., 2015. Southeast African records reveal a coherent shift from high- to low-latitude forcing mechanisms along the east African margin across last glacial-interglacial transition. Quat. Sci. Rev. 125, 117-130, 32.
Chevalier, M., Chase, B.M., Quick, L.J., Dupont, L.M., Johnson, T.C., 2021a. Temperature change in subtropical southeastern Africa during the past 790,000 yr. Geology 49, 71–75. https://doi.org/10.1130/G47841.1.

Chevalier, M., Chase, B.M., Quick, L.J., Scott, L., 2021b. An atlas of southern African pollen types and their climatic affinities. In: Quaternary Vegetation Dynamics – the African Pollen Database. CRC Press, London, pp. 239–258. https://doi.org/10.1201/9781003162766-15.

Chevalier, M., Cheddadi, R., Chase, B.M., 2014. CREST (Climate REconstruction Software): a probability density function (PDF)-based quantitative climate reconstruction method. Clim. Past 10, 2081–2098. https://doi.org/10.5194/cp-10-2081-2014.

Chevalier, M., Davis, B.A.S., Heiri, O., Seppä, H., Chase, B.M., Gajewski, K., Lacourse, T., Telford, R.J., Finsinger, W., Guiot, J., Kühl, N., Maezumi, S.Y., Tipton, J.R., Carter, V.A., Brussel, T., Phelps, L.N., Dawson, A., Zanon, M., Vallé, F., Nolan, C., Mauri, A., de Vernal, A., Izumi, K., Holström, L., Marsicek, J., Goring, S.J., Sommer, P.S., Chaput, M., Kupriyanov, D., 2020. Pollen-based climate reconstruction techniques for late Quaternary studies. Earth Sci. Rev. 210, 103384. https://doi.org/10.1016/j.earscirev.2020.103384.

Coelho, M.T.P., Barreto, E., Rangel, T.F., Diniz-Filho, J.A.F., Wüest, R.O., Bach, W., Skeels, A., McFadden, I.R., Roberts, D.W., Pellissier, L., Zimmermann, N.E., Graham, C.H., 2023. The geography of climate and the global patterns of species diversity. Nature 622, 537–544. https://doi.org/10.1038/s41586-023-06577-5.

Cruz, F.W., Burns, S.J., Karmann, I., Sharp, W.D., Vuille, M., Cardoso, A.O., Ferrari, J.A., Silva Dias, P.L., Viana, O., 2005. Insolation-driven changes in atmospheric circulation over the past 116,000 years in subtropical Brazil. Nature 434, 63–66. https://doi.org/10.1038/nature03365.

Deininger, M., Ward, B.M., Novello, V.F., Cruz, F.W., 2019. Late Quaternary variations in the south American monsoon system as inferred by speleothems—new perspectives using the SISAL database. Quaternary 2, 6. https://doi.org/10.3390/quat2010006.

Deplazes, G., Lückge, A., Peterson, L.C., Timmermann, A., Hamann, Y., Hughen, K.A., Röhl, U., Laj, C., Cane, M.A., Sigman, D.M., Haug, G.H., 2013. Links between tropical rainfall and North Atlantic climate during the last glacial period. Nat. Geosci. 6, 213–217. https://doi.org/10.1038/ngeo1712.

Fick, S.E., Hijmans, R.J., 2017. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315. https://doi.org/10.1002/joc.5086.

Flantua, S.G., Hooghiemstra, H., 2018. Historical connectivity and mountain biodiversity. In: Mountains, Climate and Biodiversity. John Wiley & Sons, Hoboken, NJ, pp. 171–185.

Flantua, S.G.A., O'Dea, A., Onstein, R.E., Giraldo, C., Hooghiemstra, H., 2019. The flickering connectivity system of the north Andean páramos. J. Biogeogr. 46, 1808–1825. https://doi.org/10.1111/jbi.13617.

Fritz, S.C., Baker, P.A., Ekdahl, E., Seltzer, G.O., Stevens, L.R., 2010. Millennial-scale climate variability during the Last Glacial period in the tropical Andes. Quat. Sci. Rev. 29, 1017–1024. https://doi.org/10.1016/j.quascirev.2010.01.001.

GBIF.Org, 2024a. Occurrence data downloaded on 25 August 2024 (Part 1). https://doi.org/10.15468/DLZGMNQ9.

GBIF.Org, 2024b. Occurrence data downloaded on 25 August 2024 (Part 2). https://doi.org/10.15468/DL9YKWC.

GBIF.Org, 2024c. Occurrence data downloaded on 25 August 2024 (Part 3). https://doi.org/10.15468/DL.UK2XV6.

GBIF.Org, 2024d. Occurrence data downloaded on 25 August 2024 (Part 4). https://doi.org/10.15468/DLQWCS68.

GBIF.Org, 2024e. Occurrence data downloaded on 25 August 2024 (Part 5). https://doi.org/10.15468/DLQ8ZUHH.

GBIF.Org, 2024f. Occurrence data downloaded on 25 August 2024 (Part 6). https://doi.org/10.15468/DLNUQ5TN.

GBIF.Org, 2024g. Occurrence data downloaded on 25 August 2024 (Part 7). https://doi.org/10.15468/DLMPFC47.

GBIF.Org, 2024h. Occurrence data downloaded on 25 August 2024 (Part 8). https://doi.org/10.15468/DL68HXG.

GBIF.Org, 2024i. Occurrence data downloaded on 25 August 2024 (Part 9). https://doi.org/10.15468/DL7BVEJK.

Groot, M.H.M., Bogotá, R.G., Lourens, L.J., Hooghiemstra, H., Vriend, M., Berrio, J.C., Tuenter, E., Van Der Plicht, J., Van Geel, B., Ziegler, M., Weber, S.L., Betancourt, A., Coneras, L., Gaviria, S., Giraldo, C., González, N., Jansen, J.H.F., Konert, M., Ortega, D., Rangel, O., Sarmiento, G., Vandenberghe, J., Van Der Hammen, T., Van Der Linden, M., Westerhoff, W., 2011. Ultra-high resolution pollen record from the northern Andes reveals rapid shifts in montane climates within the last two glacial cycles. Clim. Past 7, 299–316. https://doi.org/10.5194/cp-7-299-2011.

Groot, M.H.M., Hooghiemstra, H., Berrio, J.C., Giraldo, C., 2013. North Andean environmental and climatic change at orbital to submillennial time-scales: vegetation, water levels and sedimentary regimes from Lake Fúquene 130–27ka. Rev. Palaeobot. Palynol. 197, 186–204. https://doi.org/10.1016/j.revpalbo.2013.04.005.

Groot, M.H.M., Van Der Plicht, J., Hooghiemstra, H., Lourens, L.J., Rowe, H.D., 2014. Age modelling for Pleistocene lake sediments: a comparison of methods from the Andean Fúquene Basin (Colombia) case study. Quat. Geochronol. 22, 144–154. https://doi.org/10.1016/j.qageo.2014.01.002.

Guiot, J., 1990. Methodology of the last climatic cycle reconstruction in France from pollen data. Palaeogeogr. Palaeoclimatol. Palaeoecol. 80, 49–69. https://doi.org/10.1016/0031-0182(90)90033-4.

Hanselman, J.A., Bush, M.B., Gosling, W.D., Collins, A., Knox, C., Baker, P.A., Fritz, S.C., 2011. A 370,000-year record of vegetation and fire history around Lake Titicaca (Bolivia/Peru). Palaeogeogr. Palaeoclimatol. Palaeoecol. 305, 201–214. https://doi.org/10.1016/j.palaeo.2011.03.002.

Hooghiemstra, H., Cleef, A.M., Flantua, S.G.A., 2022a. A paleoecological context to assess the development of oak forest in Colombia: a comment on Zorilla-Azcué, S., González-Rodríguez, A., Oyama, K., González, M.A., & Rodríguez-Correa, H., the DNA history of a lonely oak: Quercus humboldtii phylogeography in the Colombian Andes. Ecology and Evolution 2021 12, e8702. https://doi.org/10.1002/ece3.8702, 10.1002/ece3.7529. Ecology and Evolution.

Hooghiemstra, H., Pérez, G.S., Torres Torres, V., Berrio, J.-C., Lourens, L., Flantua, S.G. A., 2022b. 60 years of scientific deep drilling in Colombia: the north Andean guide to the Quaternary. Sci. Drill. 30, 1–15. https://doi.org/10.5194/sd-30-1-2022.

Hoorn, C., Mosbrugger, V., Mulch, A., Antonelli, A., 2013. Biodiversity from mountain building. Nat. Geosci. 6, 154. https://doi.org/10.1038/ngeo1742, 154.

Hui, Z., Wei, X., Xue, Z., Zhao, X., Chevalier, M., Lu, X., Zhang, J., Peng, T., Chen, Y., Chen, P., 2023. Middle Miocene evolution of East Asian summer monsoon precipitation in the northeast part of the Tibetan Plateau based on a quantitative analysis of palynological records. Palaeogeogr. Palaeoclimatol. Palaeoecol. 630, 111808. https://doi.org/10.1016/j.palaeo.2023.111808.

Jackson, S.T., 2012. Representation of flora and vegetation in Quaternary fossil assemblages: known and unknown knowns and unknowns. Quat. Sci. Rev. 49, 1–15. https://doi.org/10.1016/j.quascirev.2012.05.020.

Kanner, L.C., Burns, S.J., Cheng, H., Edwards, R.L., 2012. High-latitude forcing of the south American summer monsoon during the last glacial. Science 335, 570–573. https://doi.org/10.1126/science.1213397.

Landais, A., Dreyfus, G., Capron, E., Masson-Delmotte, V., Sanchez-Goni, M.F., Desprat, S., Hoffmann, G., Jouzel, J., Leuenberger, M., Johnsen, S., 2010. What drives the millennial and orbital variations of 618Oatm? Quat. Sci. Rev. 29, 235–246. https://doi.org/10.1016/j.quascirev.2009.07.005.

Laskar, J., Robutel, P., Joutel, F., Gastineau, M., Correia, A.C.M., Levrard, B., 2004. A long-term numerical solution for the insolation quantities of the Earth. Astron. Astrophys. 428, 261–285. https://doi.org/10.1051/0004-6361:20041335.

Legendre, P., Legendre, L.F.J., 2012. Numerical Ecology. Elsevier.

Lisiecki, L.E., Raymo, M.E., 2005. A Pliocene-Pleistocene stack of 57 globally distributed benthic δ 18 O records. Paleogeography 20, PA1003. https://doi.org/10.1029/2004PA001071.

Martínez, L., Keigwin, L., Barrows, T.T., Yokoyama, Y., Southon, J., 2003. La Nina-like conditions in the eastern equatorial Pacific and a stronger Choco jet in the northern Andes during the last glaciation. Paleogeography 18, 2002PA000877. https://doi.org/10.1029/2002PA000877.

Matrat, B., Grimalt, J.O., Shackleton, N.J., de Abreu, L., Hutterli, M.A., Stocker, T.F., 2007. Four climate cycles of recurring deep and surface water destabilizations on the Iberian margin. Science 317, 502–507.

Myers, N., Mittermeier, R.A., Mittermeier, C.G., da Fonseca, G.A.B., Kent, J., 2000. Biodiversity hotspots for conservation priorities. Nature 403, 853–858. https://doi.org/10.1038/35002501.

North Greenland Ice Core Project members, 2004. High-resolution record of Northern Hemisphere climate extending into the last interglacial period. Nature 431, 147–151. https://doi.org/10.1038/nature02805.

Olson, D.M., Dinerstein, E., Wikramanayake, E.D., Burgess, N.D., Powell, G.V., Underwood, E.C., D'amico, J.A., Itoua, I., Strand, H.E., Morrison, J.C., others, 2001. Terrestrial Ecoregions of the World: A new Map of Life on Earth: a new global map of terrestrial ecoregions provides an innovative tool for conserving biodiversity. Bioscience 51, 933–938. https://doi.org/10.1641/0006-3568(2001)051[0933:TEOTWAJ2.0.CO;2].

Partridge, T.C., DeMenocal, P.B., Lorentz, S.A., Paiker, M.J., Vogel, J.C., 1997. Orbital forcing of climate over South Africa: a 200,000-year rainfall record from the Pretoria Saltpan. Quat. Sci. Rev. 16, 1125–1133. https://doi.org/10.1016/S0277-3791(97)00005-X.

Rahbek, C., Borregaard, M.K., Antonelli, A., Colwell, R.K., Holt, B.G., Nogues-Bravo, D., Rasmussen, C.M.Ø., Richardson, K., Rosing, M.T., Whittaker, R.J., Fjeldså, J., 2019a. Building mountain biodiversity: geological and evolutionary processes. Science 365, 1114–1119. https://doi.org/10.1126/science.aax0151.

Rahbek, C., Borregaard, M.K., Colwell, R.K., Dalsgård, B., Holt, B.G., Morueta-Holme, N., Nogues-Bravo, D., Whittaker, R.J., Fjeldså, J., 2019b. Humboldt's enigma: what causes global patterns of mountain biodiversity? Science 365, 1108–1113. https://doi.org/10.1126/science.aax0149.

Ramirez, V.M., Cruz, F.W., Vuille, M., Novello, V.F., Strikis, N.M., Cheng, H., Zhang, H. W., Bernal, J.P., Du, W.J., Ampuero, A., Deininger, M., Chiessi, C.M., Tejedor, E., Campos, J.L., Ait Brahim, Y., Edwards, R.L., 2023. Summer insolation controlled movements of Intertropical Convergence Zone during last glacial cycle in northern South America. Commun Earth Environ 4, 495. https://doi.org/10.1038/s43247-023-01124-6.

Rangel, T.F., Edwards, N.R., Holden, P.B., Diniz-Filho, J.A.F., Gosling, W.D., Coelho, M. T.P., Cassemiro, F.A.S., Rahbek, C., Colwell, R.K., 2018. Modeling the ecology and evolution of biodiversity: biogeographical cradles, museums, and graves. Science 361, eaar5452. https://doi.org/10.1126/science.aar5452.

Rehfeld, K., Marwan, N., Heitzig, J., Kurths, J., 2011. Comparison of correlation analysis techniques for irregularly sampled time series. Nonlinear Process Geophys. 18, 389–404. https://doi.org/10.5194/npg-18-389-2011.

Rodbell, D.T., Hatfield, R.G., Abbott, M.B., Chen, C.Y., Woods, A., Stoner, J.S., McGee, D., Tapia, P.M., Bush, M., Valero-Garcés, B.L., Lehmann, S.B., Mark, S.Z., Weidhaas, N.C., Hillman, A.L., Larsen, D.J., Delgado, G., Katz, S.A., Solada, K.E., Morey, A.E., Finkenbinder, M., Valencia, B., Rozas-Davila, A., Wattrus, N., Colman, S.M., Bustamante, M.G., Kück, J., Pierdominici, S., 2022. 700,000 years of tropical Andean glaciation. Nature 607, 301–306. https://doi.org/10.1038/s41586-022-04873-0.

Rousseau, D.-D., Antoine, P., Boers, N., Lagroix, F., Ghil, M., Lomax, J., Fuchs, M., Debret, M., Hatté, C., Moine, O., Gauthier, C., Jordanova, D., Jordanova, N., 2020.
Dansgaard-Oeschger-like events of the penultimate climate cycle: the loess point of view. Clim. Past 16, 713-727. https://doi.org/10.5194/cp-16-713-2020.

Schiferl, J., Kingston, M., Åkesson, C.M., Valencia, B.G., Rozas-Davila, A., McGee, D., Woods, A., Chen, C.Y., Hatfield, R.G., Rodbell, D.T., Abbott, M.B., Bush, M.B., 2023. A neotropical perspective on the uniqueness of the Holocene among interglacials. Nat. Commun. 14, 7404. https://doi.org/10.1038/s41467-023-43231-0.

ter Braak, C.J.F., Juggins, S., 1993. Weighted averaging partial least squares regression (WA-PLS): an improved method for reconstructing environmental variables from species assemblages. Hydrobiologia 269-270, 485-502. https://doi.org/10.1007/BF00028046.

Trabucco, A., Zomer, R.J., Bossio, D.A., van Straaten, O., Verchot, L.V., 2008. Climate change mitigation through afforestation/reforestation: a global analysis of hydrologic impacts with four case studies. Agric. Ecosyst. Environ. 126, 81-97. https://doi.org/10.1016/j.agee.2008.01.015.

Trauth, M.H., 2021. Spectral analysis in Quaternary sciences. Quat. Sci. Rev. 270, 107157. https://doi.org/10.1016/j.quascirev.2021.107157.

Turner, M.G., Wei, D., Prentice, I.C., Harrison, S.P., 2020. The impact of methodological decisions on climate reconstructions using WA-PLS. Quat. Res. 1-16. https://doi.org/10.1017/qua.2020.44.

Van Boxel, J.H., González-Carranza, Z., Hooghiemstra, H., Bierkens, M., Vélez, M.I., 2014. Reconstructing past precipitation from lake levels and inverse modelling for Andean Lake La Cocha. J. Paleolimnol. 51, 63-77. https://doi.org/10.1007/s10933-013-9755-1.

Van Geel, B., Van der Hammen, T., 1973. Upper quaternary vegetational and climatic sequence of the fuquene area (Eastern Cordillera, Colombia). Palaeogeogr. Palaeoclimatol. Palaeoecol. 14, 9-92. https://doi.org/10.1016/0031-0182(73)90064-3.

Vriend, M., Groot, M.H.M., Hooghiemstra, H., Bogotá-Angel, R.G., Berrio, J.C., 2012. Changing depositional environments in the Colombian Fúquene Basin at submillennial time-scales during 284-27 ka from unmixed grain-size distributions and aquatic pollen. Neth. J. Geosci. 91, 199-214. https://doi.org/10.1017/S0016774600001591.

Zomer, R.J., Xu, J., Trabucco, A., 2022. Version 3 of the global aridity index and potential evapotranspiration database. Sci. Data 9, 409. https://doi.org/10.1038/s41597-022-01493-1.