Temperature and Humidity as Drivers of Elevational Shifts in the Montane Forests-Paramo Boundary Over the Last 15,000 Years in the Northwest Andes

Temperature and Humidity as Drivers of Elevational Shifts in the Montane Forests-Paramo Boundary Over the Last 15,000 Years in the Northwest Andes | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Article Temperature and Humidity as Drivers of Elevational Shifts in the Montane Forests-Paramo Boundary Over the Last 15,000 Years in the Northwest Andes Andrés Silva-Duque, Sebastian González-Caro, César Velásquez, and 3 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-7368410/v1 This work is licensed under a CC BY 4.0 License Status: Under Review Version 1 posted 17 You are reading this latest preprint version Abstract In this study, conducted in the Belmira paramo in the northwestern Andes of Colombia, we used pollen records to track vegetation changes over the past 15,000 years before present (BP). Employing ordination (NMDS), species indicator, and wavelet analysis, we found a pattern of ecosystem change that shifted from a paramo to a tree-dominated community, which became a forest/paramo transitional ecosystem again. At the whole pollen community level, our results revealed spatially synchronous processes (i.e., temporal co-evolution) at regional scales between changes in pollen composition and mean annual temperature (MAT) and Fe concentrations (proxy of precipitation). Independent univariate wavelet analyses showed similar periodicities between 500 and 1000 years for MAT and the orthogonal compositional NMDS axes, while Fe exhibited significant variance at longer periodicities, ranging from 2000 to 4000 years. The strong and significant covariation assessed using wavelet coherence between NMDS II and Fe further confirmed the role of precipitation in shaping tree species turnover on cycles of approximately 1,000-500 years. However, our study points to the existence of a lag in vegetation tracking to climate change of centuries, which could constitute a major challenge for the implementation of effective nature-based strategies of conservation to ameliorate the ongoing global change. Earth and environmental sciences/Climate sciences Biological sciences/Ecology Earth and environmental sciences/Ecology Earth and environmental sciences/Environmental sciences Global climatic events Paramo Upper Forest Line Wavelet coherence Species turnover Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Introduction Palynological records have been used to analyze past shifts in elevational boundaries between highland forest and paramo ecosystems in response to temperature variation over geological time 1 – 3 . In the tropical Andes, global climatic events such as the Bølling–Allerød (BA) (14,600–12,900 y BP) and the Younger Dryas (YD) (12,900–11,700 y BP) 4 played a key role in shaping changes in the elevational boundary between forest and paramo ecosystems 5 , 6 . However, temperature changes may have been mirrored by fluctuations in other climatological variables, such as humidity 7 , precipitation, and the partial pressure of carbon dioxide (pCO 2 8 ), hampering to identify the extent to which each factor independently contributes to control shifts in the upper and lower boundary of forest and paramo ecosystems, respectively 9 . Assessing how temperature could interact with humidity or precipitation as drivers of the upper continuous forest boundary in tropical mountain ecosystems will improve our understanding of the potential response of highland tropical natural ecosystems to the ongoing climate change. In tropical paleoecology, a common approach to inferring changes in forest cover and plant community composition in response to climatic fluctuations has been the qualitative (visual) interpretation of the paired variation between pollen abundance of specific taxa and environmental proxies for temperature or humidity changes 7 . The use of environmental proxies is due to the lack of past temperature and precipitation records, which forces us to lay hands-on surrogate variables to reconstruct past climate change. For this reason, in tropical forests, the upper forest line (UFL) inferred from pollen abundance has often been used as a proxy for temperature 3 , 10 . Similarly, variation in iron ( Fe ) and titanium ( Ti ) concentrations, organic matter content, and the abundance of some specific pollen taxa (e.g., Isoetes , Plantago , etc.), have commonly been employed as proxies of precipitation and humidity 11 , 12 . In the tropical Andes, the upper montane forests and paramos provide essential ecosystem services, including water provision 13 , 14 and carbon storage 15 , 16 . Upper montane forests typically occur from approximately 2300–2500 masl up to around 3200–3500 masl, where the upper forest line marks the transition to the shrub and herb dominated vegetation characteristic of paramo ecosystems 7 . This elevational forest band has served as a corridor for the migration of temperate-affiliated taxa, such as Alnus accuminata , Drymis sp. , Podocarpus sp. , Quercus humboldtii , and Weinmania sp. among others 17 , 18 . In contrast, the paramo ecosystem is a highly distinctive highland environment of the Northern Andes, characterized by low to freezing temperatures, intense UV radiation, and iconic plant species like Espeletia 19 . Paramos, which have been estimated to have emerged around 6–10 million years ago 20 , 21 , are considered as the coldest biodiversity hotspot on earth 19 . However, due to land use and climate change, the highland Andean ecosystems are experiencing rapid shifts in species distributions 22 , 23 , reductions in their geographic extent 24 , 25 , and alterations in hydrological patterns 26 . Gaining insights into past environmental drivers of changes in the distribution of upper montane forests and paramos may improve our ability to predict their response to future climate scenarios, such as the projected increase of temperature and precipitation seasonality. Over the past 15,000 years before present (BP), tropical Andean ecosystems have undergone significant structural and plant compositional changes in response to major climatic events like the already mentioned BA and YD. For instance, the YD is known to have caused a decrease in precipitation and humidity, which in turn, influenced the distribution of upper montane forests and tropical alpine ecosystems such as paramos 27 , 28 . Likewise, during the middle and late Holocene (between 8,200 and 1,200 cal y BP), solar radiation and oceanic activity have been identified as potential drivers of external climatic forcing 29 . While solar variation proxy records (e.g., 14C, 10 Be) suggest cycles of 1–2 thousand-years, the likely influence of stochastic processes raises questions about the existence of unknown mechanisms operating at these temporal scales 30 . Building on previous findings 5 , we predict a strong association between plant community structure/composition and temperature/precipitation variation. However, the temporal changes in plant communities may not align exactly with environmental shifts due to the time required for species to migrate, adapt or go locally extinct 31 , 32 . This highlights the importance of assessing the potential temporal lags between these changes to better quantify the response time of plant communities to climate change. The traditional visual and descriptive interpretation of correlations between pollen abundance/composition and environmental proxies as a signal of climate change, can be significantly improved by the application of spectral analysis 31 , 33 . Wavelet analysis is a technique employed to decompose a signal, such as a non-stationary time series, into subsets of dilated and translated functions known as mother wavelets. This process enables the identification of oscillating components at different frequencies. By partitioning the variance of a time series into multiple components, wavelet analysis can detect energy peaks within the series 34 , 35 . In most cases, this methodology has been mainly used to identify the univariate temporal shifts in pollen composition 3 , 36 . Therefore, it is striking that the use of the bivariate spectral correlation using wavelet coherence, which enables to quantify the paired correlation between time series 37 , has been largely ignored. The present study was conducted in the Belmira region of the Northwestern Andes of Colombia (Fig. S1 ), a landscape characterized by a mosaic of upper montane forests and paramo ecosystems dominated by Q. humboldtii and Espeletia sp , respectively. The Belmira paramo is one of the lowest and wettest paramos in Colombia (lower boundary around ca. 3000 masl (Fig. S1 ), presenting a unique setting to examine the establishment of the upper forest limit and the dominance of grassland vegetation. In this paramo, which was covered by a large lake during the Pleistocene-Holocene transition, we use pollen records, common statistical tools employed in ecological studies, and spectral decomposition analysis (i.e., wavelet coherence), to improve the quantitative interpretation of coupled changes between pollen abundance/composition change and climate variability across different time scales. The main research questions are: a) How have the upper forest line (UFL) and lower paramo limit (LPL) changed over the last 15,000 y BP, and how have these shifts influenced plant community composition? b) Are changes in forest pollen abundance and composition significantly associated with variations in temperature and/or humidity? c) To what extent can we identify periodic changes in the paleoclimatic record associated with major climatic fluctuations, such as glacial and interglacial periods? d) Are there temporal lags in the response of plant communities to environmental changes? Results Sample dates Based on the stratigraphic inspection and the Accelerator Mass Spectrometry (AMS) 14 C samples, we identified changes in sediment composition over time. The BACON model 38 employed to analyze the sedimentary dynamics of the core, showed that the whole record spanned the last 14,735 cal y BP (Fig. S2, Table S1 ). From 187 cm to 161 cm depth, the sediment was composed predominantly of organic mud, followed by a 10 cm layer of peat. Between 150 and 50 cm, the sediment consisted of organic slit, which transitioned into peat in the uppermost 50 cm of the core. Between 14,700 and 13,500 cal y BP, the sedimentation rate reached 0.0135 cm·yr - 1 at 14,300 cal y BP. Thereafter, it remained relatively constant at approximately 0.0015 cm·yr - 1 between 13,500 and 10,500 cal y BP. During 10,500 and 9,000 cal y BP, a peak rate of 0.028 cm·yr - 1 was recorded at 9,400 cal y BP. Between 9,000 and 4,000 cal y BP, the accumulation rate fluctuated between 0.008 and 0.005 cm yr - 1 . The highest accumulation rate occurred from 4,000 cal y BP to the present, reaching approximately 0.02 cm yr - 1 (Fig. S3). Pollen compositional patterns A total of 91 palynological morphotypes were identified across the 187 samples, encompassing pollen types from 59 genera and 50 families. However, only 46 taxa were present in three or more samples and were selected to run the analyses. Of these, only Quercus humboldtii and Poaceae were consistently present in all samples throughout the entire record. The mean number of taxa per sample was of 24 ± 4. The NDMS analysis based on the 46 morphotypes yielded a convergent solution (stress = 0.05) after 20 iterations. The NMDS I ordination axis showed a temporal gradient ranging from a modern forest-paramo ecosystem to a forest dominated phase in the middle of the record, followed by a forest/paramo transition. The NMDS II showed a trend of change from forest to paramo from the bottom to the top of the ordination diagram (y-axis), respectively (Fig. 1 ), but emphasizing temporal changes of arboreal composition through time. Based on the time-constrained cluster classification applied on the NMDS ordination, we identified five compositional indicator groups of pollen taxa (Table S2). The two most distinct compositional groups were the first and the fifth, representing compositionally different paramo ecosystems located at opposites extremes of the first NMDS axis. The older paramo ecosystem spans a period from 14,700 to 11,000 cal y BP and is characterized by indicator taxa typical of herbaceous communities, such as Valeriana, Acaena, Puya and Vallea . In contrast, the modern paramo (extending from 2,300 cal y BP to the present) shows an increase of some paramo-related taxa, such as Iridiaceae, Rosaceae, Ericaceae and Hypericum , along with a rise in sub-Andean taxa, such as Morella and Clusia. The remaining three groups, which span the interval from 10,400 to 2,300 cal y BP, represent transitional shifts in forest species composition. Although distinct, these groups shared several of the most abundant forest taxa. During this period, Andean Forest taxa reached a peak in tree pollen abundance. The transition from the preceding zone (the older paramo) is marked by the increasing presence of arboreal taxa such as Gordonia , Juglans, Podocarpus and Drymis (Fig. 2 ). Univariate periodic decomposition The wavelet decomposition of the Fe time series, employed as a proxy of precipitation/humidity, showed a maximum average power spectral density in the 3000-year periodicity band between 12,000 and 8,000 cal y BP. During this period, additional regions of high variance were detected in the 2048-year periodicity band. Within the high-resolution sedimentary core, Fe fluctuations during the last glacial interval exhibited oscillations ranging from 512 to 1000 years, all of them encompassed within a broader 2048-year cycle. After 8,000 cal y BP, both variance and periodicity patterns showed a reduction in precipitation compared to earlier periods. Specifically, the power spectrum within the 2048-year periodicity band declined, and the 512-year cycles nearly disappeared, representing the pattern of shift in precipitation regimes from the last glacial maximum to the Holocene (Fig. 3 a). Temperature changes inferred from the pollen record exhibited significant variability throughout the Holocene. Between 12,000 and 8,000 cal y BP, temperatures ranged from 9.4 to 15.4°C, with the wavelet transform showing a statistically significant region dominated by a 2048-year periodicity. From 8,000 to 6,000 cal y BP, the time series exhibited a shift in frequency, with oscillations occurring predominantly at intervals between 1,000 and 500 years, and temperatures fluctuating between 12.7 and 14.2°C. Between approximately 6,000 to 2,500 cal y BP, temperatures peaked at 15.7°C. However, during this interval, the power spectra corresponding to the 500 and 1000-year bands disappeared from the spectrogram, indicating a loss of regular periodicity in temperature variation (Fig. 3 b). The wavelet power spectra applied to the first two NMDS axes revealed a dominant 2,048-year periodicity throughout the entire time series for both axes. However, distinct patterns emerged at shorter timescales. The first axis (NMDS I) exhibited significant spectral power in the 1,000 -512-year band during two intervals: from 10,000 to 8,000 cal y BP, and from 4,000 cal y BP to the present (Fig. 3 c). In contrast, NMDS II showed peak variance at the 2,048 years between 10,000 and 4,000 cal y BP, along with a concurrent 512-year cycle appearing between 10,000–8,000 cal y BP (Fig. 5 d). Notably, NMDS II axis displayed statistically significant high-power spectra between 10,000–8,000 cal y BP, which aligned with periodicities in the abundance of taxa such as Euphorbiaceae, Juglans, Podocarpus , and Morella . This pattern suggests a potential turnover in arboreal community composition during this period (Fig. 3 d). Environmental and pollen compositional covariation The wavelet coherence analysis between Fe concentrations and NMDS I showed a strong statistical association at shorter periodicities (~ 512 years), particularly between 6,500 and 8,000 cal y BP (Fig. 4 a). Additionally, a statistically significant region of high wavelet coherence between Fe and NMDS II time series was observed between 6,000 and 8,500 cal y BP, spanning the 2000 to 512-year periodicity bands (Fig. 4 b). In both instances, a temporal lag was evident between Fe and the ordination axis, implying that fluctuation in Fe levels is a likely driver of plant community dynamics. However, between 6,000 to 8,000 cal y BP, NMDS II and Fe were out of phase with a lag of -π/2 or ¼ of a quasi-period within the 1000- and 512-year bands with changes in Fe associated with compositional variation following a lag of approximately one quarter of quasi-period (Fig. 5 b). Regarding the relationship between temperature (MAT) and NMDS I, we identified a statistically significant coherence region between 14,000 and 10,000 cal y BP, with dominant cycles at approximately 512 years and 1000 years (Fig. 4 c). Similarly, a broader band of high coherence values was found between NMDS II and temperature over the same interval, with a predominant periodicity of 500 years (Fig. 4 d). During most of the Holocene (from 8,000 to 2,000 cal y BP), we also detected a high coherence region with a predominant cycle of 1000 years approximately. The phase difference shows that, for most of the record, the two series were out of phase by π/2. This indicates that changes in temperature were generally followed by shifts in plant community composition (Fig. 5 c, d). The influence of either MAT or Fe as drivers of indicator taxa, assessed through wavelet coherence analysis on a 500-year filtered series, showed different taxa-specific responses to each environmental variable (Fig. 6 ). For example, taxa such as Valeriana and Cariophyllaceae, representative of a paramo flora community that shifted around 11,000 cal y BP, showed a positive relationship with increased humidity (as indicated by Fe concentrations). Likewise, the peaks of abundance of taxa such as Juglans , Podocarpus and Morella mirrored the high frequency patterns in the Fe series, explaining the strong coherence found on wavelet decomposition of both series between 8,500 and 6,000 cal y BP (Figure S4). Regarding temperature, indicator taxa such as Gordonia , Aragoa and Juglans tracked the main fluctuation patterns of the MAT series, consistent with the observed wavelet coherence results (Fig. S5). Discussion Our study aligns with previous studies in which historical climate variability during the Holocene operated as a key driver of elevational shifts in the boundary between upper montane forests (UMF) and paramos in the northern Andes 2 . Using a quantitative approach, such as wavelet coherence, we were able to statistically evaluate the covariation between time series associated with environmental drivers and both plant communities and individual indicator species over the past 15,000 cal y BP in the northern Andes. We found that both temperature and humidity played a key role in shaping the historical transitions between forests to paramo ecosystems, as well as on driving species turnover through time. However, in contrast to former studies carried out in the Eastern cordillera of the Colombian Andes, our findings suggest that, in the Central cordillera of the Colombian Andes, climatic periods as the Bølling–Allerød and the Younger Dryas were characterized by wet climate 6 , rather than dry as reported in the eastern flank 3 . Our results showed a pattern of ecosystem change that shifted from a paramo dominated by Poaceae, Asteraceae and Puya, to a tree-dominated community characterized by Morella sp, Podocarpus sp , and high abundance of large trees such as Q. humboldtii, Juglans sp, Weinmannia sp and Alchornea sp , many of which decreased in abundance or disappeared from the palynological record to let the system to become a forest/paramo transitional ecosystem again (see Fig. 2 ). Similar transitions and changes in the position of the upper forest line (UFL) have been reported in previous studies in northwest Andes 3 , 5 , 6 , 36 , and have been primarily attributed to temperature changes associated with major climatic oscillations, such as the Last Glacial Maximum and the Younger Dryas 27 , 28 . Our findings, however, emphasize a differential response of plant taxa to either temperature or precipitation, improving our understanding of the likely response of this ecosystem to droughts and/or temperature increases. As evidenced along the second NMDS axis, distinct compositional changes also occurred within the forest dominated system during periods with relatively stable temperature. For example, here we show how Alnus positively responded to changes in precipitation, while abundance of Quercus humboldtii rose along with the increase of temperature (Fig. 5 ). Yet the potential influence of unmeasured biological mechanisms on shaping tree turnover during the forest-dominated phase, such as negative density dependence 39 or ecological drift 40 , 41 remains untested. To more effectively unravel the likely effect of biotic interactions and ecological drift on long-term plant species turnover in relatively stable biomes, we would require a higher taxonomic resolution than that employed here and in most palaeoecological studies (i.e., mostly genus or family), which seems a difficult technical challenge to overcome. Our time-series decomposition of the pollen community in the Belmira paramo revealed huge variance in both environmental and compositional variation signals during the late Holocene. According to the wavelet spectrum analysis (see Fig. 3 ), a significant and pronounced shift in plant composition occurred around 10,000 cal y BP, as reflected in both the transitional-paramo axis (NMDS-I) and the forest-paramo axis (NMDS-II). While previous studies have linked changes in the upper forest line of the northern Andes to long-term astronomical orbital cycles of 100,000-year (eccentricity) and 41,000-year (obliquity) 36 , our findings support the existence of internal cyclical variations at shorter temporal scales within interstadial periods of the major climatic oscillations 3 , 31 . Independent univariate wavelet analyses showed similar periodicities between 500 and 1000 years for MAT and the orthogonal compositional NMDS axes, while Fe (used as a proxy for humidity) exhibited significant variance at longer periodicities, ranging from 2000 to 4000 years. Using pollen data to reconstruct past temperature assumes that species distributions remained in equilibrium with environmental conditions over time. This approach overlooks species dispersal limitations 42 and assumes a stable adiabatic rate throughout the entire study period, assumptions that may not hold across millennial timescales 2 . We definitively need more studies with finer temporal pollen resolution as well as the use of complementary proxies highly sensitive to temperature, such as chironomid assemblages 43 , to enhance our understanding of the extent to which subtle changes in climate dynamics can determine species turnover at different temporal and spatial scales. One of the main goals of our study was to use wavelet coherence (WCO) to provide a quantitative objective interpretation of covariation based on time series decomposition 34 , rather than a visual comparison between abundance of pollen and climatic proxies, as it has been done in most traditional palynological studies. At the whole pollen community level, our results revealed spatially synchronous processes (i.e., temporal co-evolution) at regional scales between changes in pollen composition (i.e., NMDS axes I and II) and environmental variables such as MAT and Fe concentrations (proxy of humidity). In the case of the NMDS I, which represented the major changes from one biome to another (paramo-forest-paramo), we observed significant synchrony with MAT at periodicities of 500–1000-year period mainly around 10000–12,000 y cal BP. This corresponds to the periods during which landscape transitioned from a paramo to a forested dominated ecosystem, implying an upward shift of the upper forest line in response to rising temperatures. That said, forest cover, and thus the UFL, expanded upslope in synchrony with warming. In contrast, humidity showed a stronger synchrony with NMDS II, especially at shorter periodicities (i.e., ≤ 500 years) between 6,000–8,000 y cal BP, a time in which the region was likely fully covered by forest. However, our study points to the existence of a lag in vegetation tracking to climate change of even centuries 31 , which constitutes a major challenge for the implementation of effective nature-based strategies of conservation to ameliorate the ongoing global change 42 . The temporal correlation between NMDS I and Fe revealed a significant synchrony at approximately 500-year periodicities, particularly between 6,000–8,000 y cal BP. This pattern closely matched the covariation observed between NMDS I and MAT, suggesting that Fe and MAT interacted synchronously to influence tree turnover within the forest dominated biome. The strong and significant covariation between NMDS II and Fe further confirmed the role of humidity in shaping tree species turnover but expanded its influence on longer cycles of approximately 1,000–2,000 years. This highlights humidity as a key driver of long-term compositional change in Andean Forest communities. Although pronounced changes in humidity might typically be expected to trigger biome-level transitions, as observed in lowland savannas (e.g., from a forest to a grassland), our findings suggest that precipitation variability in the study region did not play a role in determining ecosystemic transitions at the biome level as strong as temperature variability. This supports the view that significant changes in temperature play a dominant role in driving systematic species turnover and upslope migration in tropical mountain ecosystems 22 , 23 . One of the most striking contributions of this study lies in its methodological approach, which analyzes pollen community-level using common tool employed in modern community ecology. Specifically, we highlight: 1) the selection of ecological groups and indicator species using statistical methods rather than visual or subjective criteria; and 2) the application of a quantitative approach, wavelet coherence analysis, to assess statistical associations between environmental drivers and the pollen abundance of indicator taxa, instead of simple visual interpretation. This approach allowed us to conclude that the abundance of some taxa in these forests, such as Drymis sp and Quercus humboldtii , was primarily associated with temperature, while changes in the abundance of other taxa, such as Cyatheaceae, Podocarpus sp. , and Puya sp. , were strongly associated with precipitation (Figure S4). Some taxa, such as Melpomene , showed no clear association with either temperature or precipitation. Overall, applying wavelet coherence and power spectral decomposition to investigate the effects of past climate change on species abundance 31 offers a powerful tool for inferring how current plant populations may respond to the ongoing climate change. Although we do not present a formal statistical test to assess the influence of major climatic oscillations on the observed changes in pollen abundance and composition due to the lack of data to do so, our results allow us to identify some of the likely effects that these climatic events may have had on shaping vegetation dynamics. The increased abundance of taxa such as Poaceae, Asteraceae, Acaena , and Puya , along with a lowering of the UFL during the late glacial interval, mirrors vegetation trends previously reported in the Central and Eastern Andes during cold periods 28 . Between ~ 15,000–13,000 cal y BP, the Fe time series exhibited a region of significant high-power spectra that may correspond to the Guantiva interstadial (~ 14,200–11,200 cal y BP) and the Abra stadial events (~ 11,200–9,700 cal y BP), which are considered the South American analogs of the Boiling-Allerod and Younger-Dryas, respectively 27 . In our study area, a subtle reduction in the UFL seems to be associated with the Abra stadial period. This interpretation is supported by the downward shift of the UFL of around ~ 300 masl below the level observed at the end of the Guantiva interstadial. This shift likely resulted from a decrease in rainfall, accompanied by a modest increase in the abundance of cold-adapted taxa such as Poaceae and Asteraceae. The findings of this research corroborate the notion that YD event was characterized by wet conditions in the Colombian Andes, consistent with evidenced from other palynological data 5 , 6 , 12 . During the BA period, the patterns of precipitation were influenced by a northward shift of the Intertropical Convergence Zone (ITCZ); this shift may have increased moisture influx from the Amazon basin contributing to the prevailing wet conditions during the YD. Around 10,000 cal yr BP, a transition to drier conditions is indicated by precipitation proxies. This climatic shift has been attributed to the weakening of the South American Summer Monsoon (SASM), which reduced moisture transport from the Amazon basin to the Andean region 12 , 44 . Although our dataset is not extensive enough to directly assess the effect of recent anthropogenic warming, the long-term trend observed over the past 15,000 years suggests that continued temperature increase must lead to a renewed shift toward a forest dominated ecosystem. However, this expectation of an upslope migration of the upper forest line in response to rising temperatures contrasts with the transitional trend observed over the past 2,500 years, which indicates a reversal from forest to paramo. This suggests that ongoing human-induced climate change may override the natural ecological trajectory of recent millennia, potentially reversing the recent expansion of páramo ecosystems. Such a shift could have direct consequences for the provision of key ecosystem services, including water supply and regulation, which are critically linked to paramo vegetation. Methods Study area The sediment core was collected from El Morro mire (06°40´11.9˝N, 75°40´4.4˝W), part of Belmira (Antioquia, Colombia) paramo in the northern Andes. The local bedrock primarily consists of quartz diorite formed during the late Cretaceous period, alongside a metamorphic complex that dates from the Precambrian era to the middle Cretaceous period. The regional low fertility and acidic soils are characterized by high organic matter (OM) and aluminum content. The paramo spans an altitudinal range between 3000 and 3400 masl, encompassing steep slopes and a flat zone characterized by diverse glacial landforms. The area experiences a tropical alpine climate, with an average annual temperature of 12°C and diurnal temperature fluctuations of approximately 10°C. The annual precipitation in this region reaches around 2710 mm, and is closely influenced by mesoscale processes, particularly the movement of the Inter Tropical Convergence Zone (ITCZ). Consequently, two distinct rainy periods occur: one from April to May, when the ITCZ moves northward, and another from September to October, when the ITCZ shifts southwards 6 (Fig. 1 ). Data collection At El Morro mire, we recovered a 300 cm long sediment core using a hand-operated Russian sampler, which extracts 50cm long, half-cylindrical sections with a 5cm diameter. A first stratigraphical description of the sediments was photographed in the field to document relevant core data, such as color and texture. The cores were wrapped in plastic foils and protected by PVC guttering for transport. Samples were taken to the paleoecology Laboratory of the National University in Medellín, where they were stored in a dark room at 4°C. A total of 300 samples were extracted from the sediment core, maintaining a one-centimeter distance between the closest samples. These samples were treated with standard procedures for pollen preparation, including sodium pyrophosphate and acetolysis. Then, the pollen samples were prepared with glycerin gelatin and fixed in microscope slides with wax. Using pollen morphological studies and the pollen reference collection of the Paleoecology Laboratory of the National University of Colombia (Medellín) 6 , pollen grains of each sample were counted and identified to the highest possible taxonomic resolution (i.e., species, genus, or family). Lithological and stratigraphic characteristics were assessed through visual inspection of the core, recording details such as grain size, color, and strata contacts. A total of 17 samples were dated at Beta Analytic Laboratory (Miami). However, only eight samples (approximately located between 60 and 190 cm) were suitable for use in the age model. The remaining samples exhibited high uncertainties in their estimated reversal ages and were therefore excluded (Fig. S6). We used the RBacon R package 38 to calibrate the AMS radiocarbon ( 14 C) ages, estimate sediment accumulation rates, and to construct the age-depth model (Table S1 ). To estimate the age of the pollen samples, we combined radiocarbon dating and depth measurements. First, the radiocarbon dates provided by the laboratory were calibrated using the IntCal20 calibration curve, which is appropriate for samples from the North Hemisphere. Next, we established a relationship between the radiocarbon-based dates and their corresponding depths using a Bayesian autoregressive gamma process model. This model accounted for the autocorrelation among samples resulting from sedimentary processes. Furthermore, our model accounted for uncertainties associated with the radiocarbon dates, including age uncertainty and 14 C calibration curve uncertainty. It also incorporated sedimentation rates to enhance the accuracy of age estimation 45 , 46 . This analysis was performed using the Rbacon package in R 38 . Elemental analysis was performed on the core using the EDAX Eagle III XPL mProbe at the Earth Science Department of the University of Geneva. Material for X-ray microfluorescence (mXRF) analyses were sampled at the National University of Colombia in Medellín in 50 cm-long U-channels from the cores recovered with the Russian corer and transported to the University of Geneva. The U-channels were cut into 10 cm long sections to be analyzed in the Eagle III spectrometer. The parameters employed for these measurements were set at 40 kV, 450 mA, 50/55 Dtm at atmospheric pressure. Measured values are expressed in counts per second (cps). Titanium ( Ti ) and Iron ( Fe ) concentrations, which represent watershed sedimentary and deposit trajectories due to their relationship with rainfall and runoff processes (i.e., precipitation proxies) were used as chemical proxies to assess the effect of climatic variation on pollen composition during the late Quaternary 47 . Upper forest limit definition Using the pollen data, we built a community matrix in which rows represent each dated pollen sample (i.e., time), and columns correspond to the identified taxa (e.g., species, genus, or family). We classified plant taxa by growth form as tree, shrub, or herb, and use the percentage of arboreal pollen (%AP) to define the elevational position of the Upper Forest Limit (UFL) 2 (Flantua et al., 2019). We assume that each 5% change in %AP corresponds to a 100 m shift in the UFL. The equation employed to estimate the UFL was as follows 2 , 5 : $$\:UFL=\:ETL+\left(\frac{\%AP-C}{R}\right)\times\:\:100$$ Where ETL is the current elevation of the tree line (3076m asl) for Belmira Paramo based on a Paramo delimitation project; C is the current percentage of arboreal pollen in the study area; and R represents the rate of change of the UFL 2 , 5 . The variation of mean annual temperature (MAT) is related to the UFL as follows 2 , 5 : $$\:MAT=\:CMT+\left(UFL-ETL\right)\times\:\:AR$$ Where CMT represents the current elevation range of the Belmira Paramo and AR denotes the adiabatic rate. In our study, we used an upward displacement in elevation. This relationship indicates that the mean annual temperature tends to decrease with increasing elevation, with the Upper Forest Line serving as a critical point where significant temperature changes occur. Patterns of pollen compositional change We used Bray-Curtis distance to calculate pollen compositional dissimilarity between samples. To estimate this metric, we included only those taxa present in three or more samples and applied a Non-Metric Dimensional Scaling (NMDS) ordination to reduce the pollen community matrix into orthogonal axes. Taxa found in only one or two samples were excluded to: a) reduce statistical noise due to under sampling; b) identify groups of compositional change based on well represented taxa 48 ; and c) detect temporal periodicity to be incorporated into the wavelet analysis (see below). We used the first two NMDS axes (NMDS I and NMDS II) to assess plant community over the past 15,000 years BP. Next, we applied a time-constrained CONISS (Constrained Incremental Sum of Squares) clustering algorithm 49 to identify homogeneous floristic groups. Once these groups were defined, we performed a species indicator analysis 50 , 51 to identify the most representative taxa within each plant community. This analysis was conducted using the R packages vegan 52 and indicspecies 51 . Univariate changes of pollen composition and environmental drivers We applied the continuous wavelet transform 34 , 35 to detect periodic changes in plant composition, wet/dry periods, cycles, and seasonality in time series spanning the last 15,000 y. Specifically, we computed the convolution of each time series using the Morlet mother wavelet and estimated the mean power spectrum by averaging the wavelet coefficients within each frequency component. The Morlet mother wavelet is defined as follows: $$\:{\psi\:}_{0}\left(\eta\:\right)=\:{\pi\:}^{1/4}{e}^{-i{\omega\:}_{0}\eta\:}{e}^{-{\eta\:}^{2}/2}$$ Where \(\:i=\:\sqrt{-1}\) , \(\:{\omega\:}_{0}\:\) is ω 0 dimensionless frequency, and η is dimensionless time. The Morlet wavelet is particularly useful to analyze periodicities in time series. We assumed the standard value of ω 0 = 6 37 . We applied univariate wavelet transforms to the Fe, MAT, NMDS I and NMDS II time series to analyze temporal periodicities in precipitation, temperature, and pollen community composition. To assess the statistical significance of regions with high power spectra, we used a null model that preserves the short-term autocorrelation observed in the original time series. This approach allows us to evaluate the significance of long-term patterns in the data 33 . These analyses were performed using the set of MATLAB® functions provided by Cazelles et al (2008). Additionally, based on the wavelet transform of each taxon`s abundance time series, we applied hierarchical clustering to the resulting frequency decomposition matrix to identify taxa with similar temporal periodicities. Patterns of covariation using wavelet coherence analysis We used wavelet coherence analysis (WCO) to assess the paired variation between environmental proxies (Fe and MAT) and pollen compositional (NMDS I and NMDS II) time series. We also analyzed WCO between the environmental proxies and the pollen abundance of those taxa identified as indicators for each independent group defined by the cluster analysis, to evaluate the influence of temperature and/or precipitation on shaping the signal of each indicator taxon in the palynological record. This analysis allows us to identify regions of convergence in both time and frequency, where the power spectra of the two series exhibit a linear correlation. Furthermore, it quantifies the phase difference between the two series, providing information about potential time lags 53 . To statistically evaluate the significance of coherence, we compared the observed cross-wavelet spectra with those generated by a null model cross-wavelet spectrum. The null model was created by simulating independent red-noise processes with the same first order autoregressive coefficient as the original time series. This comparison allows us to determine whether the observed cross-wavelet spectra significantly deviated from expectations under the null hypothesis 37 . To enhance the interpretation of the correlation between environmental and compositional variation variables, we employed a low-pass Gaussian filter to both the proxy data (MAT and Fe) and the abundance series of indicator taxa. This filter was calibrated using a bandpass value derived from regions of high coherence, such as the 500-year filtered series. These analyses were also performed using the MATLAB® functions provided by Cazelles et al (2008). Declarations Additional information Competing financial interests: The authors declare no competing financial interests. Funding The authors received no funding for this work Author Contribution A.S-D., A.D., and C.V. conceived and designed the experiments. A.S-D, I.M.C-R, and R.E.V-M performed the experiments. A.S-D analyzed the data with advice from S.G-C. and A.D. A.S-D. and A.D. wrote the main text. All authors reviewed and commented on the final version of the manuscript. Acknowledgement We thank Professor Bernard Cazelles for providing the code to run univariate wavelet decomposition, wavelet coherence, and phase difference analysis. Data Availability The datasets used and/or analysed during the current study available from the corresponding author on reasonable request. References Körner, C. Alpine Treelines . (Springer, Basel, 2012). doi:10.1007/978-3-0348-0396-0. Flantua, S. G. A., O’Dea, A., Onstein, R. E., Giraldo, C. & Hooghiemstra, H. The flickering connectivity system of the north Andean páramos. Journal of Biogeography 46 , 1808–1825 (2019). Groot, M. H. M. et al. Ultra-high resolution pollen record from the northern Andes reveals rapid shifts in montane climates within the last two glacial cycles. Climate of the Past 7 , 299–316 (2011). Alley, R. B. Ice-core evidence of abrupt climate changes. Proc. Natl. Acad. Sci. U.S.A. 97 , 1331–1334 (2000). Muñoz, P. et al. Holocene climatic variations in the Western Cordillera of Colombia: A multiproxy high-resolution record unravels the dual influence of ENSO and ITCZ. Quaternary Science Reviews 155 , 159–178 (2017). Velásquez-R., C. A. & Hooghiemstra, H. Pollen-based 17-kyr forest dynamics and climate change from the Western Cordillera of Colombia; no-analogue associations and temporarily lost biomes. Review of Palaeobotany and Palynology 194 , 38–49 (2013). Hooghiemstra, H., Wijninga, V. M. & Cleef, A. M. THE PALEOBOTANICAL RECORD OF COLOMBIA: IMPLICATIONS FOR BIOGEOGRAPHY AND BIODIVERSITY1. mobt 93 , 297–325 (2006). Marchant, R. et al. Colombian vegetation at the Last Glacial Maximum: a comparison of model- and pollen-based biome reconstructions. Journal of Quaternary Science 19 , 721–732 (2004). Bader, M. Y. et al. A global framework for linking alpine-treeline ecotone patterns to underlying processes. Ecography 44 , 265–292 (2021). Bush, M. B. Distributional change and conservation on the Andean flank: a palaeoecological perspective. Global Ecology and Biogeography 11 , 463–473 (2002). Martínez, C. et al. Neogene precipitation, vegetation, and elevation history of the Central Andean Plateau. Sci. Adv. 6 , eaaz4724 (2020). Patiño, L. et al. Late Pleistocene–Holocene environmental and climatic history of a freshwater paramo ecosystem in the northern Andes. Journal of Quaternary Science 35 , 1046–1056 (2020). Buytaert, W. & De Bièvre, B. Water for cities: The impact of climate change and demographic growth in the tropical Andes. Water Resources Research 48 , (2012). Immerzeel, W. W. et al. Importance and vulnerability of the world’s water towers. Nature 577 , 364–369 (2020). Duque, A. et al. Mature Andean forests as globally important carbon sinks and future carbon refuges. Nat Commun 12 , 2138 (2021). Peña, M. A., Feeley, K. J. & Duque, A. Effects of endogenous and exogenous processes on aboveground biomass stocks and dynamics in Andean forests. Plant Ecology 219 , 1481–1492 (2018). Segovia, R. A. et al. Freezing and water availability structure the evolutionary diversity of trees across the Americas. Sci. Adv. 6 , eaaz5373 (2020). González-Caro, S. et al. The legacy of biogeographic history on the composition and structure of Andean forests. Ecology 101 , e03131 (2020). Madriñán, S., Cortés, A. J. & Richardson, J. E. Páramo is the world’s fastest evolving and coolest biodiversity hotspot. Front. Genet. 4 , (2013). Pouchon, C. et al. Phylogenomic Analysis of the Explosive Adaptive Radiation of the Espeletia Complex (Asteraceae) in the Tropical Andes. Systematic Biology 67 , 1041–1060 (2018). Pouchon, C. et al. Phylogenetic signatures of ecological divergence and leapfrog adaptive radiation in Espeletia. American Journal of Botany 108 , 113–128 (2021). Duque, A., Stevenson, P. R. & Feeley, K. J. Thermophilization of adult and juvenile tree communities in the northern tropical Andes. Proceedings of the National Academy of Sciences 112 , 10744–10749 (2015). Fadrique, B. et al. Widespread but heterogeneous responses of Andean forests to climate change. Nature 564 , 207–212 (2018). Peyre, G. What Does the Future Hold for Páramo Plants? A Modelling Approach. Front. Ecol. Evol. 10 , (2022). Peyre, G. et al. The fate of páramo plant assemblages in the sky islands of the northern Andes. Journal of Vegetation Science 31 , 967–980 (2020). Buytaert, W. et al. Human impact on the hydrology of the Andean páramos. Earth-Science Reviews 79 , 53–72 (2006). Van Der Hammen, T. & Hooghiemstra, H. The EL abra stadial, a younger dryas equivalent in Colombia. Quaternary Science Reviews 14 , 841–851 (1995). Veer, R. van ’t, Islebe, G. A. & Hooghiemstra, H. Climatic change during the Younger Dryas chron in northern South America: a test of the evidence. Quaternary Science Reviews 19 , 1821–1835 (2000). Debret, M. et al. The origin of the 1500-year climate cycles in Holocene North-Atlantic records. Climate of the Past 3 , 569–575 (2007). Turner, T. E. et al. Solar cycles or random processes? Evaluating solar variability in Holocene climate records. Sci Rep 6 , 23961 (2016). Fastovich, D. et al. Coupled, decoupled, and abrupt responses of vegetation to climate across timescales. Science 389 , 64–68 (2025). Knight, C. A. et al. Community Assembly and Climate Mismatch in Late Quaternary Eastern North American Pollen Assemblages. The American Naturalist 195 , 166–180 (2020). Cazelles, B. et al. Wavelet Analysis of Ecological Time Series. Oecologia 156 , 287–304 (2008). Grinsted, A., Moore, J. C. & Jevrejeva, S. Application of the cross wavelet transform and wavelet coherence to geophysical time series. Nonlinear Processes in Geophysics 11 , 561–566 (2004). Torrence, C. & Webster, P. J. Interdecadal Changes in the ENSO–Monsoon System. (1999). Torres, V., Hooghiemstra, H., Lourens, L. & Tzedakis, P. C. Astronomical tuning of long pollen records reveals the dynamic history of montane biomes and lake levels in the tropical high Andes during the Quaternary. Quaternary Science Reviews 63 , 59–72 (2013). Benincà, E., Ballantine, B., Ellner, S. P. & Huisman, J. Species fluctuations sustained by a cyclic succession at the edge of chaos. Proc. Natl. Acad. Sci. U.S.A. 112 , 6389–6394 (2015). Blaauw, M. & Christen, J. A. Flexible paleoclimate age-depth models using an autoregressive gamma process. Bayesian Analysis 6 , 457–474 (2011). Volkov, I., Banavar, J. R., He, F., Hubbell, S. P. & Maritan, A. Density dependence explains tree species abundance and diversity in tropical forests. Nature 438 , 658–661 (2005). Correa-Metrio, A., Meave, J. A., Lozano-García, S. & Bush, M. B. Environmental determinism and neutrality in vegetation at millennial time scales. Journal of Vegetation Science 25 , 627–635 (2014). Hubbell, S. P. Neutral Theory in Community Ecology and the Hypothesis of Functional Equivalence. Functional Ecology 19 , 166–172 (2005). Svenning, J.-C. & Sandel, B. Disequilibrium vegetation dynamics under future climate change. American Journal of Botany 100 , 1266–1286 (2013). Matthews-Bird, F., Brooks, S. J., Holden, P. B., Montoya, E. & Gosling, W. D. Inferring late-Holocene climate in the Ecuadorian Andes using a chironomid-based temperature inference model. Climate of the Past 12 , 1263–1280 (2016). Rodríguez-Zorro, P. A. et al. Shut down of the South American summer monsoon during the penultimate glacial. Sci Rep 10 , 6275 (2020). Blaauw, M. Methods and code for ‘classical’ age-modelling of radiocarbon sequences. Quaternary Geochronology 5 , 512–518 (2010). Andrés Christen, J. & Pérez E, S. A New Robust Statistical Model for Radiocarbon Data. Radiocarbon 51 , 1047–1059 (2009). Haug, G. H., Hughen, K. A., Sigman, D. M., Peterson, L. C. & Röhl, U. Southward Migration of the Intertropical Convergence Zone through the Holocene. Science 293 , 1304–1308 (2001). Legendre, P. & Legendre, L. Numerical Ecology . (Elsevier, 2012). Grimm, E. C. CONISS: a FORTRAN 77 program for stratigraphically constrained cluster analysis by the method of incremental sum of squares. Computers & Geosciences 13 , 13–35 (1987). Cáceres, M. D. & Legendre, P. Associations between species and groups of sites: indices and statistical inference. Ecology 90 , 3566–3574 (2009). De Cáceres, M., Legendre, P., Wiser, S. K. & Brotons, L. Using species combinations in indicator value analyses. Methods in Ecology and Evolution 3 , 973–982 (2012). Oksanen, J. et al. vegan: Community Ecology Package. (2025). Maraun, D. & Kurths, J. Cross wavelet analysis: significance testing and pitfalls. Nonlinear Processes in Geophysics 11 , 505–514 (2004). Additional Declarations No competing interests reported. Supplementary Files Silvaetal20250812SI.docx Cite Share Download PDF Status: Under Review Version 1 posted Editorial decision: Revision requested 17 Nov, 2025 Reviews received at journal 17 Nov, 2025 Reviews received at journal 14 Nov, 2025 Reviews received at journal 11 Nov, 2025 Reviews received at journal 10 Nov, 2025 Reviews received at journal 13 Oct, 2025 Reviewers agreed at journal 16 Sep, 2025 Reviewers agreed at journal 14 Sep, 2025 Reviewers agreed at journal 14 Sep, 2025 Reviewers agreed at journal 13 Sep, 2025 Reviewers agreed at journal 13 Sep, 2025 Reviewers agreed at journal 12 Sep, 2025 Reviewers invited by journal 11 Sep, 2025 Editor assigned by journal 11 Sep, 2025 Editor invited by journal 01 Sep, 2025 Submission checks completed at journal 21 Aug, 2025 First submitted to journal 21 Aug, 2025 You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-7368410","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Article","associatedPublications":[],"authors":[{"id":515856227,"identity":"1720829b-e0bf-4765-99f3-2396af2fd982","order_by":0,"name":"Andrés Silva-Duque","email":"","orcid":"","institution":"Universidad Nacional de Colombia Sede Medellín","correspondingAuthor":false,"prefix":"","firstName":"Andrés","middleName":"","lastName":"Silva-Duque","suffix":""},{"id":515856228,"identity":"7e645434-c632-41a6-a153-084e72e363dc","order_by":1,"name":"Sebastian González-Caro","email":"","orcid":"","institution":"Universidad Nacional de Colombia Sede Medellín","correspondingAuthor":false,"prefix":"","firstName":"Sebastian","middleName":"","lastName":"González-Caro","suffix":""},{"id":515856229,"identity":"22188a90-ec7e-4424-9981-8cab15876dae","order_by":2,"name":"César Velásquez","email":"","orcid":"","institution":"Universidad Nacional de Colombia Sede Medellín","correspondingAuthor":false,"prefix":"","firstName":"César","middleName":"","lastName":"Velásquez","suffix":""},{"id":515856230,"identity":"f918b8ef-53ab-465b-8072-4f7221dd1525","order_by":3,"name":"Ivonne Marcela Castañeda Riascos","email":"","orcid":"","institution":"Universidad Nacional de Colombia Sede Medellín","correspondingAuthor":false,"prefix":"","firstName":"Ivonne","middleName":"Marcela Castañeda","lastName":"Riascos","suffix":""},{"id":515856231,"identity":"66d6d0b2-e3e6-4efb-b68b-a7af9137b286","order_by":4,"name":"Rosa Elena Velásquez Montoya","email":"","orcid":"","institution":"Universidad Nacional de Colombia Sede Medellín","correspondingAuthor":false,"prefix":"","firstName":"Rosa","middleName":"Elena Velásquez","lastName":"Montoya","suffix":""},{"id":515856232,"identity":"9ad9f6ce-7ad2-49c9-945c-94e6160f4104","order_by":5,"name":"Álvaro Duque","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA2klEQVRIiWNgGAWjYFCCAwYgUg6IGQ8wMDATr8UYzCRSCwNYS2ID0Vr4Gw9vk3jbZpO+4XbzgwMMFdaJDeyHD+DVInHgWJnk3La03A13jhkcYDiTntjAk5aA35oDZ8ykedsO5264kWBwgLHtcGKDBI8BXh3yEC3/0w1upH84wPiPCC0GEC0HEgxu5ABtaSBCi+GBY8WWc84lG868kVNwIOFYunEbIb/I3Ti88cabMjt5vhvpGx98qLGW7ScUYsAgY2DggXFAxrPhVw8E/A1IWkbBKBgFo2AUYAMAjVxOoaaKOsgAAAAASUVORK5CYII=","orcid":"","institution":"Universidad Nacional de Colombia Sede Medellín","correspondingAuthor":true,"prefix":"","firstName":"Álvaro","middleName":"","lastName":"Duque","suffix":""}],"badges":[],"createdAt":"2025-08-13 23:23:08","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-7368410/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-7368410/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":91616851,"identity":"ce130819-b095-4150-86e2-d922a3f0c35e","added_by":"auto","created_at":"2025-09-18 10:41:02","extension":"jpeg","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":126849,"visible":true,"origin":"","legend":"\u003cp\u003eNon-metric dimensional scaling (NMDS) scatterplot of the two first ordination axes. The color palette represents the corresponding time sampled sites. The ellipses encompass sites belonging to each temporal cluster. The size of points represents tree pollen percentage.\u003c/p\u003e","description":"","filename":"floatimage1.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7368410/v1/7ccfa67213139404ca13a6fb.jpeg"},{"id":91616854,"identity":"bfca15a5-7397-4357-899e-a8cda0440d27","added_by":"auto","created_at":"2025-09-18 10:41:02","extension":"jpeg","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":390147,"visible":true,"origin":"","legend":"\u003cp\u003eTemporal time series inferred from sedimentary core. Left to right, representation of changes of sediments units along the core. Temporal variation of upper forest line/temperature (a) and Fe (b), community matrix ordination axis (NMDS I and NMDS II). Relative abundance of representative taxa indicator for each cluster: \u003cem\u003eValeriana\u003c/em\u003e, Solanaceae, \u003cem\u003eAcaena\u003c/em\u003e, Cariophyllaceae, \u003cem\u003ePuya, Aragoa, Iridiaceae, Gordonia, Juglans, Melpomene\u003c/em\u003e, Euphorbiaceae, \u003cem\u003eDrymis, Geranium, Bromellia\u003c/em\u003e, Myrtaceae, \u003cem\u003eRumex, Morella and Clusia.\u003c/em\u003e\u003c/p\u003e","description":"","filename":"floatimage2.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7368410/v1/5aa5e659f1039ae0b915b5d6.jpeg"},{"id":91616855,"identity":"96706d52-515c-4842-a170-130031c6d844","added_by":"auto","created_at":"2025-09-18 10:41:02","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":424290,"visible":true,"origin":"","legend":"\u003cp\u003eFe (a), MAT (b), NMDS I (c) and NMDS II (d) timeseries decomposition. Spectrogram represents the scale of mother wavelet period (y axis) relate with time. Heat map pixels indicate the power spectra of series. Black contorts lines identify the regions of statistical significance in the signal.\u003c/p\u003e","description":"","filename":"floatimage3.png","url":"https://assets-eu.researchsquare.com/files/rs-7368410/v1/cc665e4bd52c319760a7c856.png"},{"id":91616862,"identity":"0c2322bc-a65e-4f31-aae0-7e878ba55864","added_by":"auto","created_at":"2025-09-18 10:41:02","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":331334,"visible":true,"origin":"","legend":"\u003cp\u003eWavelet coherence between environmental variables and ordination axis a) Fe-NMDS I, b) Fe – NMDS II, c) MAT – NMDS I, d) MAT – NMDSII. Spectrogram represents the scale of mother wavelet (y axis) relate with time. Heat map pixels indicate the coherence between each pair of time series. Black contorts lines identify the regions of statistical significance in the signal (alpha = 5%).\u003c/p\u003e","description":"","filename":"floatimage4.png","url":"https://assets-eu.researchsquare.com/files/rs-7368410/v1/24f8ea0db785ede19caff78e.png"},{"id":91616857,"identity":"2ba62a44-c62f-4f14-912b-812fffbbf58f","added_by":"auto","created_at":"2025-09-18 10:41:02","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":374581,"visible":true,"origin":"","legend":"\u003cp\u003ePhases of ordination time series and environmental surrogated variables a) Fe – NMDS I, b). Fe-NMDS II, c) MAT – NMDS I, d) MAT – NMDS II. The solid lines represent the phases of each pair of series for 500-to-1000-year period bands. Dotted line indicates the phase difference between both series. The histogram of left illustrated the distribution of phase difference between series along the time span.\u003c/p\u003e","description":"","filename":"floatimage5.png","url":"https://assets-eu.researchsquare.com/files/rs-7368410/v1/af237222afe51b65a65660bf.png"},{"id":91616866,"identity":"96fddac3-64bb-4b4a-957f-a48488f8532c","added_by":"auto","created_at":"2025-09-18 10:41:02","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":291894,"visible":true,"origin":"","legend":"\u003cp\u003eNormalized time series of MAT, Fe and some statistically significative indicator taxa. Each line represents a low filtered series using 500-year periodicities. In each panel, black line refers to one specific taxon while red and blue lines indicate the Fe and MAT filtered series respectively. The shaded polygon indicates regions where statistically significant WCO values were identified between each pair of variables.\u003c/p\u003e","description":"","filename":"floatimage6.png","url":"https://assets-eu.researchsquare.com/files/rs-7368410/v1/ac3f04968411dfc0a559e214.png"},{"id":91618353,"identity":"7a3d65ca-aab6-4d29-abc6-020398b18707","added_by":"auto","created_at":"2025-09-18 11:05:04","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":2885113,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7368410/v1/bc9d966e-1b46-4ad7-b52f-93788cc0012d.pdf"},{"id":91616852,"identity":"b1b30a00-a41e-4e2d-b2b9-7e72a8f8a4cb","added_by":"auto","created_at":"2025-09-18 10:41:02","extension":"docx","order_by":0,"title":"","display":"","copyAsset":false,"role":"supplement","size":1847875,"visible":true,"origin":"","legend":"","description":"","filename":"Silvaetal20250812SI.docx","url":"https://assets-eu.researchsquare.com/files/rs-7368410/v1/bfc442c95148e8b3e2bb9364.docx"}],"financialInterests":"No competing interests reported.","formattedTitle":"Temperature and Humidity as Drivers of Elevational Shifts in the Montane Forests-Paramo Boundary Over the Last 15,000 Years in the Northwest Andes","fulltext":[{"header":"Introduction","content":"\u003cp\u003ePalynological records have been used to analyze past shifts in elevational boundaries between highland forest and paramo ecosystems in response to temperature variation over geological time \u003csup\u003e\u003cspan additionalcitationids=\"CR2\" citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u003c/sup\u003e. In the tropical Andes, global climatic events such as the B\u0026oslash;lling\u0026ndash;Aller\u0026oslash;d (BA) (14,600\u0026ndash;12,900 y BP) and the Younger Dryas (YD) (12,900\u0026ndash;11,700 y BP) \u003csup\u003e\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e\u003c/sup\u003e played a key role in shaping changes in the elevational boundary between forest and paramo ecosystems\u003csup\u003e\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e,\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e\u003c/sup\u003e. However, temperature changes may have been mirrored by fluctuations in other climatological variables, such as humidity\u003csup\u003e\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e\u003c/sup\u003e, precipitation, and the partial pressure of carbon dioxide (pCO\u003csub\u003e2\u003c/sub\u003e\u003csup\u003e8\u003c/sup\u003e), hampering to identify the extent to which each factor independently contributes to control shifts in the upper and lower boundary of forest and paramo ecosystems, respectively\u003csup\u003e\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e\u003c/sup\u003e. Assessing how temperature could interact with humidity or precipitation as drivers of the upper continuous forest boundary in tropical mountain ecosystems will improve our understanding of the potential response of highland tropical natural ecosystems to the ongoing climate change.\u003c/p\u003e\u003cp\u003eIn tropical paleoecology, a common approach to inferring changes in forest cover and plant community composition in response to climatic fluctuations has been the qualitative (visual) interpretation of the paired variation between pollen abundance of specific taxa and environmental proxies for temperature or humidity changes\u003csup\u003e\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e\u003c/sup\u003e. The use of environmental proxies is due to the lack of past temperature and precipitation records, which forces us to lay hands-on surrogate variables to reconstruct past climate change. For this reason, in tropical forests, the upper forest line (UFL) inferred from pollen abundance has often been used as a proxy for temperature\u003csup\u003e\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e,\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e\u003c/sup\u003e. Similarly, variation in iron (\u003cem\u003eFe\u003c/em\u003e) and titanium (\u003cem\u003eTi\u003c/em\u003e) concentrations, organic matter content, and the abundance of some specific pollen taxa (e.g., \u003cem\u003eIsoetes\u003c/em\u003e, \u003cem\u003ePlantago\u003c/em\u003e, etc.), have commonly been employed as proxies of precipitation and humidity\u003csup\u003e\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e,\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eIn the tropical Andes, the upper montane forests and paramos provide essential ecosystem services, including water provision\u003csup\u003e\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e,\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e\u003c/sup\u003e and carbon storage\u003csup\u003e\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e,\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e\u003c/sup\u003e. Upper montane forests typically occur from approximately 2300\u0026ndash;2500 masl up to around 3200\u0026ndash;3500 masl, where the upper forest line marks the transition to the shrub and herb dominated vegetation characteristic of paramo ecosystems \u003csup\u003e\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e\u003c/sup\u003e. This elevational forest band has served as a corridor for the migration of temperate-affiliated taxa, such as \u003cem\u003eAlnus accuminata\u003c/em\u003e, \u003cem\u003eDrymis sp.\u003c/em\u003e, \u003cem\u003ePodocarpus sp.\u003c/em\u003e, \u003cem\u003eQuercus humboldtii\u003c/em\u003e, and \u003cem\u003eWeinmania sp.\u003c/em\u003e among others \u003csup\u003e\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e,\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e\u003c/sup\u003e. In contrast, the paramo ecosystem is a highly distinctive highland environment of the Northern Andes, characterized by low to freezing temperatures, intense UV radiation, and iconic plant species like \u003cem\u003eEspeletia\u003c/em\u003e\u003csup\u003e\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e\u003c/sup\u003e. Paramos, which have been estimated to have emerged around 6\u0026ndash;10\u0026nbsp;million years ago\u003csup\u003e\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e,\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e\u003c/sup\u003e, are considered as the coldest biodiversity hotspot on earth\u003csup\u003e\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e\u003c/sup\u003e. However, due to land use and climate change, the highland Andean ecosystems are experiencing rapid shifts in species distributions\u003csup\u003e\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e,\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e\u003c/sup\u003e, reductions in their geographic extent\u003csup\u003e\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e,\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e\u003c/sup\u003e, and alterations in hydrological patterns\u003csup\u003e\u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e\u003c/sup\u003e. Gaining insights into past environmental drivers of changes in the distribution of upper montane forests and paramos may improve our ability to predict their response to future climate scenarios, such as the projected increase of temperature and precipitation seasonality.\u003c/p\u003e\u003cp\u003eOver the past 15,000 years before present (BP), tropical Andean ecosystems have undergone significant structural and plant compositional changes in response to major climatic events like the already mentioned BA and YD. For instance, the YD is known to have caused a decrease in precipitation and humidity, which in turn, influenced the distribution of upper montane forests and tropical alpine ecosystems such as paramos\u003csup\u003e\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e,\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e\u003c/sup\u003e. Likewise, during the middle and late Holocene (between 8,200 and 1,200 cal y BP), solar radiation and oceanic activity have been identified as potential drivers of external climatic forcing\u003csup\u003e\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e\u003c/sup\u003e. While solar variation proxy records (e.g., 14C, 10 Be) suggest cycles of 1\u0026ndash;2 thousand-years, the likely influence of stochastic processes raises questions about the existence of unknown mechanisms operating at these temporal scales \u003csup\u003e\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e\u003c/sup\u003e. Building on previous findings\u003csup\u003e\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u003c/sup\u003e, we predict a strong association between plant community structure/composition and temperature/precipitation variation. However, the temporal changes in plant communities may not align exactly with environmental shifts due to the time required for species to migrate, adapt or go locally extinct\u003csup\u003e\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e,\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e\u003c/sup\u003e. This highlights the importance of assessing the potential temporal lags between these changes to better quantify the response time of plant communities to climate change.\u003c/p\u003e\u003cp\u003eThe traditional visual and descriptive interpretation of correlations between pollen abundance/composition and environmental proxies as a signal of climate change, can be significantly improved by the application of spectral analysis\u003csup\u003e\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e,\u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e\u003c/sup\u003e. Wavelet analysis is a technique employed to decompose a signal, such as a non-stationary time series, into subsets of dilated and translated functions known as mother wavelets. This process enables the identification of oscillating components at different frequencies. By partitioning the variance of a time series into multiple components, wavelet analysis can detect energy peaks within the series\u003csup\u003e\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e,\u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e\u003c/sup\u003e. In most cases, this methodology has been mainly used to identify the univariate temporal shifts in pollen composition\u003csup\u003e\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e,\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e\u003c/sup\u003e. Therefore, it is striking that the use of the bivariate spectral correlation using wavelet coherence, which enables to quantify the paired correlation between time series\u003csup\u003e\u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e37\u003c/span\u003e\u003c/sup\u003e, has been largely ignored.\u003c/p\u003e\u003cp\u003eThe present study was conducted in the Belmira region of the Northwestern Andes of Colombia (Fig. \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e), a landscape characterized by a mosaic of upper montane forests and paramo ecosystems dominated by \u003cem\u003eQ. humboldtii\u003c/em\u003e and \u003cem\u003eEspeletia sp\u003c/em\u003e, respectively. The Belmira paramo is one of the lowest and wettest paramos in Colombia (lower boundary around ca. 3000 masl (Fig. \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e), presenting a unique setting to examine the establishment of the upper forest limit and the dominance of grassland vegetation. In this paramo, which was covered by a large lake during the Pleistocene-Holocene transition, we use pollen records, common statistical tools employed in ecological studies, and spectral decomposition analysis (i.e., wavelet coherence), to improve the quantitative interpretation of coupled changes between pollen abundance/composition change and climate variability across different time scales. The main research questions are: a) How have the upper forest line (UFL) and lower paramo limit (LPL) changed over the last 15,000 y BP, and how have these shifts influenced plant community composition? b) Are changes in forest pollen abundance and composition significantly associated with variations in temperature and/or humidity? c) To what extent can we identify periodic changes in the paleoclimatic record associated with major climatic fluctuations, such as glacial and interglacial periods? d) Are there temporal lags in the response of plant communities to environmental changes?\u003c/p\u003e"},{"header":"Results","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e\u003ch2\u003eSample dates\u003c/h2\u003e\u003cp\u003eBased on the stratigraphic inspection and the Accelerator Mass Spectrometry (AMS) \u003csup\u003e\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e\u003c/sup\u003eC samples, we identified changes in sediment composition over time. The BACON model\u003csup\u003e\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e\u003c/sup\u003e employed to analyze the sedimentary dynamics of the core, showed that the whole record spanned the last 14,735 cal y BP (Fig. S2, Table \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e). From 187 cm to 161 cm depth, the sediment was composed predominantly of organic mud, followed by a 10 cm layer of peat. Between 150 and 50 cm, the sediment consisted of organic slit, which transitioned into peat in the uppermost 50 cm of the core.\u003c/p\u003e\u003cp\u003eBetween 14,700 and 13,500 cal y BP, the sedimentation rate reached 0.0135 cm\u0026middot;yr\u003csup\u003e-\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e\u003c/sup\u003e at 14,300 cal y BP. Thereafter, it remained relatively constant at approximately 0.0015 cm\u0026middot;yr\u003csup\u003e-\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e\u003c/sup\u003e between 13,500 and 10,500 cal y BP. During 10,500 and 9,000 cal y BP, a peak rate of 0.028 cm\u0026middot;yr\u003csup\u003e-\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e\u003c/sup\u003e was recorded at 9,400 cal y BP. Between 9,000 and 4,000 cal y BP, the accumulation rate fluctuated between 0.008 and 0.005 cm yr\u003csup\u003e-\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e\u003c/sup\u003e. The highest accumulation rate occurred from 4,000 cal y BP to the present, reaching approximately 0.02 cm yr\u003csup\u003e-\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e\u003c/sup\u003e (Fig. S3).\u003c/p\u003e\u003c/div\u003e\n\u003ch3\u003ePollen compositional patterns\u003c/h3\u003e\n\u003cp\u003eA total of 91 palynological morphotypes were identified across the 187 samples, encompassing pollen types from 59 genera and 50 families. However, only 46 taxa were present in three or more samples and were selected to run the analyses. Of these, only \u003cem\u003eQuercus humboldtii\u003c/em\u003e and Poaceae were consistently present in all samples throughout the entire record. The mean number of taxa per sample was of 24 \u0026plusmn; 4. The NDMS analysis based on the 46 morphotypes yielded a convergent solution (stress\u0026thinsp;=\u0026thinsp;0.05) after 20 iterations. The NMDS I ordination axis showed a temporal gradient ranging from a modern forest-paramo ecosystem to a forest dominated phase in the middle of the record, followed by a forest/paramo transition. The NMDS II showed a trend of change from forest to paramo from the bottom to the top of the ordination diagram (y-axis), respectively (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e), but emphasizing temporal changes of arboreal composition through time.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eBased on the time-constrained cluster classification applied on the NMDS ordination, we identified five compositional indicator groups of pollen taxa (Table S2). The two most distinct compositional groups were the first and the fifth, representing compositionally different paramo ecosystems located at opposites extremes of the first NMDS axis. The older paramo ecosystem spans a period from 14,700 to 11,000 cal y BP and is characterized by indicator taxa typical of herbaceous communities, such as \u003cem\u003eValeriana, Acaena, Puya\u003c/em\u003e and \u003cem\u003eVallea\u003c/em\u003e. In contrast, the modern paramo (extending from 2,300 cal y BP to the present) shows an increase of some paramo-related taxa, such as Iridiaceae, Rosaceae, Ericaceae and \u003cem\u003eHypericum\u003c/em\u003e, along with a rise in sub-Andean taxa, such as \u003cem\u003eMorella\u003c/em\u003e and \u003cem\u003eClusia.\u003c/em\u003e The remaining three groups, which span the interval from 10,400 to 2,300 cal y BP, represent transitional shifts in forest species composition. Although distinct, these groups shared several of the most abundant forest taxa. During this period, Andean Forest taxa reached a peak in tree pollen abundance. The transition from the preceding zone (the older paramo) is marked by the increasing presence of arboreal taxa such as \u003cem\u003eGordonia\u003c/em\u003e, \u003cem\u003eJuglans, Podocarpus\u003c/em\u003e and \u003cem\u003eDrymis\u003c/em\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\n\u003ch3\u003eUnivariate periodic decomposition\u003c/h3\u003e\n\u003cp\u003eThe wavelet decomposition of the Fe time series, employed as a proxy of precipitation/humidity, showed a maximum average power spectral density in the 3000-year periodicity band between 12,000 and 8,000 cal y BP. During this period, additional regions of high variance were detected in the 2048-year periodicity band. Within the high-resolution sedimentary core, Fe fluctuations during the last glacial interval exhibited oscillations ranging from 512 to 1000 years, all of them encompassed within a broader 2048-year cycle. After 8,000 cal y BP, both variance and periodicity patterns showed a reduction in precipitation compared to earlier periods. Specifically, the power spectrum within the 2048-year periodicity band declined, and the 512-year cycles nearly disappeared, representing the pattern of shift in precipitation regimes from the last glacial maximum to the Holocene (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ea).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eTemperature changes inferred from the pollen record exhibited significant variability throughout the Holocene. Between 12,000 and 8,000 cal y BP, temperatures ranged from 9.4 to 15.4\u0026deg;C, with the wavelet transform showing a statistically significant region dominated by a 2048-year periodicity. From 8,000 to 6,000 cal y BP, the time series exhibited a shift in frequency, with oscillations occurring predominantly at intervals between 1,000 and 500 years, and temperatures fluctuating between 12.7 and 14.2\u0026deg;C. Between approximately 6,000 to 2,500 cal y BP, temperatures peaked at 15.7\u0026deg;C. However, during this interval, the power spectra corresponding to the 500 and 1000-year bands disappeared from the spectrogram, indicating a loss of regular periodicity in temperature variation (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eb).\u003c/p\u003e\u003cp\u003eThe wavelet power spectra applied to the first two NMDS axes revealed a dominant 2,048-year periodicity throughout the entire time series for both axes. However, distinct patterns emerged at shorter timescales. The first axis (NMDS I) exhibited significant spectral power in the 1,000 -512-year band during two intervals: from 10,000 to 8,000 cal y BP, and from 4,000 cal y BP to the present (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ec). In contrast, NMDS II showed peak variance at the 2,048 years between 10,000 and 4,000 cal y BP, along with a concurrent 512-year cycle appearing between 10,000\u0026ndash;8,000 cal y BP (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e5\u003c/span\u003ed). Notably, NMDS II axis displayed statistically significant high-power spectra between 10,000\u0026ndash;8,000 cal y BP, which aligned with periodicities in the abundance of taxa such as Euphorbiaceae, \u003cem\u003eJuglans, Podocarpus\u003c/em\u003e, and \u003cem\u003eMorella\u003c/em\u003e. This pattern suggests a potential turnover in arboreal community composition during this period (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ed).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\n\u003ch3\u003eEnvironmental and pollen compositional covariation\u003c/h3\u003e\n\u003cp\u003eThe wavelet coherence analysis between Fe concentrations and NMDS I showed a strong statistical association at shorter periodicities (~\u0026thinsp;512 years), particularly between 6,500 and 8,000 cal y BP (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e4\u003c/span\u003ea). Additionally, a statistically significant region of high wavelet coherence between Fe and NMDS II time series was observed between 6,000 and 8,500 cal y BP, spanning the 2000 to 512-year periodicity bands (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e4\u003c/span\u003eb). In both instances, a temporal lag was evident between Fe and the ordination axis, implying that fluctuation in Fe levels is a likely driver of plant community dynamics. However, between 6,000 to 8,000 cal y BP, NMDS II and Fe were out of phase with a lag of -π/2 or \u0026frac14; of a quasi-period within the 1000- and 512-year bands with changes in Fe associated with compositional variation following a lag of approximately one quarter of quasi-period (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e5\u003c/span\u003eb).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eRegarding the relationship between temperature (MAT) and NMDS I, we identified a statistically significant coherence region between 14,000 and 10,000 cal y BP, with dominant cycles at approximately 512 years and 1000 years (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e4\u003c/span\u003ec). Similarly, a broader band of high coherence values was found between NMDS II and temperature over the same interval, with a predominant periodicity of 500 years (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e4\u003c/span\u003ed). During most of the Holocene (from 8,000 to 2,000 cal y BP), we also detected a high coherence region with a predominant cycle of 1000 years approximately. The phase difference shows that, for most of the record, the two series were out of phase by π/2. This indicates that changes in temperature were generally followed by shifts in plant community composition (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e5\u003c/span\u003ec, d).\u003c/p\u003e\u003cp\u003eThe influence of either MAT or Fe as drivers of indicator taxa, assessed through wavelet coherence analysis on a 500-year filtered series, showed different taxa-specific responses to each environmental variable (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003e). For example, taxa such as \u003cem\u003eValeriana\u003c/em\u003e and Cariophyllaceae, representative of a paramo flora community that shifted around 11,000 cal y BP, showed a positive relationship with increased humidity (as indicated by Fe concentrations). Likewise, the peaks of abundance of taxa such as \u003cem\u003eJuglans\u003c/em\u003e, \u003cem\u003ePodocarpus\u003c/em\u003e and \u003cem\u003eMorella\u003c/em\u003e mirrored the high frequency patterns in the Fe series, explaining the strong coherence found on wavelet decomposition of both series between 8,500 and 6,000 cal y BP (Figure S4). Regarding temperature, indicator taxa such as \u003cem\u003eGordonia\u003c/em\u003e, \u003cem\u003eAragoa\u003c/em\u003e and \u003cem\u003eJuglans\u003c/em\u003e tracked the main fluctuation patterns of the MAT series, consistent with the observed wavelet coherence results (Fig. S5).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e"},{"header":"Discussion","content":"\u003cp\u003eOur study aligns with previous studies in which historical climate variability during the Holocene operated as a key driver of elevational shifts in the boundary between upper montane forests (UMF) and paramos in the northern Andes\u003csup\u003e\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e\u003c/sup\u003e. Using a quantitative approach, such as wavelet coherence, we were able to statistically evaluate the covariation between time series associated with environmental drivers and both plant communities and individual indicator species over the past 15,000 cal y BP in the northern Andes. We found that both temperature and humidity played a key role in shaping the historical transitions between forests to paramo ecosystems, as well as on driving species turnover through time. However, in contrast to former studies carried out in the Eastern cordillera of the Colombian Andes, our findings suggest that, in the Central cordillera of the Colombian Andes, climatic periods as the Bølling–Allerød and the Younger Dryas were characterized by wet climate\u003csup\u003e\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e\u003c/sup\u003e, rather than dry as reported in the eastern flank\u003csup\u003e\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eOur results showed a pattern of ecosystem change that shifted from a paramo dominated by Poaceae, Asteraceae and Puya, to a tree-dominated community characterized by \u003cem\u003eMorella sp, Podocarpus sp\u003c/em\u003e, and high abundance of large trees such as \u003cem\u003eQ. humboldtii, Juglans sp, Weinmannia sp\u003c/em\u003e and \u003cem\u003eAlchornea sp\u003c/em\u003e, many of which decreased in abundance or disappeared from the palynological record to let the system to become a forest/paramo transitional ecosystem again (see Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e). Similar transitions and changes in the position of the upper forest line (UFL) have been reported in previous studies in northwest Andes \u003csup\u003e\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e,\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e,\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e,\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e\u003c/sup\u003e, and have been primarily attributed to temperature changes associated with major climatic oscillations, such as the Last Glacial Maximum and the Younger Dryas \u003csup\u003e\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e,\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e\u003c/sup\u003e. Our findings, however, emphasize a differential response of plant taxa to either temperature or precipitation, improving our understanding of the likely response of this ecosystem to droughts and/or temperature increases.\u003c/p\u003e\u003cp\u003eAs evidenced along the second NMDS axis, distinct compositional changes also occurred within the forest dominated system during periods with relatively stable temperature. For example, here we show how \u003cem\u003eAlnus\u003c/em\u003e positively responded to changes in precipitation, while abundance of \u003cem\u003eQuercus humboldtii\u003c/em\u003e rose along with the increase of temperature (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e5\u003c/span\u003e). Yet the potential influence of unmeasured biological mechanisms on shaping tree turnover during the forest-dominated phase, such as negative density dependence \u003csup\u003e\u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e\u003c/sup\u003e or ecological drift \u003csup\u003e\u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e40\u003c/span\u003e,\u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e\u003c/sup\u003e remains untested. To more effectively unravel the likely effect of biotic interactions and ecological drift on long-term plant species turnover in relatively stable biomes, we would require a higher taxonomic resolution than that employed here and in most palaeoecological studies (i.e., mostly genus or family), which seems a difficult technical challenge to overcome.\u003c/p\u003e\u003cp\u003eOur time-series decomposition of the pollen community in the Belmira paramo revealed huge variance in both environmental and compositional variation signals during the late Holocene. According to the wavelet spectrum analysis (see Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e), a significant and pronounced shift in plant composition occurred around 10,000 cal y BP, as reflected in both the transitional-paramo axis (NMDS-I) and the forest-paramo axis (NMDS-II). While previous studies have linked changes in the upper forest line of the northern Andes to long-term astronomical orbital cycles of 100,000-year (eccentricity) and 41,000-year (obliquity)\u003csup\u003e\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e\u003c/sup\u003e, our findings support the existence of internal cyclical variations at shorter temporal scales within interstadial periods of the major climatic oscillations\u003csup\u003e\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e,\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eIndependent univariate wavelet analyses showed similar periodicities between 500 and 1000 years for MAT and the orthogonal compositional NMDS axes, while Fe (used as a proxy for humidity) exhibited significant variance at longer periodicities, ranging from 2000 to 4000 years. Using pollen data to reconstruct past temperature assumes that species distributions remained in equilibrium with environmental conditions over time. This approach overlooks species dispersal limitations\u003csup\u003e\u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e42\u003c/span\u003e\u003c/sup\u003e and assumes a stable adiabatic rate throughout the entire study period, assumptions that may not hold across millennial timescales\u003csup\u003e\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e\u003c/sup\u003e. We definitively need more studies with finer temporal pollen resolution as well as the use of complementary proxies highly sensitive to temperature, such as chironomid assemblages\u003csup\u003e\u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e43\u003c/span\u003e\u003c/sup\u003e, to enhance our understanding of the extent to which subtle changes in climate dynamics can determine species turnover at different temporal and spatial scales.\u003c/p\u003e\u003cp\u003eOne of the main goals of our study was to use wavelet coherence (WCO) to provide a quantitative objective interpretation of covariation based on time series decomposition\u003csup\u003e\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e\u003c/sup\u003e, rather than a visual comparison between abundance of pollen and climatic proxies, as it has been done in most traditional palynological studies. At the whole pollen community level, our results revealed spatially synchronous processes (i.e., temporal co-evolution) at regional scales between changes in pollen composition (i.e., NMDS axes I and II) and environmental variables such as MAT and Fe concentrations (proxy of humidity). In the case of the NMDS I, which represented the major changes from one biome to another (paramo-forest-paramo), we observed significant synchrony with MAT at periodicities of 500–1000-year period mainly around 10000–12,000 y cal BP. This corresponds to the periods during which landscape transitioned from a paramo to a forested dominated ecosystem, implying an upward shift of the upper forest line in response to rising temperatures. That said, forest cover, and thus the UFL, expanded upslope in synchrony with warming. In contrast, humidity showed a stronger synchrony with NMDS II, especially at shorter periodicities (i.e., ≤ 500 years) between 6,000–8,000 y cal BP, a time in which the region was likely fully covered by forest. However, our study points to the existence of a lag in vegetation tracking to climate change of even centuries\u003csup\u003e\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e\u003c/sup\u003e, which constitutes a major challenge for the implementation of effective nature-based strategies of conservation to ameliorate the ongoing global change\u003csup\u003e\u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e42\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eThe temporal correlation between NMDS I and Fe revealed a significant synchrony at approximately 500-year periodicities, particularly between 6,000–8,000 y cal BP. This pattern closely matched the covariation observed between NMDS I and MAT, suggesting that Fe and MAT interacted synchronously to influence tree turnover within the forest dominated biome. The strong and significant covariation between NMDS II and Fe further confirmed the role of humidity in shaping tree species turnover but expanded its influence on longer cycles of approximately 1,000–2,000 years. This highlights humidity as a key driver of long-term compositional change in Andean Forest communities. Although pronounced changes in humidity might typically be expected to trigger biome-level transitions, as observed in lowland savannas (e.g., from a forest to a grassland), our findings suggest that precipitation variability in the study region did not play a role in determining ecosystemic transitions at the biome level as strong as temperature variability. This supports the view that significant changes in temperature play a dominant role in driving systematic species turnover and upslope migration in tropical mountain ecosystems\u003csup\u003e\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e,\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eOne of the most striking contributions of this study lies in its methodological approach, which analyzes pollen community-level using common tool employed in modern community ecology. Specifically, we highlight: 1) the selection of ecological groups and indicator species using statistical methods rather than visual or subjective criteria; and 2) the application of a quantitative approach, wavelet coherence analysis, to assess statistical associations between environmental drivers and the pollen abundance of indicator taxa, instead of simple visual interpretation. This approach allowed us to conclude that the abundance of some taxa in these forests, such as \u003cem\u003eDrymis sp\u003c/em\u003e and \u003cem\u003eQuercus humboldtii\u003c/em\u003e, was primarily associated with temperature, while changes in the abundance of other taxa, such as Cyatheaceae, \u003cem\u003ePodocarpus sp.\u003c/em\u003e, and \u003cem\u003ePuya sp.\u003c/em\u003e, were strongly associated with precipitation (Figure S4). Some taxa, such as \u003cem\u003eMelpomene\u003c/em\u003e, showed no clear association with either temperature or precipitation. Overall, applying wavelet coherence and power spectral decomposition to investigate the effects of past climate change on species abundance\u003csup\u003e\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e\u003c/sup\u003e offers a powerful tool for inferring how current plant populations may respond to the ongoing climate change.\u003c/p\u003e\u003cp\u003eAlthough we do not present a formal statistical test to assess the influence of major climatic oscillations on the observed changes in pollen abundance and composition due to the lack of data to do so, our results allow us to identify some of the likely effects that these climatic events may have had on shaping vegetation dynamics. The increased abundance of taxa such as Poaceae, Asteraceae, \u003cem\u003eAcaena\u003c/em\u003e, and \u003cem\u003ePuya\u003c/em\u003e, along with a lowering of the UFL during the late glacial interval, mirrors vegetation trends previously reported in the Central and Eastern Andes during cold periods \u003csup\u003e\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e\u003c/sup\u003e. Between ~ 15,000–13,000 cal y BP, the Fe time series exhibited a region of significant high-power spectra that may correspond to the Guantiva interstadial (~ 14,200–11,200 cal y BP) and the Abra stadial events (~ 11,200–9,700 cal y BP), which are considered the South American analogs of the Boiling-Allerod and Younger-Dryas, respectively\u003csup\u003e\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e\u003c/sup\u003e. In our study area, a subtle reduction in the UFL seems to be associated with the Abra stadial period. This interpretation is supported by the downward shift of the UFL of around ~ 300 masl below the level observed at the end of the Guantiva interstadial. This shift likely resulted from a decrease in rainfall, accompanied by a modest increase in the abundance of cold-adapted taxa such as Poaceae and Asteraceae.\u003c/p\u003e\u003cp\u003eThe findings of this research corroborate the notion that YD event was characterized by wet conditions in the Colombian Andes, consistent with evidenced from other palynological data\u003csup\u003e\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e,\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e,\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e\u003c/sup\u003e. During the BA period, the patterns of precipitation were influenced by a northward shift of the Intertropical Convergence Zone (ITCZ); this shift may have increased moisture influx from the Amazon basin contributing to the prevailing wet conditions during the YD. Around 10,000 cal yr BP, a transition to drier conditions is indicated by precipitation proxies. This climatic shift has been attributed to the weakening of the South American Summer Monsoon (SASM), which reduced moisture transport from the Amazon basin to the Andean region\u003csup\u003e\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e,\u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e44\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eAlthough our dataset is not extensive enough to directly assess the effect of recent anthropogenic warming, the long-term trend observed over the past 15,000 years suggests that continued temperature increase must lead to a renewed shift toward a forest dominated ecosystem. However, this expectation of an upslope migration of the upper forest line in response to rising temperatures contrasts with the transitional trend observed over the past 2,500 years, which indicates a reversal from forest to paramo. This suggests that ongoing human-induced climate change may override the natural ecological trajectory of recent millennia, potentially reversing the recent expansion of páramo ecosystems. Such a shift could have direct consequences for the provision of key ecosystem services, including water supply and regulation, which are critically linked to paramo vegetation.\u003c/p\u003e"},{"header":"Methods","content":"\u003ch2\u003eStudy area\u003c/h2\u003e\u003cp\u003eThe sediment core was collected from El Morro mire (06°40´11.9˝N, 75°40´4.4˝W), part of Belmira (Antioquia, Colombia) paramo in the northern Andes. The local bedrock primarily consists of quartz diorite formed during the late Cretaceous period, alongside a metamorphic complex that dates from the Precambrian era to the middle Cretaceous period. The regional low fertility and acidic soils are characterized by high organic matter (OM) and aluminum content. The paramo spans an altitudinal range between 3000 and 3400 masl, encompassing steep slopes and a flat zone characterized by diverse glacial landforms. The area experiences a tropical alpine climate, with an average annual temperature of 12°C and diurnal temperature fluctuations of approximately 10°C. The annual precipitation in this region reaches around 2710 mm, and is closely influenced by mesoscale processes, particularly the movement of the Inter Tropical Convergence Zone (ITCZ). Consequently, two distinct rainy periods occur: one from April to May, when the ITCZ moves northward, and another from September to October, when the ITCZ shifts southwards\u003csup\u003e\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e\u003c/sup\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e).\u003c/p\u003e\u003ch3\u003eData collection\u003c/h3\u003e\u003cp\u003eAt El Morro mire, we recovered a 300 cm long sediment core using a hand-operated Russian sampler, which extracts 50cm long, half-cylindrical sections with a 5cm diameter. A first stratigraphical description of the sediments was photographed in the field to document relevant core data, such as color and texture. The cores were wrapped in plastic foils and protected by PVC guttering for transport. Samples were taken to the paleoecology Laboratory of the National University in Medellín, where they were stored in a dark room at 4°C. A total of 300 samples were extracted from the sediment core, maintaining a one-centimeter distance between the closest samples. These samples were treated with standard procedures for pollen preparation, including sodium pyrophosphate and acetolysis. Then, the pollen samples were prepared with glycerin gelatin and fixed in microscope slides with wax. Using pollen morphological studies and the pollen reference collection of the Paleoecology Laboratory of the National University of Colombia (Medellín)\u003csup\u003e\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e\u003c/sup\u003e, pollen grains of each sample were counted and identified to the highest possible taxonomic resolution (i.e., species, genus, or family).\u003c/p\u003e\u003cp\u003eLithological and stratigraphic characteristics were assessed through visual inspection of the core, recording details such as grain size, color, and strata contacts. A total of 17 samples were dated at Beta Analytic Laboratory (Miami). However, only eight samples (approximately located between 60 and 190 cm) were suitable for use in the age model. The remaining samples exhibited high uncertainties in their estimated reversal ages and were therefore excluded (Fig. S6). We used the \u003cem\u003eRBacon\u003c/em\u003e R package\u003csup\u003e\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e\u003c/sup\u003e to calibrate the AMS radiocarbon (\u003csup\u003e\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e\u003c/sup\u003eC) ages, estimate sediment accumulation rates, and to construct the age-depth model (Table \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eTo estimate the age of the pollen samples, we combined radiocarbon dating and depth measurements. First, the radiocarbon dates provided by the laboratory were calibrated using the IntCal20 calibration curve, which is appropriate for samples from the North Hemisphere. Next, we established a relationship between the radiocarbon-based dates and their corresponding depths using a Bayesian autoregressive gamma process model. This model accounted for the autocorrelation among samples resulting from sedimentary processes. Furthermore, our model accounted for uncertainties associated with the radiocarbon dates, including age uncertainty and \u003csup\u003e\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e\u003c/sup\u003eC calibration curve uncertainty. It also incorporated sedimentation rates to enhance the accuracy of age estimation\u003csup\u003e\u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e45\u003c/span\u003e,\u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e46\u003c/span\u003e\u003c/sup\u003e. This analysis was performed using the \u003cem\u003eRbacon\u003c/em\u003e package in R\u003csup\u003e\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eElemental analysis was performed on the core using the EDAX Eagle III XPL mProbe at the Earth Science Department of the University of Geneva. Material for X-ray microfluorescence (mXRF) analyses were sampled at the National University of Colombia in Medellín in 50 cm-long U-channels from the cores recovered with the Russian corer and transported to the University of Geneva. The U-channels were cut into 10 cm long sections to be analyzed in the Eagle III spectrometer. The parameters employed for these measurements were set at 40 kV, 450 mA, 50/55 Dtm at atmospheric pressure. Measured values are expressed in counts per second (cps). Titanium (\u003cem\u003eTi\u003c/em\u003e) and Iron (\u003cem\u003eFe\u003c/em\u003e) concentrations, which represent watershed sedimentary and deposit trajectories due to their relationship with rainfall and runoff processes (i.e., precipitation proxies) were used as chemical proxies to assess the effect of climatic variation on pollen composition during the late Quaternary\u003csup\u003e\u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e47\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003ch2\u003eUpper forest limit definition\u003c/h2\u003e\u003cp\u003eUsing the pollen data, we built a community matrix in which rows represent each dated pollen sample (i.e., time), and columns correspond to the identified taxa (e.g., species, genus, or family). We classified plant taxa by growth form as tree, shrub, or herb, and use the percentage of arboreal pollen (%AP) to define the elevational position of the Upper Forest Limit (UFL)\u003csup\u003e\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e\u003c/sup\u003e (Flantua et al., 2019). We assume that each 5% change in %AP corresponds to a 100 m shift in the UFL. The equation employed to estimate the UFL was as follows\u003csup\u003e\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e,\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u003c/sup\u003e:\u003c/p\u003e\u003cdiv id=\"Equa\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equa\" name=\"EquationSource\"\u003e\n$$\\:UFL=\\:ETL+\\left(\\frac{\\%AP-C}{R}\\right)\\times\\:\\:100$$\u003c/div\u003e\u003c/div\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eWhere ETL is the current elevation of the tree line (3076m asl) for Belmira Paramo based on a Paramo delimitation project; C is the current percentage of arboreal pollen in the study area; and R represents the rate of change of the UFL\u003csup\u003e\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e,\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eThe variation of mean annual temperature (MAT) is related to the UFL as follows\u003csup\u003e\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e,\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u003c/sup\u003e:\u003c/p\u003e\u003cdiv id=\"Equb\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equb\" name=\"EquationSource\"\u003e\n$$\\:MAT=\\:CMT+\\left(UFL-ETL\\right)\\times\\:\\:AR$$\u003c/div\u003e\u003c/div\u003e\u003cp\u003eWhere CMT represents the current elevation range of the Belmira Paramo and AR denotes the adiabatic rate. In our study, we used an upward displacement in elevation. This relationship indicates that the mean annual temperature tends to decrease with increasing elevation, with the Upper Forest Line serving as a critical point where significant temperature changes occur.\u003c/p\u003e\u003ch2\u003ePatterns of pollen compositional change\u003c/h2\u003e\u003cp\u003eWe used Bray-Curtis distance to calculate pollen compositional dissimilarity between samples. To estimate this metric, we included only those taxa present in three or more samples and applied a Non-Metric Dimensional Scaling (NMDS) ordination to reduce the pollen community matrix into orthogonal axes. Taxa found in only one or two samples were excluded to: a) reduce statistical noise due to under sampling; b) identify groups of compositional change based on well represented taxa\u003csup\u003e\u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e48\u003c/span\u003e\u003c/sup\u003e; and c) detect temporal periodicity to be incorporated into the wavelet analysis (see below). We used the first two NMDS axes (NMDS I and NMDS II) to assess plant community over the past 15,000 years BP. Next, we applied a time-constrained CONISS (Constrained Incremental Sum of Squares) clustering algorithm\u003csup\u003e\u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e49\u003c/span\u003e\u003c/sup\u003e to identify homogeneous floristic groups. Once these groups were defined, we performed a species indicator analysis\u003csup\u003e\u003cspan citationid=\"CR50\" class=\"CitationRef\"\u003e50\u003c/span\u003e,\u003cspan citationid=\"CR51\" class=\"CitationRef\"\u003e51\u003c/span\u003e\u003c/sup\u003e to identify the most representative taxa within each plant community. This analysis was conducted using the R packages \u003cem\u003evegan\u003c/em\u003e\u003csup\u003e\u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e52\u003c/span\u003e\u003c/sup\u003e and \u003cem\u003eindicspecies\u003c/em\u003e\u003csup\u003e\u003cspan citationid=\"CR51\" class=\"CitationRef\"\u003e51\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003ch2\u003eUnivariate changes of pollen composition and environmental drivers\u003c/h2\u003e\u003cp\u003eWe applied the continuous wavelet transform\u003csup\u003e\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e,\u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e\u003c/sup\u003e to detect periodic changes in plant composition, wet/dry periods, cycles, and seasonality in time series spanning the last 15,000 y. Specifically, we computed the convolution of each time series using the Morlet mother wavelet and estimated the mean power spectrum by averaging the wavelet coefficients within each frequency component. The Morlet mother wavelet is defined as follows:\u003c/p\u003e\u003cdiv id=\"Equc\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equc\" name=\"EquationSource\"\u003e\n$$\\:{\\psi\\:}_{0}\\left(\\eta\\:\\right)=\\:{\\pi\\:}^{1/4}{e}^{-i{\\omega\\:}_{0}\\eta\\:}{e}^{-{\\eta\\:}^{2}/2}$$\u003c/div\u003e\u003c/div\u003e\u003cp\u003eWhere \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:i=\\:\\sqrt{-1}\\)\u003c/span\u003e\u003c/span\u003e, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\omega\\:}_{0}\\:\\)\u003c/span\u003e\u003c/span\u003eis ω\u003csub\u003e0\u003c/sub\u003e dimensionless frequency, and η is dimensionless time. The Morlet wavelet is particularly useful to analyze periodicities in time series. We assumed the standard value of ω\u003csub\u003e0\u003c/sub\u003e = 6\u003csup\u003e37\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eWe applied univariate wavelet transforms to the Fe, MAT, NMDS I and NMDS II time series to analyze temporal periodicities in precipitation, temperature, and pollen community composition. To assess the statistical significance of regions with high power spectra, we used a null model that preserves the short-term autocorrelation observed in the original time series. This approach allows us to evaluate the significance of long-term patterns in the data\u003csup\u003e\u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e\u003c/sup\u003e. These analyses were performed using the set of MATLAB® functions provided by Cazelles et al (2008). Additionally, based on the wavelet transform of each taxon`s abundance time series, we applied hierarchical clustering to the resulting frequency decomposition matrix to identify taxa with similar temporal periodicities.\u003c/p\u003e\u003ch2\u003ePatterns of covariation using wavelet coherence analysis\u003c/h2\u003e\u003cp\u003eWe used wavelet coherence analysis (WCO) to assess the paired variation between environmental proxies (Fe and MAT) and pollen compositional (NMDS I and NMDS II) time series. We also analyzed WCO between the environmental proxies and the pollen abundance of those taxa identified as indicators for each independent group defined by the cluster analysis, to evaluate the influence of temperature and/or precipitation on shaping the signal of each indicator taxon in the palynological record. This analysis allows us to identify regions of convergence in both time and frequency, where the power spectra of the two series exhibit a linear correlation. Furthermore, it quantifies the phase difference between the two series, providing information about potential time lags\u003csup\u003e\u003cspan citationid=\"CR53\" class=\"CitationRef\"\u003e53\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eTo statistically evaluate the significance of coherence, we compared the observed cross-wavelet spectra with those generated by a null model cross-wavelet spectrum. The null model was created by simulating independent red-noise processes with the same first order autoregressive coefficient as the original time series. This comparison allows us to determine whether the observed cross-wavelet spectra significantly deviated from expectations under the null hypothesis\u003csup\u003e\u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e37\u003c/span\u003e\u003c/sup\u003e. To enhance the interpretation of the correlation between environmental and compositional variation variables, we employed a low-pass Gaussian filter to both the proxy data (MAT and Fe) and the abundance series of indicator taxa. This filter was calibrated using a bandpass value derived from regions of high coherence, such as the 500-year filtered series. These analyses were also performed using the MATLAB® functions provided by Cazelles et al (2008).\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003ch2\u003eAdditional information\u003c/h2\u003e\u003cp\u003eCompeting financial interests: The authors declare no competing financial interests.\u003c/p\u003e\u003c/p\u003e\u003ch2\u003eFunding\u003c/h2\u003e\u003cp\u003eThe authors received no funding for this work\u003c/p\u003e\u003ch2\u003eAuthor Contribution\u003c/h2\u003e\u003cp\u003eA.S-D., A.D., and C.V. conceived and designed the experiments. A.S-D, I.M.C-R, and R.E.V-M performed the experiments. A.S-D analyzed the data with advice from S.G-C. and A.D. A.S-D. and A.D. wrote the main text. All authors reviewed and commented on the final version of the manuscript.\u003c/p\u003e\u003ch2\u003eAcknowledgement\u003c/h2\u003e\u003cp\u003eWe thank Professor Bernard Cazelles for providing the code to run univariate wavelet decomposition, wavelet coherence, and phase difference analysis.\u003c/p\u003e\u003ch2\u003eData Availability\u003c/h2\u003e\u003cp\u003eThe datasets used and/or analysed during the current study available from the corresponding author on reasonable request.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eK\u0026ouml;rner, C. \u003cem\u003eAlpine Treelines\u003c/em\u003e. (Springer, Basel, 2012). doi:10.1007/978-3-0348-0396-0.\u003c/li\u003e\n\u003cli\u003eFlantua, S. G. A., O\u0026rsquo;Dea, A., Onstein, R. E., Giraldo, C. \u0026amp; Hooghiemstra, H. The flickering connectivity system of the north Andean p\u0026aacute;ramos. \u003cem\u003eJournal of Biogeography\u003c/em\u003e \u003cstrong\u003e46\u003c/strong\u003e, 1808\u0026ndash;1825 (2019).\u003c/li\u003e\n\u003cli\u003eGroot, M. H. M. \u003cem\u003eet al.\u003c/em\u003e Ultra-high resolution pollen record from the northern Andes reveals rapid shifts in montane climates within the last two glacial cycles. \u003cem\u003eClimate of the Past\u003c/em\u003e \u003cstrong\u003e7\u003c/strong\u003e, 299\u0026ndash;316 (2011).\u003c/li\u003e\n\u003cli\u003eAlley, R. B. Ice-core evidence of abrupt climate changes. \u003cem\u003eProc. Natl. Acad. Sci. U.S.A.\u003c/em\u003e \u003cstrong\u003e97\u003c/strong\u003e, 1331\u0026ndash;1334 (2000).\u003c/li\u003e\n\u003cli\u003eMu\u0026ntilde;oz, P. \u003cem\u003eet al.\u003c/em\u003e Holocene climatic variations in the Western Cordillera of Colombia: A multiproxy high-resolution record unravels the dual influence of ENSO and ITCZ. \u003cem\u003eQuaternary Science Reviews\u003c/em\u003e \u003cstrong\u003e155\u003c/strong\u003e, 159\u0026ndash;178 (2017).\u003c/li\u003e\n\u003cli\u003eVel\u0026aacute;squez-R., C. A. \u0026amp; Hooghiemstra, H. Pollen-based 17-kyr forest dynamics and climate change from the Western Cordillera of Colombia; no-analogue associations and temporarily lost biomes. \u003cem\u003eReview of Palaeobotany and Palynology\u003c/em\u003e \u003cstrong\u003e194\u003c/strong\u003e, 38\u0026ndash;49 (2013).\u003c/li\u003e\n\u003cli\u003eHooghiemstra, H., Wijninga, V. M. \u0026amp; Cleef, A. M. THE PALEOBOTANICAL RECORD OF COLOMBIA: IMPLICATIONS FOR BIOGEOGRAPHY AND BIODIVERSITY1. \u003cem\u003emobt\u003c/em\u003e \u003cstrong\u003e93\u003c/strong\u003e, 297\u0026ndash;325 (2006).\u003c/li\u003e\n\u003cli\u003eMarchant, R. \u003cem\u003eet al.\u003c/em\u003e Colombian vegetation at the Last Glacial Maximum: a comparison of model- and pollen-based biome reconstructions. \u003cem\u003eJournal of Quaternary Science\u003c/em\u003e \u003cstrong\u003e19\u003c/strong\u003e, 721\u0026ndash;732 (2004).\u003c/li\u003e\n\u003cli\u003eBader, M. Y. \u003cem\u003eet al.\u003c/em\u003e A global framework for linking alpine-treeline ecotone patterns to underlying processes. \u003cem\u003eEcography\u003c/em\u003e \u003cstrong\u003e44\u003c/strong\u003e, 265\u0026ndash;292 (2021).\u003c/li\u003e\n\u003cli\u003eBush, M. B. Distributional change and conservation on the Andean flank: a palaeoecological perspective. \u003cem\u003eGlobal Ecology and Biogeography\u003c/em\u003e \u003cstrong\u003e11\u003c/strong\u003e, 463\u0026ndash;473 (2002).\u003c/li\u003e\n\u003cli\u003eMart\u0026iacute;nez, C. \u003cem\u003eet al.\u003c/em\u003e Neogene precipitation, vegetation, and elevation history of the Central Andean Plateau. \u003cem\u003eSci. Adv.\u003c/em\u003e \u003cstrong\u003e6\u003c/strong\u003e, eaaz4724 (2020).\u003c/li\u003e\n\u003cli\u003ePati\u0026ntilde;o, L. \u003cem\u003eet al.\u003c/em\u003e Late Pleistocene\u0026ndash;Holocene environmental and climatic history of a freshwater paramo ecosystem in the northern Andes. \u003cem\u003eJournal of Quaternary Science\u003c/em\u003e \u003cstrong\u003e35\u003c/strong\u003e, 1046\u0026ndash;1056 (2020).\u003c/li\u003e\n\u003cli\u003eBuytaert, W. \u0026amp; De Bi\u0026egrave;vre, B. Water for cities: The impact of climate change and demographic growth in the tropical Andes. \u003cem\u003eWater Resources Research\u003c/em\u003e \u003cstrong\u003e48\u003c/strong\u003e, (2012).\u003c/li\u003e\n\u003cli\u003eImmerzeel, W. W. \u003cem\u003eet al.\u003c/em\u003e Importance and vulnerability of the world\u0026rsquo;s water towers. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e577\u003c/strong\u003e, 364\u0026ndash;369 (2020).\u003c/li\u003e\n\u003cli\u003eDuque, A. \u003cem\u003eet al.\u003c/em\u003e Mature Andean forests as globally important carbon sinks and future carbon refuges. \u003cem\u003eNat Commun\u003c/em\u003e \u003cstrong\u003e12\u003c/strong\u003e, 2138 (2021).\u003c/li\u003e\n\u003cli\u003ePe\u0026ntilde;a, M. A., Feeley, K. J. \u0026amp; Duque, A. Effects of endogenous and exogenous processes on aboveground biomass stocks and dynamics in Andean forests. \u003cem\u003ePlant Ecology\u003c/em\u003e \u003cstrong\u003e219\u003c/strong\u003e, 1481\u0026ndash;1492 (2018).\u003c/li\u003e\n\u003cli\u003eSegovia, R. A. \u003cem\u003eet al.\u003c/em\u003e Freezing and water availability structure the evolutionary diversity of trees across the Americas. \u003cem\u003eSci. Adv.\u003c/em\u003e \u003cstrong\u003e6\u003c/strong\u003e, eaaz5373 (2020).\u003c/li\u003e\n\u003cli\u003eGonz\u0026aacute;lez-Caro, S. \u003cem\u003eet al.\u003c/em\u003e The legacy of biogeographic history on the composition and structure of Andean forests. \u003cem\u003eEcology\u003c/em\u003e \u003cstrong\u003e101\u003c/strong\u003e, e03131 (2020).\u003c/li\u003e\n\u003cli\u003eMadri\u0026ntilde;\u0026aacute;n, S., Cort\u0026eacute;s, A. J. \u0026amp; Richardson, J. E. P\u0026aacute;ramo is the world\u0026rsquo;s fastest evolving and coolest biodiversity hotspot. \u003cem\u003eFront. Genet.\u003c/em\u003e \u003cstrong\u003e4\u003c/strong\u003e, (2013).\u003c/li\u003e\n\u003cli\u003ePouchon, C. \u003cem\u003eet al.\u003c/em\u003e Phylogenomic Analysis of the Explosive Adaptive Radiation of the Espeletia Complex (Asteraceae) in the Tropical Andes. \u003cem\u003eSystematic Biology\u003c/em\u003e \u003cstrong\u003e67\u003c/strong\u003e, 1041\u0026ndash;1060 (2018).\u003c/li\u003e\n\u003cli\u003ePouchon, C. \u003cem\u003eet al.\u003c/em\u003e Phylogenetic signatures of ecological divergence and leapfrog adaptive radiation in Espeletia. \u003cem\u003eAmerican Journal of Botany\u003c/em\u003e \u003cstrong\u003e108\u003c/strong\u003e, 113\u0026ndash;128 (2021).\u003c/li\u003e\n\u003cli\u003eDuque, A., Stevenson, P. R. \u0026amp; Feeley, K. J. Thermophilization of adult and juvenile tree communities in the northern tropical Andes. \u003cem\u003eProceedings of the National Academy of Sciences\u003c/em\u003e \u003cstrong\u003e112\u003c/strong\u003e, 10744\u0026ndash;10749 (2015).\u003c/li\u003e\n\u003cli\u003eFadrique, B. \u003cem\u003eet al.\u003c/em\u003e Widespread but heterogeneous responses of Andean forests to climate change. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e564\u003c/strong\u003e, 207\u0026ndash;212 (2018).\u003c/li\u003e\n\u003cli\u003ePeyre, G. What Does the Future Hold for P\u0026aacute;ramo Plants? A Modelling Approach. \u003cem\u003eFront. Ecol. Evol.\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e, (2022).\u003c/li\u003e\n\u003cli\u003ePeyre, G. \u003cem\u003eet al.\u003c/em\u003e The fate of p\u0026aacute;ramo plant assemblages in the sky islands of the northern Andes. \u003cem\u003eJournal of Vegetation Science\u003c/em\u003e \u003cstrong\u003e31\u003c/strong\u003e, 967\u0026ndash;980 (2020).\u003c/li\u003e\n\u003cli\u003eBuytaert, W. \u003cem\u003eet al.\u003c/em\u003e Human impact on the hydrology of the Andean p\u0026aacute;ramos. \u003cem\u003eEarth-Science Reviews\u003c/em\u003e \u003cstrong\u003e79\u003c/strong\u003e, 53\u0026ndash;72 (2006).\u003c/li\u003e\n\u003cli\u003eVan Der Hammen, T. \u0026amp; Hooghiemstra, H. The EL abra stadial, a younger dryas equivalent in Colombia. \u003cem\u003eQuaternary Science Reviews\u003c/em\u003e \u003cstrong\u003e14\u003c/strong\u003e, 841\u0026ndash;851 (1995).\u003c/li\u003e\n\u003cli\u003eVeer, R. van \u0026rsquo;t, Islebe, G. A. \u0026amp; Hooghiemstra, H. Climatic change during the Younger Dryas chron in northern South America: a test of the evidence. \u003cem\u003eQuaternary Science Reviews\u003c/em\u003e \u003cstrong\u003e19\u003c/strong\u003e, 1821\u0026ndash;1835 (2000).\u003c/li\u003e\n\u003cli\u003eDebret, M. \u003cem\u003eet al.\u003c/em\u003e The origin of the 1500-year climate cycles in Holocene North-Atlantic records. \u003cem\u003eClimate of the Past\u003c/em\u003e \u003cstrong\u003e3\u003c/strong\u003e, 569\u0026ndash;575 (2007).\u003c/li\u003e\n\u003cli\u003eTurner, T. E. \u003cem\u003eet al.\u003c/em\u003e Solar cycles or random processes? Evaluating solar variability in Holocene climate records. \u003cem\u003eSci Rep\u003c/em\u003e \u003cstrong\u003e6\u003c/strong\u003e, 23961 (2016).\u003c/li\u003e\n\u003cli\u003eFastovich, D. \u003cem\u003eet al.\u003c/em\u003e Coupled, decoupled, and abrupt responses of vegetation to climate across timescales. \u003cem\u003eScience\u003c/em\u003e \u003cstrong\u003e389\u003c/strong\u003e, 64\u0026ndash;68 (2025).\u003c/li\u003e\n\u003cli\u003eKnight, C. A. \u003cem\u003eet al.\u003c/em\u003e Community Assembly and Climate Mismatch in Late Quaternary Eastern North American Pollen Assemblages. \u003cem\u003eThe American Naturalist\u003c/em\u003e \u003cstrong\u003e195\u003c/strong\u003e, 166\u0026ndash;180 (2020).\u003c/li\u003e\n\u003cli\u003eCazelles, B. \u003cem\u003eet al.\u003c/em\u003e Wavelet Analysis of Ecological Time Series. \u003cem\u003eOecologia\u003c/em\u003e \u003cstrong\u003e156\u003c/strong\u003e, 287\u0026ndash;304 (2008).\u003c/li\u003e\n\u003cli\u003eGrinsted, A., Moore, J. C. \u0026amp; Jevrejeva, S. Application of the cross wavelet transform and wavelet coherence to geophysical time series. \u003cem\u003eNonlinear Processes in Geophysics\u003c/em\u003e \u003cstrong\u003e11\u003c/strong\u003e, 561\u0026ndash;566 (2004).\u003c/li\u003e\n\u003cli\u003eTorrence, C. \u0026amp; Webster, P. J. Interdecadal Changes in the ENSO\u0026ndash;Monsoon System. (1999).\u003c/li\u003e\n\u003cli\u003eTorres, V., Hooghiemstra, H., Lourens, L. \u0026amp; Tzedakis, P. C. Astronomical tuning of long pollen records reveals the dynamic history of montane biomes and lake levels in the tropical high Andes during the Quaternary. \u003cem\u003eQuaternary Science Reviews\u003c/em\u003e \u003cstrong\u003e63\u003c/strong\u003e, 59\u0026ndash;72 (2013).\u003c/li\u003e\n\u003cli\u003eBeninc\u0026agrave;, E., Ballantine, B., Ellner, S. P. \u0026amp; Huisman, J. Species fluctuations sustained by a cyclic succession at the edge of chaos. \u003cem\u003eProc. Natl. Acad. Sci. U.S.A.\u003c/em\u003e \u003cstrong\u003e112\u003c/strong\u003e, 6389\u0026ndash;6394 (2015).\u003c/li\u003e\n\u003cli\u003eBlaauw, M. \u0026amp; Christen, J. A. Flexible paleoclimate age-depth models using an autoregressive gamma process. \u003cem\u003eBayesian Analysis\u003c/em\u003e \u003cstrong\u003e6\u003c/strong\u003e, 457\u0026ndash;474 (2011).\u003c/li\u003e\n\u003cli\u003eVolkov, I., Banavar, J. R., He, F., Hubbell, S. P. \u0026amp; Maritan, A. Density dependence explains tree species abundance and diversity in tropical forests. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e438\u003c/strong\u003e, 658\u0026ndash;661 (2005).\u003c/li\u003e\n\u003cli\u003eCorrea-Metrio, A., Meave, J. A., Lozano-Garc\u0026iacute;a, S. \u0026amp; Bush, M. B. Environmental determinism and neutrality in vegetation at millennial time scales. \u003cem\u003eJournal of Vegetation Science\u003c/em\u003e \u003cstrong\u003e25\u003c/strong\u003e, 627\u0026ndash;635 (2014).\u003c/li\u003e\n\u003cli\u003eHubbell, S. P. Neutral Theory in Community Ecology and the Hypothesis of Functional Equivalence. \u003cem\u003eFunctional Ecology\u003c/em\u003e \u003cstrong\u003e19\u003c/strong\u003e, 166\u0026ndash;172 (2005).\u003c/li\u003e\n\u003cli\u003eSvenning, J.-C. \u0026amp; Sandel, B. Disequilibrium vegetation dynamics under future climate change. \u003cem\u003eAmerican Journal of Botany\u003c/em\u003e \u003cstrong\u003e100\u003c/strong\u003e, 1266\u0026ndash;1286 (2013).\u003c/li\u003e\n\u003cli\u003eMatthews-Bird, F., Brooks, S. J., Holden, P. B., Montoya, E. \u0026amp; Gosling, W. D. Inferring late-Holocene climate in the Ecuadorian Andes using a chironomid-based temperature inference model. \u003cem\u003eClimate of the Past\u003c/em\u003e \u003cstrong\u003e12\u003c/strong\u003e, 1263\u0026ndash;1280 (2016).\u003c/li\u003e\n\u003cli\u003eRodr\u0026iacute;guez-Zorro, P. A. \u003cem\u003eet al.\u003c/em\u003e Shut down of the South American summer monsoon during the penultimate glacial. \u003cem\u003eSci Rep\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e, 6275 (2020).\u003c/li\u003e\n\u003cli\u003eBlaauw, M. Methods and code for \u0026lsquo;classical\u0026rsquo; age-modelling of radiocarbon sequences. \u003cem\u003eQuaternary Geochronology\u003c/em\u003e \u003cstrong\u003e5\u003c/strong\u003e, 512\u0026ndash;518 (2010).\u003c/li\u003e\n\u003cli\u003eAndr\u0026eacute;s Christen, J. \u0026amp; P\u0026eacute;rez E, S. A New Robust Statistical Model for Radiocarbon Data. \u003cem\u003eRadiocarbon\u003c/em\u003e \u003cstrong\u003e51\u003c/strong\u003e, 1047\u0026ndash;1059 (2009).\u003c/li\u003e\n\u003cli\u003eHaug, G. H., Hughen, K. A., Sigman, D. M., Peterson, L. C. \u0026amp; R\u0026ouml;hl, U. Southward Migration of the Intertropical Convergence Zone through the Holocene. \u003cem\u003eScience\u003c/em\u003e \u003cstrong\u003e293\u003c/strong\u003e, 1304\u0026ndash;1308 (2001).\u003c/li\u003e\n\u003cli\u003eLegendre, P. \u0026amp; Legendre, L. \u003cem\u003eNumerical Ecology\u003c/em\u003e. (Elsevier, 2012).\u003c/li\u003e\n\u003cli\u003eGrimm, E. C. CONISS: a FORTRAN 77 program for stratigraphically constrained cluster analysis by the method of incremental sum of squares. \u003cem\u003eComputers \u0026amp; Geosciences\u003c/em\u003e \u003cstrong\u003e13\u003c/strong\u003e, 13\u0026ndash;35 (1987).\u003c/li\u003e\n\u003cli\u003eC\u0026aacute;ceres, M. D. \u0026amp; Legendre, P. Associations between species and groups of sites: indices and statistical inference. \u003cem\u003eEcology\u003c/em\u003e \u003cstrong\u003e90\u003c/strong\u003e, 3566\u0026ndash;3574 (2009).\u003c/li\u003e\n\u003cli\u003eDe C\u0026aacute;ceres, M., Legendre, P., Wiser, S. K. \u0026amp; Brotons, L. Using species combinations in indicator value analyses. \u003cem\u003eMethods in Ecology and Evolution\u003c/em\u003e \u003cstrong\u003e3\u003c/strong\u003e, 973\u0026ndash;982 (2012).\u003c/li\u003e\n\u003cli\u003eOksanen, J. \u003cem\u003eet al.\u003c/em\u003e vegan: Community Ecology Package. (2025).\u003c/li\u003e\n\u003cli\u003eMaraun, D. \u0026amp; Kurths, J. Cross wavelet analysis: significance testing and pitfalls. \u003cem\u003eNonlinear Processes in Geophysics\u003c/em\u003e \u003cstrong\u003e11\u003c/strong\u003e, 505\u0026ndash;514 (2004).\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"scientific-reports","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"scirep","sideBox":"Learn more about [Scientific Reports](http://www.nature.com/srep/)","snPcode":"","submissionUrl":"","title":"Scientific Reports","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"stoa","reportingPortfolio":"Scientific Reports","inReviewEnabled":true,"inReviewRevisionsEnabled":true},"keywords":"Global climatic events, Paramo, Upper Forest Line, Wavelet coherence, Species turnover","lastPublishedDoi":"10.21203/rs.3.rs-7368410/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-7368410/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eIn this study, conducted in the Belmira paramo in the northwestern Andes of Colombia, we used pollen records to track vegetation changes over the past 15,000 years before present (BP). Employing ordination (NMDS), species indicator, and wavelet analysis, we found a pattern of ecosystem change that shifted from a paramo to a tree-dominated community, which became a forest/paramo transitional ecosystem again. At the whole pollen community level, our results revealed spatially synchronous processes (i.e., temporal co-evolution) at regional scales between changes in pollen composition and mean annual temperature (MAT) and Fe concentrations (proxy of precipitation). Independent univariate wavelet analyses showed similar periodicities between 500 and 1000 years for MAT and the orthogonal compositional NMDS axes, while Fe exhibited significant variance at longer periodicities, ranging from 2000 to 4000 years. The strong and significant covariation assessed using wavelet coherence between NMDS II and Fe further confirmed the role of precipitation in shaping tree species turnover on cycles of approximately 1,000-500 years. However, our study points to the existence of a lag in vegetation tracking to climate change of centuries, which could constitute a major challenge for the implementation of effective nature-based strategies of conservation to ameliorate the ongoing global change.\u003c/p\u003e","manuscriptTitle":"Temperature and Humidity as Drivers of Elevational Shifts in the Montane Forests-Paramo Boundary Over the Last 15,000 Years in the Northwest Andes","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-09-18 10:40:57","doi":"10.21203/rs.3.rs-7368410/v1","editorialEvents":[{"type":"communityComments","content":0},{"type":"decision","content":"Revision requested","date":"2025-11-18T03:37:32+00:00","index":"","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2025-11-17T18:02:03+00:00","index":"hide","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2025-11-14T23:12:08+00:00","index":"hide","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2025-11-11T12:28:25+00:00","index":"hide","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2025-11-10T14:52:33+00:00","index":"hide","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2025-10-13T20:35:01+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"85682283274602756407014370139668308059","date":"2025-09-16T09:53:53+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"150012897360733087455130690643491050508","date":"2025-09-14T21:14:15+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"239144738826968533771905291351435618953","date":"2025-09-14T08:45:59+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"168135240632996028571440324319042209168","date":"2025-09-13T11:33:38+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"207175199962394286184111479070867068755","date":"2025-09-13T08:40:43+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"137420865148769644170480696372843215723","date":"2025-09-12T14:45:07+00:00","index":"hide","fulltext":""},{"type":"reviewersInvited","content":"","date":"2025-09-11T07:55:50+00:00","index":"","fulltext":""},{"type":"editorAssigned","content":"","date":"2025-09-11T07:31:23+00:00","index":"","fulltext":""},{"type":"editorInvited","content":"","date":"2025-09-01T12:24:42+00:00","index":"","fulltext":""},{"type":"checksComplete","content":"","date":"2025-08-21T13:33:42+00:00","index":"","fulltext":""},{"type":"submitted","content":"Scientific Reports","date":"2025-08-21T12:07:31+00:00","index":"","fulltext":""}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"scientific-reports","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"scirep","sideBox":"Learn more about [Scientific Reports](http://www.nature.com/srep/)","snPcode":"","submissionUrl":"","title":"Scientific Reports","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"stoa","reportingPortfolio":"Scientific Reports","inReviewEnabled":true,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"8d55405d-604c-4802-bb92-140b6953d032","owner":[],"postedDate":"September 18th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"under-review","subjectAreas":[{"id":54801007,"name":"Earth and environmental sciences/Climate sciences"},{"id":54801008,"name":"Biological sciences/Ecology"},{"id":54801009,"name":"Earth and environmental sciences/Ecology"},{"id":54801010,"name":"Earth and environmental sciences/Environmental sciences"}],"tags":[],"updatedAt":"2026-05-11T14:38:39+00:00","versionOfRecord":[],"versionCreatedAt":"2025-09-18 10:40:57","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-7368410","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-7368410","identity":"rs-7368410","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}