A global analysis of pollen-based reconstructions of land climate changes during Dansgaard–Oeschger events

Dansgaard–Oeschger (D–O) warming events are comparable in magnitude and rate to the anticipated 21st century warming. As such, they provide a good target for evaluation of the ability of state-of-the-art climate models to simulate rapid climate changes. Despite the wealth of qualitative information about climate changes during the D–O events, there has been no attempt to date to make quantitative reconstructions globally. Here we use frequency-corrected Tolerance-weighted Weighted Averaging Partial Least Squares regression (fxTWA-PLS) to reconstruct mean temperature of the coldest month, mean temperature of the warmest month, and a plant-available moisture index across multiple D–O events between 50 and 30 ka based on available pollen records across the globe. The reconstruction of plant-available moisture is corrected for the impact of changing atmospheric CO₂ concentrations on plant water use efficiency. These reconstructions show that the largest warming occurred in northern extratropics, especially Eurasia, while western North America and the southern extratropics were characterised by cooling. The change in winter temperature was significantly larger than the change in summer temperature in the northern extratropics and the tropics, indicating that the D–O warming events were characterised by reduced seasonality, but there was no significant difference between the summer and winter temperature changes in the southern extratropics. The antiphasing between northern and southern extratropical changes, and the west-east pattern of cooling and warming in North America were generally consistent across the eight D–O events examined, although coherency is greatest during the strongest events. There was no globally consistent pattern between changes in moisture and changes in temperature. These reconstructions can be used to evaluate the spatial patterns of changes in temperature and moisture in the transient simulations of the D–O events planned as part of the Palaeoclimate Modelling Intercomparison Project.

1 Introduction

Dansgaard–Oeschger (D–O) events are characterised in Greenland by a transition from cold Greenland Stadial (GS) to warmer Greenland Interstadial (GI) conditions (Dansgaard et al., 1993). The surface air temperature in Greenland increased by 10–16 °C during the warming phases; these warming events occurred over an interval of between 50 and 200 years (Huber et al., 2006; Kindler et al., 2014). Thus, the D–O events offer a parallel in terms of speed to projected future warming, although both the baseline state and the mechanism inducing this warming differ from anticipated 21st century climate changes. D–O events could therefore provide an opportunity to determine how well climate models that are used for future projections can simulate rapid climate changes (Malmierca-Vallet et al., 2023), particularly regional patterns of warming (and cooling) that are regarded as a challenge for modelling (Doblas-Reyes et al., 2021; Lee et al., 2021) and are highly important in assessing the vulnerability of human societies to future climate changes (IPCC, 2022).

Although D–O events are found throughout the last glacial period, the largest number and the most regular patterning occurred during Marine Isotope Stage 3 (MIS 3; 57 to 29 ka) when there were 11 separate events (D–O 15 to D–O 5), while MIS 4 (71 to 57 ka) only had 3 separate events (D–O 18 to 16) (Kindler et al., 2014). The typical duration of a cycle as manifested in Greenland is ca. 1500 years and is characterised by an initial short slow warming, followed by an abrupt large warming in matter of decades, followed by a long slow cooling over centuries to millennia, with a terminal phase of fast cooling (e.g. D–O 8, D–O 12); however, there are also cycles in which the warming and cooling phases took roughly the same time (e.g. D–O 5, D–O 6, D–O 9) (Kindler et al., 2014). The magnitude of changes also differ, with both strong events (e.g. D–O 8, D–O 12) and weak events (e.g. D–O 9) (Kindler et al., 2014).

The D–O signals are not just seen in Greenland – they are registered globally (Adolphi et al., 2018; Corrick et al., 2020; Harrison and Sanchez Goñi, 2010; Sánchez Goñi et al., 2017; Voelker, 2002) and are reflected in changes in both temperature and precipitation. Both oceanic and ice-core records indicate that temperature changes are out-of-phase between the northern and southern hemispheres, and the southern hemisphere response both in terms of warming and cooling phases is generally less abrupt (Dima et al., 2018; Vettoretti and Peltier, 2015). There is a comparative lack of information from the continents about the manifestation of D–O events. Shifts in vegetation types between GI and GS states have been interpreted as primarily a temperature signal in the extratropics and a moisture signal in the tropics (Harrison and Sanchez Goñi, 2010). Speleothem records provide a good time-control on the synchroneity of climate changes globally with the D–O events registered in Greenland (Adolphi et al., 2018; Corrick et al., 2020), but the driver of this signal can either be temperature or precipitation depending on the region. There are quantitative climate reconstructions based on terrestrial pollen records from La Grande Pile (Guiot et al., 1993), Lago Grande di Monticchio (Huntley et al., 1999), Padul (Camuera et al., 2022), El Cañizar de Villarquemado (Camuera et al., 2022; Wei et al., 2021) and Lake Ohrid (Sinopoli et al., 2019), marine cores in the western Mediterranean and offshore from Portugal (Sánchez Goñi et al., 2002), diatom assemblages at Les Echets, France (Ampel et al., 2010), chironomids from Lake Bergsee in central Europe (Lapellegerie et al., 2024), bacterial membrane lipid records from the Eifel region (Zander et al., 2024), isotopic measurements of earthworm calcite from the Rhine Valley (Prud'homme et al., 2022) and clumped isotope measurements on snails in Hungary (Újvári et al., 2021). Aside from the lack of comparable quantitative estimates from outside Europe, differences in the methodology employed and in the specific climate variables reconstructed in each of these studies limits their usefulness for model evaluation. In particular, given that there is still uncertainty as to whether the D–O cycles are characterised by changes in seasonality such that warming events are primarily driven by changes in winter (Flückiger et al., 2008; Zander et al., 2024; Zumaque et al., 2025), in the regional strength of the warming (Harrison and Sanchez Goñi, 2010) and how warming relates to changes in moisture (Wei et al., 2021), there is a need for more systematic reconstruction of seasonal climate changes.

In the paper, we provide reconstructed changes in seasonal temperatures and plant-available moisture during the intervals corresponding to D–O warming events in Greenland between 50 and 30 ka based on available pollen records globally, using a standard regression-based approach, fx-corrected Tolerance-weighted Weighted Averaging Partial Least Squares (fxTWA-PLS: Liu et al., 2020, 2023). We also correct the reconstructions of plant-available moisture to take account of the impact of lower CO₂ on plant water-use efficiency following Prentice et al. (2022a). We analyse the regional patterns to identify key targets for model evaluation.

2 Methods

2.1 Data sources

2.1.1 Modern pollen dataset

Modern pollen data were obtained from version 3 of the SPECIAL Modern Pollen Dataset (SMPDSv3) (Harrison et al., 2025a). This global dataset was constructed by amalgamating and standardising records from public repositories (e.g. Neotoma, Pangaea), existing regional databases (e.g. European Modern Pollen Database, African Pollen Database), individual publications and records provided by the original authors. The records were carefully screened to remove duplicates that were present in more than one source. The modern samples were obtained from multiple types of record, including pollen traps, surface samples, moss polsters and different types of sediment, including cores from lakes and peatbogs, and section through e.g. fluvial or loess deposits. In cases where the record was radiometrically dated, the database preserves all samples younger than 50 yr BP. However, some samples were undated and are therefore recorded as modern if a collection date was given or assumed modern if not.

The dataset contains 26 704 samples from 18 202 different locations, and was created after removing taxa that are not climatically diagnostic (e.g. obligate aquatics, carnivorous species, cultivated plants). The dataset provides several levels of taxonomic aggregation; here we use the most aggregated level, where woody species were generally combined at genus level and herbaceous species at sub-family or family level unless they were palynologically distinctive, occupied distinctive ecological niches and were sufficiently geographically widespread. This “amalgamated” dataset contains relative abundance information for 1367 taxa. These samples were aggregated by location (longitude, latitude and elevation) in order to remove duplicates. Counts for Quercus, Quercus (deciduous) and Quercus (evergreen) were combined because of inconsistent differentiation of Quercus pollen in different regional records. Deciduous and evergreen oaks occupy different areas of climate space, particularly in terms of seasonal moisture; specifically, evergreen oaks are typically found in areas characterised by winter rainfall such as the Mediterranean. Nevertheless, since there are other plant taxa that are similarly diagnostic of such regimes, the amalgamation of Quercus (deciduous) and Quercus (evergreen) should not have a major effect on the robustness of our climate reconstructions. We have tested this assumption by making reconstructions based on all taxa except Quercus (Supplement, Sect. S4). Taxa that occurred in less than 10 samples in the training dataset were not used to make reconstructions because it is unlikely that the available samples provided a reasonable estimate of the climate space occupied by these rare taxa (Liu et al., 2020). After the location aggregation and the taxa filter, the dataset contains information on 18 202 samples with relative abundance information for 609 taxa (Fig. 1a).

Figure 1 Map showing the locations of (a) modern pollen records used to derive the pollen-climate transfer functions, and (b) fossil pollen records covering the interval 50–30 ka used for the climate reconstructions.

We focus on three climate variables: mean temperature of the coldest month (MTCO), mean temperature of the warmest month (MTWA), and a plant-available moisture index (α_plant) defined as the estimated ratio of actual to equilibrium evapotranspiration. These three variables reflect ecophysiological controls on plant distribution (Harrison, 2020; Woodward, 1987) that have been shown to independently influence the distribution and abundance of plant species (Boucher-Lalonde et al., 2012; Wang et al., 2013; Wei et al., 2020). α_plant is a transformation of the commonly used moisture index MI (defined as the estimated ratio of annual precipitation to annual potential evapotranspiration) that emphasizes differences at the dry end of the climate range, which have a more pronounced effect on vegetation distribution than differences at the wet end (Prentice et al., 2017). Thus, α_plant can be better reconstructed from the pollen records than MI.

The climate values at each SMPDSv3 site were obtained using a geographically-weighted regression (GWR) of climatological values of mean monthly temperature, precipitation, and fractional sunshine hours from the Climatic Research Unit Time-Series version 4.04 (CRU TS4.04; Harris et al., 2020) dataset averaged over the period 1961–1990, which corresponds to the interval from which most of the pollen samples were derived. GWR was to correct for elevation differences between the CRU grid cells and the pollen sites. MTCO and MTWA were taken directly from the GWR. MI was calculated for each site using SPLASH v1.0 (Davis et al., 2017) based on daily values of precipitation, temperature and sunshine hours obtained using a mean-conserving interpolation of the monthly values of each. MI was then transformed to α_plant using the parametric Fu-Zhang formulation of the Budyko relationship (Supplement, Sect. S2). The climate space occupied by SMPDSv3 (Fig. S1 in the Supplement) samples a reasonable range of global climate space and therefore should provide robust reconstructions of climate changes under glacial conditions.

We use a global modern dataset for calibration of the pollen-climate relationships. The use of a global dataset, rather than region-specific training data, relies on the principle of phylogenetic niche conservatism (Harvey and Pagel, 1991; Qian and Ricklefs, 2004; Wang et al., 2025), which states that traits tend to remain constant over time and that the climatic niches of specific genera are also conservative (Harrison et al., 2025c). The use of a global dataset for calibration makes it possible to sample a large range of climates, and thus makes it more likely that the reconstructions of glacial climates are realistic and not confined to the limited climate range sampled in any one region in modern times (Turner et al., 2020).

2.1.2 Fossil pollen dataset

The Abrupt Climate Changes and Environmental Responses (ACER) database (Sánchez Goñi et al., 2017) was originally created to provide a source of pollen and charcoal data for Marine Isotope Stage 3 (MIS 3), which includes 93 records with sufficient resolution and dating control to detect sub-millennial scale variability. Many more records covering MIS 3 have become available since the compilation of the ACER database, such as the synthetic pollen databases available for Siberia (Cao et al., 2019, 2020) and China (Zhou et al., 2023) and the global Legacy 2 dataset (Li et al., 2025), which can substantially cover the spatial gaps in the original ACER database. We obtained these data from public sources or directly from the authors and used them to create an update: ACER2 (Harrison et al., 2025b), which contains 233 additional records covering some part or all of MIS 3 (note that the original ACER records are not included in ACER2 due to licensing issue). The two datasets are combined in our analyses to serve as the fossil pollen dataset (Supplement, Sect. S3) to reconstruct the past climates. We focus on the 279 records (253 terrestrial records and 26 marine records) between 50 and 30 ka (Fig. 1b; Table 1). The fossil pollen data are taxonomically harmonised to be consistent with the SMPDSv3.

Table 1 Site information of the fossil pollen records covering the interval 50–30 ka. The references with full citation are given in the Supplement, Sect. S3. lat is the latitude, lon is the longitude, elv is the elevation (unit: m). TERR means terrestrial record, MARI means marine record. n_due is the number of D–O events that should be found based on the time interval covered by the record. n_miss is the number of D–O events that were not identified. n_low is the number of D–O events missed because of low resolution of that part of the record. Some sites provide records in 50–30 ka but do not cover the intervals of the D–O events; some marine sites are too far from the land to extract GWR modern MTGR to apply CO₂ correction; these sites are all indicated by NA in n_due, n_miss and n_low.

2.2 Climate reconstruction method

We use fx-corrected Tolerance-weighted Weighted Averaging Partial Least Squares (fxTWA-PLS: Liu et al., 2020, 2023) regression to derive the pollen-climate relationships in the modern training dataset, and then apply these relationships to reconstruct past climates from the fossil pollen records (Fig. 2). fxTWA-PLS reduces the tendency of regression methods to compress reconstructions towards the centre of the sampled climate range by applying a sampling frequency correction to reduce the influence of uneven sampling of climate space and weighting the contribution of individual taxa according to their climate tolerances (Liu et al., 2020). Version 2 of fxTWA-PLS (fxTWA-PLSv2: Liu et al., 2023) uses P-splines smoothing to derive the frequency correction and applies this correction both in estimating the climate optima and tolerances, and in the regression itself, producing a further improvement in model performance compared to fxTWA-PLSv1 (Liu et al., 2020).

Figure 2 Flow chart showing the methodology.

We use fxTWA-PLSv2 here. The evaluation is made by comparing the reconstructions made using modern pollen data with modern climates using leave-out cross-validation, where one site at a time is randomly selected as a test site and sites that are both geographically close (within 50 km horizontal distance from the site) and climatically close (within 2 % of the full range of each climate variable in the dataset) are also removed from the training set, to prevent redundancy in the climate information from inflating the cross-validation goodness of fit, following Liu et al. (2020). This ensures that we are not just tuning to the training dataset, and that we can reconstruct climates even when the training set does not completely cover the climate to be reconstructed because there are gaps in the climate space. Performance is assessed using R² and RMSEP (root-mean-square error of prediction), and compression is assessed using linear regression of the leave-out cross-validated reconstructions against the climate variable. The last significant number of components (p ≤ 0.01) is selected to avoid overfitting due to the increase in the number of components. Reconstructions of MTCO, MTWA and α_plant are then made for every sample in each fossil record, using the last significant number of components. Sample-specific errors are estimated via bootstrapping (resampling the training set 1000 times) as described in Liu et al. (2020).

However, the low CO₂ at glacial period could lead to potential bias between reconstructed and actual plant-available moisture. Atmospheric CO₂ concentration has a direct impact on plant physiological processes, by modulating water-use efficiency (WUE), that is the ratio of carbon uptake to water loss through the stomata (Hatfield and Dold, 2019). The low CO₂ during the glacial period led to reduced water use efficiency (Farquhar, 1997; Gerhart and Ward, 2010; Prentice and Harrison, 2009). Strictly statistical reconstructions cannot take this into account since they are based on modern relationships between pollen assemblages and climate under recent CO₂ levels (Bartlein et al., 2011; Chevalier et al., 2020). Reconstructions based on inversion of a vegetation model (e.g. Garreta et al., 2010; Guiot et al., 2000; Izumi and Bartlein, 2016; Wu et al., 2007) implicitly account for the impact of CO₂ on vegetation composition, but the reconstructions are dependent on the reliability of the vegetation model and its sensitivity to CO₂ changes (Chevalier et al., 2020; Guiot et al., 2009). The actual conditions under low CO₂ should be wetter than the vegetation-based reconstructions of moisture variables (Prentice et al., 2017, 2022a). Prentice et al. (2022a) provide a correction method to account for variations on CO₂ based on combining eco-evolutionary optimality theory and experimental evidence on how the water-use efficiency as expressed by the ratio of leaf-internal to ambient CO₂ responds to variations in CO₂, as follows:

$$\begin{array}{}\text{(1)}& e\left({\mathrm{MTGR}}_{\mathrm{1}},{\mathrm{MI}}_{\mathrm{1}},c_{\mathrm{a}\mathrm{1}}\right)=e\left({\mathrm{MTGR}}_{\mathrm{0}},{\mathrm{MI}}_{\mathrm{0}},c_{\mathrm{a}\mathrm{0}}\right)\end{array}$$

where e is the ratio of water loss to CO₂ uptake, a function of the mean temperature of the growing season (MTGR), moisture index (MI) and atmospheric CO₂ concentration (c_a). For MTGR and c_a, the subscript “1” denotes the past value, and the subscript “0” denotes the modern value. MI₀ is the reconstructed uncorrected past value, MI₁ is the “true” past value (to be estimated). The equation means that the “true” MI under past atmospheric conditions should produce the same e with the reconstructed uncorrected MI under modern atmospheric conditions, i.e. those pertaining to the modern pollen calibration dataset.

We transfer our reconstructed past α_plant back to the uncorrected moisture index MI₀, and apply the CO₂ correction to obtain the actual moisture index MI₁, then transfer it to actual plant-available moisture α_{plant,corrected} (Figs. 2 and S2a, b). Past and modern values of CO₂ concentrations are taken from Bereiter et al. (2015), following Prentice et al. (2022a). Past MTGR values are inferred by sinusoidal interpolation of reconstructed MTCO and MTWA, assuming that the growing season corresponds to the period with temperatures > 0 °C. Modern MTGR values are obtained using a geographically-weighted regression (GWR) of climatological values (1961–1990) from the Climatic Research Unit Time-Series version 4.04 (CRU TS4.04; Harris et al., 2020) dataset averaged over the period 1961–1990, in order to correct for elevation differences between the CRU grid cells and the fossil pollen sites. The elevations of marine sites are set to 0 when applying GWR.

The CO₂ correction is implemented through the package COdos 0.0.2 (Prentice et al., 2022b) with one modification, as follows. We found when applying the correction in cases where the temperature reduction from modern was large (> 5 °C) that the use of different temperature values to calculate the stomatal sensitivity term (ξ) and the compensation point (Γ^∗) in the correction algorithm sometimes produced an unrealistically large countervailing effect due to the temperature difference alone. To avoid this problem, we calculate these physiological quantities (ξ and Γ^∗) using the mean of MTGR₁ and MTGR₀.

2.3 Age modelling of fossil records

Both the ACER and ACER2 database provide age models for each pollen record. However, the resolution of the individual records are variable (ranging from 57 years to 13 415 years) and these age models are often imperfectly aligned with the dating of D–O warming events as recorded in the Greenland ice core, and which have been shown to have a globally synchronous imprint through analysis of speleothem records (Adolphi et al., 2018; Corrick et al., 2020). To create a better alignment, we use dynamic time warping (DTW: Alshehri et al., 2019; Burstyn et al., 2021; Giorgino, 2009) to adjust the age scale for each individual record (Fig. 2). Dynamic time warping optimises the similarity between two sequences (one “query” and one “reference”) by stretching or compressing one sequence in the time dimension to match the other. It adjusts the age scale without influencing the variable values, thus retaining the original amplitude of change.

The LOVECLIM simulation (Menviel et al., 2014) (covering the interval 50–30 ka) is currently the only published simulation that has attempted to reproduce the specific timing and magnitude of successive D–O cycles. It is a coupled ocean-atmosphere-vegetation general circulation model of intermediate complexity. The model was spun up to equilibrium using an initial atmospheric CO₂ concentration of 207.5 ppm, orbital forcing appropriate for 50 ka BP, and an estimate of the 50 ka BP ice-sheet orography and albedo obtained from an off-line ice-sheet model simulation (Abe-Ouchi et al., 2007). After this initialization, the model was forced by time-varying changes in orbital parameters, atmospheric trace gas concentrations and ice-sheet configurations following Timm et al. (2008). In addition, meltwater pulses were added in the North Atlantic in such a way as to reproduce observed sea-surface temperature (SST) variations along the Iberian margin (Martrat et al., 2007). The simulations have proved adequate to capture at least the broad features of actual D–O events, and are generally consistent with the qualitative signals in Voelker (2002) compilation (Liu et al., 2022). We convert the age scale of LOVECLIM simulations to match the Antarctic Ice Core Chronology 2012 (AICC2012) time scale (Veres et al., 2013).

We treat the mean annual temperature (MAT) calculated as the average of MTCO and MTWA reconstructed from each individual fossil pollen record as the “query” time series, and find the corresponding grid cell (the location of the fossil pollen record) in the LOVECLIM simulations, and use the simulated MAT at this grid cell as the “reference” time series. We further divide each “query” and “reference” time series into discrete intervals using the mid-points between the start dates of each D–O warming event as recorded in the Greenland ice core (Wolff et al., 2010; converted into AICC2012 timescale), and normalize both time series in each interval to remove the influence of differences in absolute values and the amplitude of changes. Then we apply dynamic time warping to modify the time scale of the “query” to match the “reference” in each interval. The adjusted age model for each fossil record is then applied to the reconstructions of MTCO, MTWA, α_plant and α_{plant,corrected} from that record for subsequent analyses.

2.4 Assessment of regional climate changes during Greenland D–O warming events

The magnitude of climate change during the interval corresponding to each D–O warming event as registered in Greenland is calculated individually for each climate variable at each site. To avoid making an assumption about the sign of the climate change at a site, we use a third-order polynomial to fit the reconstructions during the interval from 300 years before to 600 years after the official start date corresponding to Greenland D–O warming for each event (Wolff et al., 2010; converted into AICC2012 timescale) to determine whether the change was positive or negative. We then find the ages where this polynomial curve reaches the minimum and maximum (t_{min,polynomial} and t_{max,polynomial}). Since the smoothed polynomial may underestimate or overestimate the amplitude of change, we use the reconstructions corresponding to t_{min,polynomial} and t_{max,polynomial} to obtain the changes (see Fig. S3 for illustration). Whether it's an increasing or decreasing signal depends on whether t_{min,polynomial} occurs before or after t_{max,polynomial}. The change of each climate variable (ΔV) is calculated as:

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }V=V_{\mathrm{end}}-V_{\mathrm{start}}\end{array}$$

where V_start is the reconstructed value at the start and V_end is the reconstructed value at the end of the event. The error of change (σ_ΔV) is calculated using the following equation assuming V_start and V_end are independent:

$$\begin{array}{}\text{(3)}& {\mathit{\sigma }}_{\mathrm{\Delta }V}=\sqrt{{\mathit{\sigma }}_{V_{\mathrm{end}}}^{\mathrm{2}}+{\mathit{\sigma }}_{V_{\mathrm{start}}}^{\mathrm{2}}}\end{array}$$

where ${\mathit{\sigma }}_{V_{\mathrm{start}}}$ is the sample-specific error of V_start and ${\mathit{\sigma }}_{V_{\mathrm{end}}}$ is the sample-specific error of V_end.

To obtain the relationships between changes in different climate variables, we use a maximum likelihood method to estimate the ratio of ΔMTCO to ΔMTWA and the ratio of Δα_{plant,corrected} to ΔMTWA to take account of the errors on both variables, following Liu et al. (2022).

As a measure of the accuracy of the DTW method to identify D–O events, we compare the number of identified events with the number of D–O events that should occur during the time covered by each record (Table 1). To assess whether events are missed in a particular record due to low sampling resolution, we examine the number of samples present in the 900-year interval covering the sampled D–O (i.e. 300 years before to 600 years after the official start date corresponding to Greenland D–O warming for each event), where low resolution is defined as ≤ 3 samples in this 900-year interval.

3 Results

fxTWA-PLS reproduces the modern climate reasonably well (Table 2; Fig. S4a, b). The performance is best for MTCO (R² = 0.74, RMSEP = 6.66, slope = 0.84) but is also good for MTWA (R² = 0.60, RMSEP = 3.63, slope = 0.72) and α_plant (R² = 0.63, RMSEP = 0.186, slope = 0.68). Assessment of the variance inflation factor scores shows that there is no problem of multicollinearity so that it is possible to reconstruct all three climate variables independently (Table S1).

Table 2 Leave-out cross-validation (with geographically and climatically close sites removed) using fxTWA-PLSv2, for mean temperature of the coldest month (MTCO), mean temperature of the warmest month (MTWA) and plant-available moisture (α_plant). P-splines smoothed fx was estimated using 200 bins. n is the number of components used; Avg.bias is the average bias; RMSEP is the root-mean-square error of prediction; ΔRMSEP% is the per cent change of RMSEP of the current number of components compared to using one component less, i.e. 100 × (RMSEP_n − RMSEP_n−1) $/$ RMSEP_n−1. The p value assesses whether using the current number of components represents a significant (p ≤ 0.01) difference over using one fewer component. To avoid over-fitting, the last significant number of components (i.e. the first insignificant number of components minus 1; p can become significant again after being insignificant with increasing n, but this is regarded meaningless) is selected for subsequent analyses and indicated in bold. The degree of overall compression is assessed by linear regression of the cross-validated reconstructions against the variable; b₁ and ${\mathit{\sigma }}_{b_{\mathrm{1}}}$ are the slope and the standard error of the slope, respectively. A slope (b₁) of 1 indicates no compression. NA means that no value could be calculated for ΔRMSEP% and p when n=1, since they are defined based on comparison with using one fewer component.

The use of dynamic time warping makes it possible to identify D–O events robustly (Fig. S5a–h; Tables 1 and S2). Some sites provide records in 50–30 ka but do not cover the intervals of the D–O events; some marine sites are too far from the land to extract GWR modern MTGR to apply the CO₂ correction. Across the remaining 179 sites which should have D–O events registered, we have identified 544 out of the 696 individual events (78 %). In the majority of cases where a D–O event should have been registered but could not be identified in an individual record (134 out of 152 cases), the resolution of that part of the record is extremely poor.

ΔMTCO is found to be significantly larger than ΔMTWA in the northern extratropics and tropics when considered across all D–O events and sites, indicating reduced seasonal contrast between winter and summer temperatures; ΔMTCO is found to be larger than ΔMTWA, but not significantly larger, in the southern extratropics (Fig. 3; Table 3). There is no globally consistent relationship between Δα_{plant,corrected} and ΔMTWA, although the positive relationship in the tropics is marginally significant (Fig. 4; Table 4).

Figure 3 Scatter plot of the change in mean temperature of the coldest month (ΔMTCO) versus the change in mean temperature of the warmest month (ΔMTWA) during individual Dansgaard–Oeschger (D–O) events at individual sites. The points are grouped into the northern extratropics (NET, north of 23.5° N), the tropics (TROP, between 23.5° N and 23.5° S) and the southern extratropics (SET, south of 23.5° S). The grey lines indicate ±1 error of the change.

Table 3 Maximum likelihood estimates of the relationship between the change in mean temperature of the coldest month (ΔMTCO) and the change in mean temperature of the warmest month (ΔMTWA) for the northern extratropics (NET, north of 23.5° N), tropics (TROP, between 23.5° N and 23.5° S) and southern extratropics (SET, south of 23.5° S). The intercepts were set to zero since both variables are changes.

Figure 4 Scatter plot of the change in CO₂-corrected plant-available moisture (Δα_{plant,corrected}) versus the change in mean temperature of the warmest month (ΔMTWA) during individual Dansgaard–Oeschger (D–O) events at individual sites. The points are grouped into the northern extratropics (NET, north of 23.5° N), the tropics (TROP, between 23.5° N and 23.5° S) and the southern extratropics (SET, south of 23.5° S). The grey lines indicate ±1 error of the change.

Table 4 Maximum likelihood estimates of the relationship between the change in CO₂-corrected plant-available moisture (Δα_{plant,corrected}) and the change in mean temperature of the warmest month (ΔMTWA) for the northern extratropics (NET, north of 23.5° N), tropics (TROP, between 23.5° N and 23.5° S) and southern extratropics (SET, south of 23.5° S). The intercepts were set to zero since both variables are changes.

The spatial patterns of ΔMTCO and ΔMTWA are generally consistent across multiple D–O events (Fig. 5), most noticeably that the largest warming occurs in Eurasia, while western North America and the southern extratropics are characterised by cooling. These patterns are also shown if only reconstructions where the change is twice the error of change are considered (Fig. S6), proving that the spatial patterns are robust to the choice of threshold. Nevertheless, both the magnitude of the changes and the spatial patterns vary between the D–O events (Fig. S7a, b). Strong events such as D–O 8 show more apparent changes (whether warming or cooling), as well as a strong antiphasing between northern and southern extratropical changes; while weak events such as D–O 9 show less apparent changes with almost no north-south antiphasing (Fig. 6).

Figure 5 Map showing the median change of site-based reconstructions for Dansgaard–Oeschger (D–O) events 5 to 12. The panels from top to bottom show the changes in mean temperature of the coldest month (ΔMTCO), mean temperature of the warmest month (ΔMTWA) and CO₂-corrected plant-available moisture (Δα_{plant,corrected}).

Figure 6 Map showing the change in mean temperature of the coldest month (ΔMTCO) and the change in mean temperature of the warmest month (ΔMTWA) for D–O 9 (a weak event) and D–O 8 (a strong event). The upper panel shows ΔMTCO, while the lower panel shows ΔMTWA.

The changes in plant-available moisture are less spatially coherent that the changes in temperature (Fig. 5). There is an increase in α_{plant,corrected} in some regions characterised by warming, for example, southeastern China and Japan; but there are mixed signals of drying and wetting in other regions characterised by warming, such as southern Europe. Furthermore, regions characterised by cooling, such as western North America and the southern extratropics, can also show both drying and wetting. Changes in Δα_{plant,corrected} also show more variability between D–O events than changes in temperature (Fig. S7d). We have applied a correction for low CO₂ values during the glacial period to plant-available moisture. The actual values (α_{plant,corrected}) are generally higher than the vegetation-based reconstructed values (α_plant) (Fig. S2a). However, the correction does not have a significant impact on the spatial patterns during D–O events (Figs. S2b and S7c).

4 Discussion and Conclusions

4.1 Comparison with previous reconstructions

We have presented a first attempt to map the spatial patterns of quantitative changes in seasonal temperatures and plant-available moisture during D–O events globally, using a consistent methodology. These analyses show that there is an anti-phasing between changes in the northern extratropics and the southern extratropics, with warming in the north and cooling in the south. The largest and most consistent warming during D–O events occurs in Eurasia. There is a significant difference between winter warming and summer warming in the northern extratropics, resulting in an overall reduction in seasonality. Site-based reconstructions from the Eifel region in central Europe, based on branched glycerol dialkyl glycerol tetraethers, indicate minimal temperature changes during summer (Zander et al., 2024) and thus support the idea that the D–O changes were driven by large changes in winter temperature. Zumaque et al., (2025) provide seasonal temperature and precipitation reconstructions for 12 of the sites from southern Europe (which are included in our fossil pollen records) but using the modern analogue technique as the reconstruction method and the Eurasian Modern Pollen Database version 2 (Davis et al., 2020) (EMPDv2; also included in SMPDSv3) as the modern training dataset. They show relatively stable summer temperatures but large change in MTCO through the MIS3 D–O events, consistent with our reconstructions (using a regression-based reconstruction method and a global modern training dataset) of a reduction in seasonality during warming events in the northern extratropics. We find no significant reduction in seasonality in the southern extratropics. Since only quantitative reconstructions of MAT (rather than MTCO and MTWA) are available from the southern extratropics (e.g. Fletcher and Thomas, 2010; Newnham et al., 2017), there is no independent confirmation of this result.

Qualitative interpretation of palaeo-records suggest that many regions are characterised by both warming and wetting, such as western Europe (Fletcher et al., 2010; Sánchez Goñi et al., 2008), eastern Europe (Fleitmann et al., 2009; Stockhecke et al., 2016), central Siberia (Grygar et al., 2006), and the Great Basin USA (Denniston et al., 2007; Jiménez-Moreno et al., 2010). Previous studies have also indicated drier conditions during D–O events, particularly in parts of the USA such as the Pacific Northwest (Grigg and Whitlock, 2002) and Florida (Grimm et al., 2006; Jiménez-Moreno et al., 2010). Our reconstructions show more mixed signals and that there is no globally consistent relationship between changes in temperature and moisture, either in regions characterised by warming or by cooling (Figs. 4 and 5).

Figure 7 Map showing the median change of LOVECLIM simulations over (ice-free) land for Dansgaard–Oeschger (D–O) events 5 to 12. The upper panel shows the change in mean temperature of the coldest month (ΔMTCO), and the lower panel shows the change in mean temperature of the warmest month (ΔMTWA).

4.2 Global training dataset vs. local training dataset

We have used a global pollen dataset for calibration of the pollen-climate relationships. In general, reconstructions of glacial climates have used region-specific data sets (e.g. Dugerdil et al., 2021, 2025; Newnham et al., 2017; Wei et al., 2021; Zumaque et al., 2025). Herzschuh et al. (2023) made this explicit in their reconstructions of northern hemisphere climate over the past 30 000 years, by restricting the modern training data to within a 2000 km radius of individual fossil sites. The use of a region-specific raining data set can be justified on the grounds that it produces better statistics for the modern-day relationship between pollen abundance and specific climate variables. Nevertheless, as pointed out by Chevalier et al. (2020), an important issue is that the modern calibration data set has a span that adequately samples the climate space experienced in the past. The use of a global dataset for calibration makes it possible to sample a larger range of climates, and specifically to reconstruct climates that might be very different from the modern range in that region. For example, reconstructions of past European climate (Fig. S8) based on region-specific training dataset would yield less extreme winter temperatures than reconstructed using the global training data set. Although the trend and spatial pattern might not be influenced greatly, the amplitude of change might be underestimated.

The use of a global dataset, rather than region-specific training data, relies on the principle of phylogenetic niche conservatism (Harvey and Pagel, 1991; Qian and Ricklefs, 2004; Wang et al., 2025), which states that traits tend to remain constant over time. This also applies to the climate niche (Crisp and Cook, 2012; Jiang et al., 2023; Peterson, 2011; Wiens et al., 2010; Wiens and Graham, 2005) as evidenced by disjunct distributions of taxa across different continents (Yin et al., 2021). Niche conservatism underpins the fact that the modern distribution of specific genera can be predicted using climate-pollen relationships developed from other regions (e.g. Huntley et al., 1989). However, the use of a global dataset can create issues because of inconsistencies in taxonomic resolution between regions. The necessity for treating all species of Quercus as a single taxon, despite the fact that evergreen and deciduous species may occupy distinct climate niches in some regions, is a consequence of this. However, we have shown (Supplement, Sect. S4) that this has little impact on our reconstructions – largely because the climatic distinction that would be conveyed through separating deciduous and evergreen Quercus is also registered by the presence of other taxa. Although the use of a global training dataset for climate reconstructions has not been a common practice, it addresses the need to ensure that the modern training data adequately represents past climate conditions and also facilitates making reconstructions for sites from regions with limited modern pollen data.

4.3 Targets for model evaluation

The reconstructions in this paper can be used as targets for model evaluation, specifically the two transient D–O experiments planned for the next phase of the Palaeoclimate Modelling Intercomparison (see Malmierca-Vallet et al., 2023 for the experimental protocol). The first of these experiments is a baseline simulation starting at 34 ka, a time with low obliquity, moderate MIS3 greenhouse gas values, and an intermediate ice sheet configuration, which appears to be most conducive to generating D–O like behaviour in climate models. The second experiment involves the addition of freshwater, to examine whether this is necessary to precondition a state conducive to generating D–O events. The anti-phasing in reconstructed temperature changes between the northern and southern hemispheres is a general feature of climate model experiments. Most models show larger warming in winter than in summer in the northern hemisphere (e.g. Flückiger et al., 2008; Izumi et al., 2023; Van Meerbeeck et al., 2011), which is also consistent with our reconstructions. However, the cooling in western North America during D–O warming events in our reconstructions is not a feature of all climate model simulations.

Models generally show an intensification of the northern hemisphere monsoons during D–O events (e.g. Izumi et al., 2023; Menviel et al., 2020), but there is less consistency about changes in plant-available moisture in the extratropics. Our reconstructions show an increase in α_{plant,corrected} in southeastern China and Japan (Fig. 5). Although α_{plant,corrected} is not a direct reflection of summer precipitation, these changes are consistent with enhanced northern hemisphere monsoons during D–O warming events, as shown by speleothem records from the Caribbean (Warken et al., 2019) and speleothem and pollen records from Asia (Fohlmeister et al., 2023; Wang et al., 2001; Zorzi et al., 2022). However, there is more spatial variability and mixed signals in our reconstructions.

The LOVECLIM model was used as a reference to adjust the age scale in the reconstructions using MAT, but this does not preclude comparison of the seasonal temperatures. Here we approximate the winter-season temperature as MTCO and summer-season temperature as MTWA, since monthly temperatures are not available (only seasonal temperatures are available) in LOVECLIM. The general spatial pattern of simulated changes in MTCO and MTWA (Fig. 7) is consistent with the reconstructions, with largest warming in Eurasia, and cooling in the southern extratropics. The simulated changes are strong during D–O 8 but weak during D–O 9 (Fig. S9a, b), again as shown by the reconstructions. However, there are important differences. For example, simulated changes generally have smaller amplitude than shown by the reconstructions, and the cooling over western North America is generally only in winter, while the reconstructions show cooling over this region in both seasons. The relationship between ΔMTCO and ΔMTWA is also different (Fig. 8; Table 5): the simulated ΔMTCO is shown to be significantly larger than ΔMTWA in the northern extratropics, but significantly smaller than ΔMTWA in the southern extratropics, a contrast that is not so marked in the reconstructions. This comparison illustrates the usefulness of the reconstructions for model evaluation and to investigate the mechanisms that may not be adequately captured by current models.

Figure 8 Scatter plot of the change in mean temperature of the coldest month (ΔMTCO) versus the change in mean temperature of the warmest month (ΔMTWA) during individual Dansgaard–Oeschger (D–O) events at individual (ice-free) land grids simulated by the LOVECLIM model, using the same way to identify changes as the reconstructions. The points are grouped into the northern extratropics (NET, north of 23.5° N), the tropics (TROP, between 23.5° N and 23.5° S) and the southern extratropics (SET, south of 23.5° S). The grey lines indicate ±1 error of the change.

Table 5 Maximum likelihood estimates of the relationship between the change in mean temperature of the coldest month (ΔMTCO) and the change in mean temperature of the warmest month (ΔMTWA) over (ice-free) land in LOVECLIM simulations, for the northern extratropics (NET, north of 23.5° N), tropics (TROP, between 23.5° N and 23.5° S) and southern extratropics (SET, south of 23.5° S). The intercepts were set to zero since both variables are changes.

4.4 Implications of the use of dynamic time warping

Identifying D–O events in pollen records is often problematic, particularly in regions where warming (especially if accompanied by dryer conditions) leads to a reduction (or an hiatus) in sedimentation as reflected in the variable resolution of the available pollen records (e.g. Camuera et al., 2022; Pini et al., 2022; Sinopoli et al., 2019; Wei et al., 2021). The use of dynamic time warping goes some way to improving the identification of potential D–O events. However, it precludes the calculation of a rate of change in climate. Thus, we have focused here on the magnitude of the changes during specific warming events. It is also likely that some of the variability in the reconstructed changes between different D–O events reflects imperfect identification of specific events because of the comparatively modest resolution of the records.