Ultra-high resolution pollen record from the northern Andes reveals rapid shifts in montane climates within the last two glacial cycles

Here we developed a composite pollen-based record of altitudinal vegetation changes from Lake F ´ (5 N) in Colombia at 2540 m elevation. We quantitatively calibrated Arboreal Pollen percentages (AP%) into mean annual temperature (MAT) changes with an unprecedented 60-year resolution over the past 284 000 years. An age model for the AP% record was constructed using frequency analysis in the depth domain and tuning of the distinct obliquity-related variations to the latest marine oxygen iso- tope stacked record. The reconstructed MAT record largely concurs with the 100 and 41-kyr (obliquity) paced glacial cycles and is superimposed by extreme changes of up to 7 to 10 Celsius within a few hundred years at the major glacial terminations and during marine isotope stage 3, sug- gesting an unprecedented North Atlantic - equatorial link. Using intermediate complexity transient climate modelling experiments, we demonstrate that ice volume and greenhouse

Sediment sequences from the deep Colombian basins have shown to present a complete Pleistocene archive of montane vegetation changes (Hooghiemstra, 1984;Mommersteeg, 1998;Torres et al., 2005;Van der Hammen, 1974;Van der Hammen and González, 1960;Van Geel and Van der Hammen, 1973;Wille et al., 2001). In particular, sediment records from Lake Fúquene in the Eastern Cordillera of Colombia (5 • 28 N, 73 • 45 W, 2540 m) have revealed detailed altitudinal vegetation distribution patterns since the last glacial maximum (LGM) ( Van der Hammen, 1974;Boom et al., 2002;Hooghiemstra and Van der Hammen, 2004). For a wider perspective on South American past climate variability at shorter time scales as the record presented here, we refer to Wright Jr. et al. (1993), Markgraf (2001) and Vimeux et al. (2009). Here we elaborate on the previous paleoecological investigations of Lake Fúquene by establishing an ultra-high resolution (i.e. ∼60 year between 4646 time slices) pollen-based record of climate change to resolve the orbital and sub-Milankovitch mean annual temperature variations in the northern Andes over the past 284 000 years. In addition, using the separate forcings of orbitally induced insolation, greenhouse gas and Northern Hemisphere ice-sheet variations (e.g., Yoshimori et al., 2001;Timmermann et al., 2009), we have carried out three transient climate modelling experiments to unravel the cause(s) of the reconstructed MAT changes during the last two glacial cycles.

Climate setting of Lake Fúquene
The modern precipitation regime at Lake Fúquene ( Fig. 1) is controlled by the annual migration of the inter tropical convergence zone (ITCZ) causing two dry seasons (December to February and from June to August) and two rainy seasons (March to May and from September to November). The seasonal temperature cycle is very weak with mean monthly temperatures of ∼13.5 ± 0.5 • C. The daily temperature range is large and during the dry season night frost may occur (Van Geel and Van der Hammen, 1973). At present the lake lies within the Andean forest belt (Fig. 2). The upper boundary of this belt or upper forest line (UFL) coincides approximately with the 9.5 • C mean annual isotherm, while the lower boundary is at an elevation where night frost no longer occurs (Hooghiemstra, 1984;Van der Hammen, 1974;Van der Hammen and González, 1960). During glacial conditions, lower temperatures cause a descent in the altitudinal position of individual plant taxa (Hooghiemstra, 1984;Van't Fig. 2. Distribution pattern of source areas of pollen taxa during full glacial conditions and modern (interglacial) times. After Van der Hammen (1974) with additions after Wille et al. (2001). Veer and Hooghiemstra, 2000;Van't Veer et al., 1995;Van der Hammen, 1974;Van der Hammen and González, 1960;Wille et al., 2001).

Sediment cores
Ten meters apart, two ∼60 m long sediment cores, Fúquene-9 (Fq-9) and Fúquene-10 (Fq-10) were retrieved from Lake Fúquene, using a floating platform with Longyear drilling equipment from Gavesa Drilling Co. Bogotá. Consolidated sediments were first approached at c. 6 m below the water surface. Sediments were retrieved in segments of 100 cm length with a diameter of 75 mm. Core samples at the two drilling sites were collected with 50 cm overlapping depth intervals to maximize sediment recovery (see Table 1 in the Supplement). Undisturbed sediments in pvc-tubes were directly transported by air freight to The Netherlands for further treatment. The fresh sediment surface was photographed in a standardized photographic room. The two cores were transported to the NIOZ laboratory (Texel, The Netherlands) to obtain along the full length of both cores XRF-based geochemical data. Subsequently, the cores were transported to the University of Amsterdam for collecting >5000 samples at 1 cm increments for pollen and grain size analysis. Grain size analysis was carried out at the Vrije Universiteit Amsterdam.

XRF analysis
Bulk chemistry was measured with an Avaatech X-ray fluorescence (XRF) core scanner at the Royal Netherlands Institute for Sea Research (NIOZ). The XRF core cortex-scanner counts the number of the chemical elements aluminum (Al, atomic number 13) to bismuth (Bi, atomic number 83) per second (cps) directly at the surface of a split sediment core, a measurement which is proportional to chemical concentrations (Jansen et al., 1998). Prior to the measurement, the split core surface was smoothed horizontally without contaminating sediment surface. Subsequently, the surface was covered with a 4 µm thin SPEXCerti Prep Ultralene foil to avoid contamination of the X-ray unit during measurement and to avoid desiccation of the sediment. Air bubbles under the foil were carefully removed. Measurements were carried out at 1 cm increments along the full length of the 60 m long cores. We used generator settings of 10 kV and 30 kV and measurement time was 30 s per cm. The standard procedure included a control measurement with a standard after every 1 m core interval. Further technical and practical details about the XRF core scanner were described in Richter et al. (2006).

Composite section (Fq-9C)
Cores Fq-9 and Fq-10 were used to build a composite record with a minimal number of gaps in the sedimentary sequence. Down core changes in the lithology represent the first information for this exercise. Distinct and repetitive layers with peat and intervals with higher concentrations of aeolian dispersed fine grained volcanic ash allowed an adequate first correlation between cores. Subsequently, we used records of Fe and Zr obtained by X-ray fluorescence at 1 cm distance over the full length of both cores to fine tune the correlation. Selection of iron (Fe) and zircon (Zn) out of the suite of measured elements is justified by their physical and chemical properties. During XRF measurements heavier elements (Fe and Zr) remained relatively unaffected by the variation of physical properties along the core. In addition, Fe and Zr content may be indicative of variations in source areas and/or variations in sedimentary environments. It is to be expected that changes in both variables  (Sarmiento et al., 2008). The Fe content may also be influenced by changes in the redox state of the water column as well as the sediment columns (Davidson, 1993). The latter may be associated with lake level changes which cause alternations between submersed lacustrine environments to shallow swampy conditions, and even to a drained status of the lake.
Zircon is a conservative element and relatively resistant to chemical weathering processes (Balan et al., 2001). Zircon is found as detrital grains in igneous, metamorphic, and sedimentary rocks. The zircon content is positively correlated with weight percent of coarse silt plus sand (Alfonso et al., 2006;Nyakairu and Koeberl, 2002;Stiles et al., 2003). Therefore, Zr may reflect high energetic sedimentary environments mostly coinciding with a proximal sediment source in relation to the drilling location of the cores. Core Fq-9 had the least technical drilling artifacts and was therefore used as the backbone for our study. This implies that the depth of Fq-10 was adjusted so that the patterns of the various proxy records from both cores aligned. The procedure was carried out as follows: (1) High resolution photographs of the freshly split sediment cores and binocularbased lithological descriptions (Sarmiento et al., 2008) were used to obtain an initial framework of stratigraphic correlation.
(2) Time series of Fe and Zr content of Fq-9 and Fq-10 were visually matched using Analyseries 1.2 software (Paillard et al., 1996). Tie points were preferably chosen at the steepest parts of the curves; occasionally maxima or minima were used (see Table 2 in the Supplement). The matched records were compared with the initial stratigraphic framework.
(3) The procedure described in points (1) and (2) was repeated until all paired proxy records, e.g. Fe, showed comparable variation at any given depth. (4) No adjustments, such as squeezing and stretching were introduced within core segments. All core segments of Fq-10 were depth-shifted and stratigraphically aligned relative to the Fq-9 core segments. Subsequently, a final Fúquene-9 Composite record (labeled Fq-9C) was built where inclusion of disturbed intervals and sediment gaps was minimized (see Table 3 in the Supplement).
The Fq-9C pollen record represents approximately 90% of the sediment infill of the uppermost 60 m of the Fúquene Basin; the remaining 10% reflect small coring gaps, inadequate sediment intervals, and not analyzed parts due to organizational problems. The average lateral offset between stratigraphic layers was 52 cm, expressing differences in sedimentation rates, erosion, and methane emissions during the drilling procedure (Fig. 3). Sediments between 26 and 21 m contain significant proportions of organic matter. Interval 22 to 20.8 m reflects compressed peat. This 5 meter of sediments contained over-pressured methane which escaped from the borehole when the corer contacted this reservoir at 21 m depth. For safety reasons, drilling activities were stopped for a full day until sediments dissipated (GAVESA, 2002). The escape of methane and subsequent compaction of the peat explains the significant shift in the offset between both cores at 24-26 m (Fig. 3). The lake surface was used as a reference. During the 5 week long drilling operation, lakelevel changes were noticed after periods of heavy rains. As a consequence, the offset between the zero calibration points of cores Fq-9 and Fq-10 is estimated offset up to 20 cm.

Pollen analysis
Pollen samples of 1 cm 3 were processed using the standard pre-treatment including sodium pyrophosphate, acetolysis, and heavy liquid (bromoform) separation. We counted pollen and spore taxa with known ecological envelopes and a clear response to climate change through altitudinal shifts (Hooghiemstra, 1984;Mommersteeg, 1998;Torres et al., 2005;Van't Veer and Hooghiemstra, 2000; Van der Hammen   Table 2 in the Supplement. Gaps between vertical red lines represent mineral-rich sediments. and González, 1960;Wille et al., 2001). Pollen types were assigned to the following ecological groups: (1) taxa of subandean forest, (2) taxa of Andean forest, (3) taxa of subpáramo vegetation, (4) taxa of grasspáramo vegetation and (5) taxa indicating dry conditions. Down core changes in the relative contribution of the pollen types in www.clim-past.net/7/299/2011/ these ecological groups reflect shifts in their altitudinal distribution. Following Van der Hammen and González (1960) and Hooghiemstra (1984), the measured AP% were used to reconstruct the position of UFL along the record. The AP% of Fq-9C fluctuates between 2% and 98% with the highest values found between 21 and 26 mcd and the lowest around 10, 29 and 53 mcd (Fig. 4a).

Spectral analysis
Spectral analyses was carried out on the Fq-9C AP% record in the depth domain using the CLEAN algorithm (Roberts et al., 1987), the REDFIT program (Schulz and Mudelsee, 2002) and a Blackman-Tukey power spectrum (Blackman and Tuckey, 1958) to identify the possible imprint of orbital signals. The CLEAN algorithm was run as a MAT-LAB routine (Heslop and Dekkers, 2002) and is particular robust in determine the frequency distribution of unevenly spaced data series and noisy signals (Baisch and Bokelmann, 1999;Roberts et al., 1987). The CLEAN procedure performs a nonlinear deconvolution in the frequency domain in order to remove any artifacts resulting from incomplete sampling in the time (depth) domain. It includes a series of Monte Carlo simulations which allow a large number of slightly different spectra to be generated for a single input signal. The differences between these spectra are utilized to determine a mean spectrum and confidence intervals both for individual frequency peaks and for the spectrum as a whole. Through the use of an Inverse Fourier Transform of the MC-CLEAN spectrum, the data can be reconstructed in the time domain, providing a "noise free" version of the input signal. Because of the large number of data points and the almost equally spaced sampling resolution, we choose to linearly detrend the AP% depth series and slightly smooth it with a 5 cm average to reduce the number of data points from 4586 to 1134 before the CLEAN analysis (Fig. 4b). The CLEAN spectrum was subsequently determined by adding 10% red noise (i.e. control parameter = 0.1), a clean/gain factor of 0.1, 500 CLEAN iterations, an interpolation step (dt) value of 1 cm and 500 simulation iterations (Fig. 5a).
We evaluated the CLEAN spectral outcome first by comparison with a Blackman-Tukey power spectrum, which was performed with AnalySeries 2.2 (Paillard et al., 1996) using 75% of the 5-cm interpolated and detrended time series, a Parzen window and a Band Width of 0.0087 (Fig. 5b). Secondly, we applied the computational spectral analysis program, REDFIT version 3.8 (Schulz and Mudelsee, 2002) on the original (non-smoothed and non-detrended) AP% data set. In this model a first-order autoregressive process (AR1), which is assumed to present a good estimate of the red-noise signature, is estimated directly from the (unevenly spaced) time series and, subsequently, transferred into the frequency domain using the Lomb-Scargel algorithm. Confidence levels based on a χ 2 distribution are calculated from the AR1noise and from percentiles of the Monte Carlo ensemble. For the analysis, we applied one Welch Overlapped Segment Averaging (WOSA) segment (n 50 = 1), 2000 Monte Carlo simulations (N sim = 2000), and a rectangular window (I win = 0). This resulted in a value of τ of 20.9 with 2 • of freedom (Fig. 5d). Third, we ran the CLEAN algorithm again, but in this case, Dave Heslop (personal communication, 2009) incorporated the method of Mudelsee to estimate the AR1 characteristic period of the input data series . With each iteration step, a new time series is calculated with the same AR1 characteristics and processed with CLEAN. These spectra are then used to calculate the confidence levels. This approach is a little different from before, because the noise is not added to the signal, but studied separately. This implies that the spectrum for the input data does not have error bounds, because it remains the same for each iteration step. Only the separate noise component is changing.
The results are plotted in Fig. 5c and largely confirm the significance levels of the spectral peaks obtained from REDFIT.

Climate model and experimental set-up
We performed three transient climate modelling experiments over the past 284 ka to serve as an explanation for the temporal variations in the AP% record. For this purpose we used a coupled model of intermediate complexity, CLIMBER-2 (version 3) (Petoukhov et al., 2000), that is very suitable for long transient simulations due to its fast turnaround time (Calov et al., 2005a,b;Claussen et al., 1999;Tuenter et al., 2005). The model consists of an atmosphere model, an ocean/sea-ice model, and a land/vegetation model. No flux adjustments are used. The atmospheric model is a 2.5dimensional statistical-dynamical model with a resolution of 10 • in latitude and approximately 51 • in longitude. The vertical resolution for the circulation, temperature and humidity is 10 levels and for the radiation, 16 levels. The time step is one day. The terrestrial vegetation model is VECODE (VEgetation COntinuous Description) (Brovkin et al., 1997). The model computes the fraction of the potential vegetation (i.e., grass, trees and bare soil) with a time-step of one year. The ocean model (Stocker et al., 1992) contains a thermodynamic sea-ice model and describes the zonally averaged temperature, salinity and velocity for three separate basins (Atlantic, Indian and Pacific oceans). The latitudinal resolution is 2.5 • and the vertical resolution is 20 unequal levels. Results of CLIMBER-2 compare favorably with data of the present day climate (Ganopolski et al., 1998;Petoukhov et al., 2000), as do the results derived from sensitivity experiments (like changes in vegetation cover and solar irradiance) with those of more comprehensive models (Ganopolski and Rahmstorf, 2001). At low latitudes the monsoonal strengths in CLIMBER-2 during the mid-Holocene lies within the range of GCMs (Braconnot et al., 2002) and it also simulates realistic transient behaviour of the monsoons (Claussen et al., 1999). The model thus adequately simulates many aspects of large-scale climate variability.  With CLIMBER-2 we performed three transient simulations for the interval from 650 ka to present, of which we show the last 284 kyr. The only forcing used in the EXP O simulation is insolation changes induced by the La04 orbital parameters (Laskar et al., 2004), while the ice sheets volume and greenhouse gas forcing were kept fixed at presentday and pre-industrial values, respectively. In the EXP OI simulation the same orbital forcing was used with varying ice-sheets on the Northern Hemisphere. In the EXP OIG simulation the same orbital forcing was used with varying ice-sheets and changes in the atmospheric greenhouse gas concentrations of carbon dioxide (CO 2 ) and methane (CH 4 ) derived from ice core data.
The greenhouse gas concentrations and ice-sheet variations were prescribed because CLIMBER-2 does not contain a carbon cycle model and interactive ice-sheet model. Prior to 500 yr BP the concentrations of CO 2 and CH 4 used were obtained from the Antarctica ice core EPICA Dome C with the EDC 3 timescale (Lüthi et al., 2008;Loulergue et al., 2008). For the last 500 years we used values from several sources (Robertson et al., 2001). For EXP OI and OIG we prescribed the ice area and height of the Eurasian and North American (Laurentide) ice sheet, which are considered to portray the largest ice sheet fluctuations during the studied time interval. Hence, for all simulations the height and surface area of Greenland and Antarctica as well as small glaciers were kept at present-day values. Ice-sheet volumes of the Laurentide and Eurasian ice sheet were obtained from a simulation with a 3-dimensional ice sheet model coupled to a deep-ocean temperature model (Bintanja et al., 2005). This model was forced by the LR04 benthic δ 18 O record (Lisiecki and Raymo, 2005) in an inversed mode resulting in time series of the NH ice-sheet volumes and air temperatures (Fig. 6). From the simulated ice volumes the time-varying areas and heights of the American and the Eurasian ice-sheet had to be reconstructed. For this purpose, the ice-sheet areas defined by the ICE-5G ice distribution from the LGM until present (Peltier, 2004) were set on the spatial scale of CLIMBER-2 as a reference, and extrapolated back to 650 kyr BP pending the ice sheet volumes. In this way, the total ice-sheet area as well as the area for each ice sheet included the same temporal behavior as the volumes reconstructed by the 3-D ice sheet model. From the areas and volumes of the ice-sheets, their heights can be reconstructed.
To identify the separate influence of variations in albedo and height due to the variations in ice-sheets, two additional experiments were carried out in which the heights of the ice sheets were changed with the fixed albedo of vegetation and soil, and where the albedo was changed with fixed presentday heights of the ice sheets. The first experiment shows that variations in height mainly affect the atmospheric circulation, causing climate variations especially above and to the east (i.e., downstream) of the ice-sheets. In contrast, the albedo-related changes influence climate both at high and low latitudes, including northern South America.
During the (prescribed) waxing and waning of the ice sheets there is no transport of water from the oceans to the ice sheets and vice versa, i.e., the sea-level in the model does not change during glacial cycles. The results are shown as averages over 100 yr as the periods of the oscillations of the orbital forcing and variations in ice-sheet volume and greenhouse gas concentration are much longer than 100 yr. See Weber and Tuenter (2011) for a detailed description of the boundary conditions and simulations.

Orbital tuning of the Fúquene Basin Composite (FqBC) record
The CLEAN spectrum of the Fq-9C AP(%) record reveals highly significant (99%) peaks at 9.07 and 22.65 m, and peaks with lower significance at 12.58, 5.96, 4.19 and 3.65 m (Fig. 5a). Although the power distribution is much more smoothed in the BT spectrum, strong peaks occur at ∼9 and 25 m; consistent with the CLEAN estimates (Fig. 5b). Also the REDFIT power estimate and distribution are in close agreement with that of the BT and CLEAN methods, although an additional peak occurs at 1193 cm ( Fig. 5b and d). The 22.65 m periodicity coincides with the largescale changes in the AP% record, which we attribute to the imprint of the late Pleistocene ∼100 kyr glacial rhythm . Accordingly, the 9.07 m cycle corresponds with a 41 kyr period, indicating a large obliquity control of the climate variability in this region. We therefore choose to use only the obliquity-related changes in the AP% record to construct an astronomical tuned time scale for Fq-9C. As a starting point for our tuning, we extracted the 9 m component of the AP% record using a Gaussian filter, as implemented in the freeware AnalySeries 2.2 (Paillard et al., 1996) centered at a frequency of 0.0011 ± 0.0004 (Fig. 3b). The filtered 9 m signal was subsequently correlated to the filtered obliquity-related 41 kyr component of the LR04 benthic δ 18 O record (Lisiecki and Raymo, 2005), to establish an age model for Fq-9C (Fig. 3c). For this purpose we applied a Gaussian filter centered at a frequency of 0.0245 ± 0.002, and tuned our data by peak-to-low matching of the filtered 9-m signal of the Fq-9C AP% depth series and the filtered 41-kyr signal of the LR04 time series (Table 1). Extrapolation of the resulting age model provides an age estimate for the bottom and top of the Fq-9C sedimentary sequence of respectively 284 and 27 kyr before present (ka), and an average sample resolution of ∼60 yr.
Data from the last 27 kyr of the 12 m long Fúquene-2 (Fq-2) core (Van Geel and Van der Hammen, 1973) was implemented on top of the Fq-9C to extend the AP% record up to the latest Holocene. The resulting Fúquene Basin Composite (FqBC) reveals a complete and ultra high-resolution AP% record for the past 284 000 years. These correlations were verified with the 18 m long Fúquene-7C (Fq-7C) core obtained from the lake shore at approximately 1.5 km distance from Fq-9 and Fq-10 (Hessler et al., 2010;Mommersteeg, 1998). For this purpose, we revised the 14 C ages for the Fq-2 and Fq-7C records (Table 1), using the CALIB REV 5.0.2. program (Stuiver et al., 2009). Tie points between the Fq-2, Fq-7C and Fq-9C AP% records were based on the following biostratigraphic markers: Arboreal Pollen, Alnus, Polylepis-Acaena, Quercus, Myrica, Podocarpus, Asteraceae tubuliflorae, Hypericum (Table 1). The comparison between the AP% records of Fq-2 and Fq-7C showed that the Younger Dryas is only present in core Fq-2 (not shown). Therefore, we linked Fq-9C directly to Fq-2 at 27.19 ka to produce the FqBC record (Fig. 6).
Correlation to the LR04 δ 18 O benthic record was chosen, because this record is used as the backbone for many late Pleistocene paleoclimate studies. The LR04 chronology follows the SPECMAP approach Imbrie and Imbrie, 1980), in which the δ 18 O benthic record is tuned to a simple ice sheet model that includes a forcing function (i.e. 21 June insolation curve for 65 • N), an average ice-sheet response time and a non-linearity coefficient that describes the slow build-up and fast terminations of the ice-sheets. For the past 300 000 years, the LR04 chronology yields time lags for the obliquity (41-kyr) and precession (23-kyr) components of the δ 18 O benthic record of 7.5 ± 0.8 and 4.5 ± 0.5 kyr, respectively. This implies that also the AP% time series includes a ∼7.5 kyr time lag for its obliquity-related component, which is supported by the compliance between the extrapolated tuned ages of positions 300 and 353 cm in Fq-9C of respectively 32.9 and 35.2 ka, and their corresponding corrected radiocarbon ages of 32.54 ± 0.32 and 35.74 ± 0.31 ka (Fig. 6).
The AP% time series improved the temporal resolution of the earlier pollen-based records from the Fúquene Basin with better than centennial resolution. Evidently, the FqBC record depicts the last three glacial terminations, T I -T I I , T IIIa and T IIIb (Fig. 5). As expected on the basis of the spectral results in the depth domain, highly significant spectral power is found at the 41 kyr and 113 kyr periods, while a clear imprint of the precession cycle is lacking (Fig. 5e). In addition, wavelet analysis shows a distinct 8-kyr period appearing only at the major terminations, which is actually part of a large range of periods, reflecting the rapid and extreme changes in the AP% record at these times (Fig. 7a).

Mean annual temperature reconstructions
From previous altitudinal pollen studies of the tropical Andes, it appeared that changes in AP% respond quasilinearly to temperature-driven vertical shifts in the UFL between 3500 m (the highest mountains at close distance) and the LGM position at ∼2000 m (Hooghiemstra, 1984;Van't Veer and Hooghiemstra, 2000;Van der Hammen, 1974;Van der Hammen and González, 1960;Wille et al., 2001). Accordingly, MAT at Lake Fúquene was approximately 7.8 • C lower during the LGM than at present, which is in close  (Thompson and Goldstein, 2006). agreement with the 5-9 • C decrease derived from the change in snowline along the Eastern Cordilleras of the central Andes (Broecker and Denton, 1990;Klein et al., 1999), and the ice core δ 18 O data of the Huascarán (Peru) and Sajama (Bolivia) glaciers in the southern Andes (Thompson et al., 1995(Thompson et al., , 1998(Thompson et al., , 2000Broecker, 1995). In addition, Wille et al. (2001) showed that the LGM lapse rate was much steeper (∼0.76 • C/100 m) than today (0.6-0.64 • C/100 m; Florez, 1986). This change in lapse rate would imply that the air temperatures at sea level decreased by 3 to 5 • C during the LGM, which is consistent with the ∼3 ± 1 • C sea surface temperature (SST) lowering derived from the Mg/Ca and U k 37 -based temperature reconstructions of cores TR163-19 (Lea et al., 2000) and MD02-2529 (Leduc et al., 2007) in the Equatorial Pacific (Figs. 1 and 6e).
The AP% record may thus be used to provide a good approximation of MAT at Lake Fúquene through time, using the modern to glacial temperature difference, considering that the seasonal variations in mean monthly temperatures remained as small as today. However, over the last 3000 years, the AP% record has been biased by the destruction of montane forests through intensive agriculture and soil erosion (Van Geel and Van der Hammen, 1973). Therefore, we assigned the AP% of 73 ± 6% at ∼3.1 ka a MAT of 13.5 ± 0.5 • C, assuming that MAT remained close to presentday values over the past 3 ka (Van Geel and Van der Hammen, 1973). With a mean AP% value of 15 ± 6% at ∼20 ka, this implies that a change in AP% of 10% corresponds with a MAT change of 1.3 ± 0.3 • C. The estimated standard error of the presented MAT record (Fig. 6) is 0.6 ± 0.4 • C, by taking into account a mean temperature difference between 20 and 3 ka of 7.8 • C and the estimated errors in the lapse rates and AP% (Hooghiemstra, 1984;Van't Veer and Hooghiemstra, 2000;Van der Hammen, 1974;Van der Hammen and González, 1960;Wille et al., 2001). The resulting MAT estimate of 15 ± 1.5 • C at 7 ka (early Holocene) is consistent with earlier reconstructions (Van Geel and Van der Hammen, 1973).
Comparison between our reconstructed MAT at Lake Fúquene and the Mg/Ca-derived SST estimates of core TR163-19 for the last 284 000 years shows that temperature variations at high altitudes in the northern Andes are larger and much more rapid, i.e. up to 10 ± 2 • C within a few hundred of years at T I I and T IIIb , than reflected in the equatorial marine record (Fig. 6e). This extreme temperature change is exceptional for the low latitudes, and has yet only been found in high latitude temperature records of Antarctica (Jouzel et al., 2007;Parrenin et al., 2007) or the estimated continental mean annual air temperature changes between 40 and 80 • N (Bintanja et al., 2005) at these terminations ( Fig. 6b and d).
The temperature changes of up to 7 • C during MIS 3 are also exceptional and are yet only observed in Greenland ice core records (Fig. 8). We further obtained a considerable longer duration for MIS 5.5 (defined by the temperatures above present-day values between 110-133 ka) than the marine and polar temperature records (120-132 ka). A longer duration for MIS 5.5 is, however, in close agreement with the radiometric dated sea level records (Blanchon et al., 2009;Gallup et al., 2002;Thompson and Goldstein, 2006) during this period (Fig. 7b).

Transient climate modelling experiments
EXP O reveals only small temperature differences of less than ∼0.8 • C, which oscillate dominantly on a precession frequency (Figs. 5f and 7c). The influence of orbital-induced insolation changes alone can therefore not explain the reconstructed MAT shifts at Lake Fúquene, since the FqBC AP% record lacks a (dominant) precession-related signal (Fig. 5e). The spectral fit between the reconstructed and the modelled MAT becomes much better in case ice volume changes are taken into account (EXP OI). Although the 41 and 100 kyr dominated ice volume changes largely describe the pattern in the MAT at ∼2.5 km altitude in the tropical Andes, alone they are insufficient to explain the whole magnitude of the reconstructed changes (Figs. 5f and 7c). Evidently, the modelled MAT compares much better with the orbital variations in the reconstructed MAT record when greenhouse gas forcing (EXP OIG) is added (Fig. 7c). In particular, the relative spectral power of the 41 and 113 kyr increased significantly with respect to the precession (19-23 kyr) related peak in EXP OIG (Fig. 5f).

Discussion
On a global scale, our new FqBC record improved the temporal resolution of long pollen-based records of climate change with an order of magnitude. The record shows a hitherto unknown level of variability in the altitudinal distribution of forest and páramo taxa and plant associations at both orbital and sub-Milankovitch time scales.

Orbital-driven MAT changes
The high sensitivity of MAT at high altitudes in the northern Andes to ice volume and greenhouse gas variations on orbital time scales is supported by recent modelling studies (Urrutia and Vuille, 2009), which project large changes in South American (sub)alpine climates by the end of the 21st century due to enhanced anthropogenic greenhouse gas emissions. Still, the simulated glacial-interglacial MAT changes of 3 to 4 • C significantly underestimate the reconstructed variations at Lake Fúquene ( Fig. 7b and c). Part of this discrepancy can be explained by the large divergence between simulated glacial-interglacial changes at a lapse rate of less than 0.005 • C/100 m and the reconstructed change at a lapse rate of up to ∼0.3 • C/100 m. Another important factor in controlling this offset is the low spatial resolution of CLIMBER, which excludes to resolve specific changes in the local hydrology and vegetation feedbacks within the studied region. Moreover, due to its statistical-dynamical approach, CLIMBER-2 runs strongly underestimate the sub-Milankovitch and millennial scale variability (i.e. <11 kyr), which clearly affected the MAT at the lake. Remarkable is the absence of a distinct precession-related signal in the FqBC AP% record over the past 284 000 years, since this astronomical cycle dominates tropical climate variations Short et al., 1991). This accounts in particular for tropical regions, where the seasonal migration of the ITCZ and the strength of the monsoons vary significantly (Kutzbach and Street-Perrot, 1985;Tuenter et al., 2005). Although possible changes in sedimentation rates at precession time scales or the few gaps in the record may have weakened the imprint of the precession signal, the close correspondence between our reconstructed and modelled MAT changes suggests that the obliquity-driven and ∼100 kyr dominated global changes in ice volume and greenhouse gas forcing were much more significant for the ecological changes at 2540 m elevation in the northern Andes than the precession-driven changes in local/regional temperature and precipitation. This is in contrast to records from the Bolivian Altiplano and northeastern Brazil, which indicate significant precession-forced precipitation changes (Baker et al., 2001;Wang et al., 2008). Possibly, the latitudinal displacement of the ITCZ during summer and winter was small in the northern Andes on precession time scales, thereby causing only minor precipitation changes that could have affected the vegetation patterns at the surroundings of Lake Fúquene.

Sub-Milankovitch related MAT changes
Ice-cores provide some of the highest temporally resolved records of late Pleistocene climatic variability. Oxygen isotope records from Greenland ice-cores record rapid transitions between cold stadial and warm interstadial conditions in the Northern Hemisphere during the last glacial and deglacial period (GRIP-Members, 1993;Grootes et al., 1993). These so-called Dansgaard-Oeschger (DO) events represent episodes of extremely rapid warming (up to 10 • C temperature rise within a few decades), followed by more gradual cooling. DO events have been linked to variations in the mode of North Atlantic deep-water formation, driven by changes in the supply of freshwater to the region of deepwater formation (Broecker et al., 1985), which may have culminated into short periods of massive iceberg discharge, named Heinrich (H) events (Fig. 8), in the North Atlantic (Broecker, 1994). Corresponding climatic shifts have been found in low-latitude records, e.g. showing southward shifts in the annual average ITCZ (Peterson et al., 2000) and reduced Asian monsoon rainfall (Wang et al., 2008) during stadials. In contrast to the changes observed in Greenland, glacial age temperature variability recorded by Antarctic ice cores is characterized by a more gradual and symmetric behaviour that is out-of-phase with the changes in the Northern Hemisphere (Blunier and Brook, 2001). This relationship provides the basis for the so-called bipolar seesaw hypothesis (Stocker, 1998), whereby changes in ocean circulation associated with cooling across the North Atlantic drive a corresponding warming across the Southern Ocean and Antarctica.
The clear signature of the Younger Dryas, constrained by 14 C dates, and the interstadial DO cycles 1 (Bølling-Allerød), 8,12,14,19 and 20 in the reconstructed MAT record of Lake Fúquene suggests an unprecedented North Atlanticequatorial link (Fig. 8). In addition, the short interval with low MAT during MIS 5.5 and the rapid MAT changes during the penultimate glacial period and Termination II mirrors the Greenland Ice Core Project (GRIP) δ 18 O record (GRIP-Members, 1993). The lower part of the GRIP core is, however, suspect to disturbance, since it shows a different pattern than the one found in the North Greenland Ice Core Project (NGRIP; Svensson et al., 2008;NGRIP Members, 2004) and the nearby Greenland Ice Sheet Program 2 (GISP2; Grootes et al., 1993), although gas measurements suggest that it contains ice from the last interglacial and penultimate glacial maximum (Landais et al., 2003). The North Greenland At present, we consider that the age constraints of our MAT record are not accurate enough to determine the exact phase relationship with the North Atlantic cold events: i.e. the warm events appear also closely linked to the inferred Antarctic Isotope Maximum (AIM; EPICA Dome C Members, 2006) (Fig. 8). It is tempting, however, to link maximum MAT conditions at Lake Fúquene to interstadial periods in line with the δ 18 O ice records of the Huascarán (Peru) and Sajama (Bolivia) glaciers in the southern Andes (Thompson et al., 1995(Thompson et al., , 1998(Thompson et al., , 2000. In addition, palynological investigations of the Cariaco Basin off northern Venezuela (Fig. 1) have revealed the highest pollen concentrations and the maximum extent of semi-deciduous and evergreen forests in the northernmost part of South America during these times (González et al., 2008). During stadials, the region around the Cariaco Basin is characterized by increases of salt marshes, herbs, and montane forests, while during Heinrich events forest abundance decreased (González et al., 2008). It has been proposed that during these events, a reduced Atlantic meridional overturning circulation resulted in extreme winter cooling of the North Atlantic (Cheng et al., 2006;Denton et al., 2005). Through an atmospheric connection, the ITCZ was, probably also with a winter bias (Ziegler et al., 2008), and shifted to a more southern position (Cane and Clement, 1999;Chiang et al., 2003;Clement and Peterson, 2008;Peterson et al., 2000;Hessler et al., 2010), and causing wetter climate conditions in the northeastern part of Brazil and the Bolivian Altiplano (Baker et al., 2001;Wang et al., 2004). A comparison with a detailed record of North Atlantic, C 37:4 alkenone record of the Iberian Margin (Martrat et al., 2007), shows that in particular during H1-2 and H6 Lake Fúquene was affected by the lowest MAT (Fig. 8).

Conclusions
A new ultra-high resolution pollen record from the Fúquene Basin in the northern Andes over the past 284 000 years has been presented and found to reveal a distinct one-to-one coupling between tropical MAT and the North Atlantic climate variability on both orbital and millennial time scales. Spectral analysis and orbital tuning techniques showed that the long-term MAT changes at the lake concur with the 41 kyr and ∼100 kyr glacial-bound ice volume changes. Transient experiments with a climate model of intermediate complexity (CLIMBER-2) revealed in addition that the large-scale, orbital-induced vegetation changes cannot be explained by ice volume variations alone, but that the modelled and reconstructed MAT fits better when greenhouse gas forcing is added. However, the simulated glacial-interglacial MAT changes are less amplified with respect to the reconstructed variations, probably because the model underestimates the glacial-interglacial changes in lapse rate. In addition, the low spatial resolution of CLIMBER-2 excludes resolving specific changes in the local hydrology and vegetation feedbacks.
Superimposed on the orbital-scale variations, the reconstructed MAT record revealed rapid (within a century) and extreme temperature changes of up to 10 ± 2 • C at the major glacial terminations and during MIS 3. Evidently, tropical montane forest and tropical alpine grassland (páramo) ecosystems are very sensitive and react with significant altitudinal migrations on global climate change on orbital to millennial and centennial time scales, thereby supporting evidence from tropical ice core records and current projections of future climate change. Hence, lake sediments have an enormous potential to reach the accuracy of marine and ice core records when traditional bottlenecks in analysing long continental records are adequately addressed.