Seasonal temperature and moisture changes in interior semi‐arid Spain from the last interglacial to the Late Holocene

Abstract

The El Cañizar de Villarquemado pollen record covers the last part of MIS 6 to the Late Holocene. We use Tolerance-Weighted Averaging Partial Least Squares (TWA-PLS) to reconstruct mean temperature of the coldest month (MTCO) and growing degree days above 0°C (GDD₀) and the ratio of annual precipitation to annual potential evapotranspiration (MI), accounting for the ecophysiological effect of changing CO₂ on water-use efficiency. Rapid summer warming occurred during the Zeifen-Kattegat Oscillation at the transition to MIS 5. Summers were cold during MIS 4 and MIS 2, but some intervals of MIS 3 had summers as warm as the warmest phases of MIS 5 or the Holocene. Winter temperatures declined from MIS 4 to MIS 2. Changes in temperature seasonality within MIS 5 and MIS 1 are consistent with insolation seasonality changes. Conditions became progressively more humid during MIS 5, and MIS 4 was also humid, although MIS 3 was more arid. Changes in MI and GDD₀ are anti-correlated, with increased MI during summer warming intervals. Comparison with other records shows glacial-interglacial changes were not unform across the circum-Mediterranean region, but available quantitative reconstructions are insufficient to determine if east-west differences reflect the circulation-driven precipitation dipole seen in recent decades.

INTRODUCTION

The modern climate of the Mediterranean region is influenced by high-pressure systems in summer and westerly storm tracks in winter, giving rise to a highly seasonal climate with dry summers, relatively wet winters, and strong seasonal temperature contrasts. Nevertheless, the climate of the region is not homogeneous and the pattern of climate changes across the region is complex (Lionello, 2012). There is a contrast between the reduction in annual precipitation during recent decades (1975–2015) across the Iberian Peninsula and in northern Italy, for example, and the generally positive trends in precipitation in Greece, the Balkans, and Turkey (Deitch et al., 2017). This east-west dipole is also characteristic of rainfall trends over the twentieth century (Zittis, 2018) and during the Holocene (Peyron et al., 2017). It is highly likely that climate changes on longer timescales, specifically orbital timescales in response to change in insolation, and on glacial-interglacial timescales in response to changes in ice-sheet volume (Magri and Tzedakis, 2000; Rohling et al., 2013), are equally complex.

The El Cañizar de Villarquemado palaeolake is located in Teruel province in the semi-arid interior of the Iberian Peninsula (40.49°N, 1.29°W, 985 m a.s.l.), and provides a pollen record that stretches from the end of the penultimate glacial (Marine Isotope Stage 6, MIS 6) through the last interglacial (MIS 5), the last glaciation (MIS 4–2), and into MIS 1 and the Holocene (Valero-Garcés et al., 2019; González-Sampériz et al., 2020). The length of the record and the quality of the age model make El Cañizar de Villarquemado an important site for documenting climate changes in the western Mediterranean, on both orbital (reflecting changes in precession, obliquity, and eccentricity on insolation) and glacial-interglacial (reflecting ice sheet and CO₂ changes) timescales. The only other record from the western Mediterranean that spans this length of time is Padul (Pons and Reille, 1988; Camuera et al., 2018, 2019), in southern Spain. Despite some similarities between the Padul and El Cañizar de Villarquemado pollen records, there are distinctive features that reflect their different climatic settings and sensitivities. There are also notable differences between the pollen record from El Cañizar de Villarquemado and other long records from the Mediterranean region, including the classic site of Tenaghi Phillippon in Greece (Tzedakis et al., 2006; Milner et al., 2013), Lago di Monticchio in central Italy (Allen et al., 2000; Allen and Huntley, 2009, 2018; Martín-Puertas et al., 2014), and Lake Ohrid on the North Macedonia/Albania border (Wagner et al., 2017; Sinopoli et al., 2018, 2019), which suggest regional differentiation of long-term climate changes across the Mediterranean.

Inference of changes in precipitation from these long pollen records is complicated by large changes in atmospheric CO₂ concentration [CO₂] on glacial-interglacial timescales. [CO₂] directly affects photosynthesis and water use (Ehleringer et al., 1997; Farquhar, 1997; Prentice and Harrison, 2009). Increasing [CO₂] allows plants that use the standard C₃ pathway of photosynthesis (including temperate grasses and forbs, and nearly all shrubs and trees) to assimilate more carbon while losing less water, implying an increase in water-use efficiency (Bramley et al., 2013). Conversely, at low [CO₂], C₃ plants assimilate less carbon per unit water supply, necessarily resulting in reduced tree cover and a shift towards more drought-tolerant plants (including C₄ plants). Shifts in the balance of steppe plants and tree taxa have been interpreted in terms of moisture-related climate variables (e.g., Sinopoli et al., 2018; Camuera et al., 2019), but this neglects the negative influence of low [CO₂] on water-use efficiency. As a result, inferred or statistically reconstructed moisture variables are expected to indicate drier conditions than actually occurred (Prentice et al., 2017).

In this paper, we present a quantitative reconstruction of three bioclimatic variables: winter temperature (mean temperature of the coldest month, MTCO), growing-season warmth (growing degree days above a base level of 0°C, GDD₀), and a moisture index (the ratio of annual precipitation to annual potential evapotranspiration, MI) using the pollen record from El Cañizar de Villarquemado and Tolerance-Weighted Averaging Partial Least Squares regression (TWA-PLS) (Liu et al., 2020) that are calibrated using modern pollen (Harrison, 2019) and climate (Harrison, 2020) datasets for Europe and northern Eurasia. These three variables were chosen to represent the independent influences of extreme winter cold, integrated summer warmth, and plant-available moisture on the success of different plants. We apply a method of correcting MI to take account of the direct physiological influence of [CO₂] on water-use efficiency. Although El Cañizar de Villarquemado is well dated, there are several intervals with poor pollen preservation, so we focus on the long-term evolution of climate rather than on abrupt climate changes. Finally, we compare our reconstructions with pollen-based reconstructions from the wider Mediterranean region.

METHODS

Modern Pollen and Climate Data

The modern pollen dataset, SPECIAL Modern Pollen Data Set (SMPDS), consists of relative abundance records of the 247 most important pollen taxa from 6458 terrestrial sites from Europe and northern Eurasia (Harrison, 2019; Turner et al., 2020; Fig. 1). The SMPDS taxon list includes several layers of specificity: species, genera, sub-families, and families. Some taxa are represented at <10 sites and are therefore unlikely to provide a reliable basis for climate reconstructions. We used the 195 taxa that occur >10 times in the SMPDS as the initial basis for the climate reconstructions (Supplementary Table 1). We ran additional analyses after removing two taxa (Poaceae, Polypodiales) that displayed anomalous peaks in abundance that appeared to be related to changes in the sedimentary environment rather than climate. The final climate reconstructions are therefore based on 193 taxa.

Figure 1. Map showing the location of El Cañizar de Villarquemado and the distribution of modern pollen samples used in the calibration. The background map shows the climatological (1961–1990) mean temperature of the coldest month (MTCO, °C). Modern pollen data from Harrison (2019) and modern climate data from the CRU CL 2.0 gridded dataset (New et al., 2002).

The modern climate at each of the SMPDS pollen sites was taken from Harrison (2020). This dataset provides three bioclimatic variables: (a) mean temperature of the coldest month (MTCO), (b) growing degree days above a baseline of 0°C (GDD₀), and (c) moisture index (MI). MI is the ratio of annual precipitation to annual potential evapotranspiration, calculated using climatological data (mean monthly temperature, precipitation, and fractional sunshine hours) derived from the CRU CL 2.0 gridded dataset of modern (1961–1990) surface climate at 10 arc minute resolution (~18 km) (New et al., 2002) and geographically weighted regression (GWR) to correct for elevation differences between each pollen site and the centroid of the corresponding grid cell. MI values were square-root transformed before analysis to emphasize the stronger vegetational differentiation at low MI.

Canonical correspondence analysis (CCA) (ter Braak, 1986; Legendre and Legendre, 2012) was used to perform a constrained ordination of the modern pollen data in response to the bioclimatic variables. CCA was implemented with the vegan package (Oksanen et al., 2017) in R (v3.3.1). We examined the variance inflation factors (VIFs) to determine whether the predictors showed multi-collinearity. The significance of the CCA model was assessed with an ANOVA-like permutation test.

Fossil pollen data

The El Cañizar de Villarquemado palaeolake is located in the Jiloca basin in the semi-arid region of north-eastern Spain (Fig. 2). The site is occupied today by a wetland and cultivated land. The surrounding vegetation is dominated by evergreen or marcescent trees (e.g., Quercus ilex subsp. rotundifolia, Q. coccifera, Q. faginea) and xerophytic shrubs (e.g., Rhamnus lycioides, Genista scorpius, Ephedra fragilis, Thymus vulgaris, T. zygis). Juniperus (J. thurifera, J. phoenicea, J. communis, J. sabina, J. oxycedrus) and Pinus (P. sylvestris, P. pinaster, P. nigra) occur at higher elevations. The sequence provides a paleoenvironmental record spanning the interval from the last part of MIS 6 to the Late Holocene (Valero-Garcés et al., 2019; González-Sampériz et al., 2020). The age model from MIS 3 onwards is well constrained. Although uncertainties in the earlier part of the record are larger, the dates for key transitions in this part of the record are consistent with accepted estimates at the 90% confidence level (Valero-Garcés et al., 2019). Sedimentation rates are low during the initial part of the record and increase from the beginning of MIS 2 onwards. Intervals with poor pollen preservation occur between 16,000–22,300, 31,200–37,500, 43,100–50,100, and 87,900–93,800 cal yr BP. The average pollen sampling interval is ca 300 yr, increasing to ca. 140 yr during MIS 1. The pollen data were re-expressed in terms of the taxon groupings used in the modern pollen dataset for the application of TWA-PLS. The data used for TWA-PLS are available from http://dx.doi.org/10.17864/1947.219 (Harrison et al., 2019).

Figure 2. Vegetation map of El Cañizar de Villarquemado area showing dominant tree species, surface area of the basin and location of the core site.

TWA-PLS

The modern bioclimatic and pollen data were used to develop pollen-climate transfer functions independently for MTCO, GDD₀, and √MI using tolerance-weighted-averaging partial least squares regression (TWA-PLS) (Liu et al., 2020). TWA-PLS is an extension of weighted-averaging partial least squares regression (WA-PLS) (ter Braak and Juggins, 1993) that reduces the tendency of WA-PLS to produce climate reconstructions compressed towards the middle of the climate range represented by the modern training dataset. In TWA-PLS, taxon abundances are weighted by the inverse square of their modern climate tolerances. The theory behind TWA-PLS (a) assumes that the abundances of each taxon are unimodal with respect to the climate variable considered, and (b) follows a multinomial distribution. Under these conditions, TWA-PLS approximately minimises the difference between theoretical and actual abundances as assessed by the likelihood criterion. The performance of the calibration models was assessed through leave-one-out cross-validation. The number of components used in each model was estimated through a randomization t-test on the results of this cross-validation (van der Voet, 1994). We selected the significant component with the lowest root mean square error of prediction (RMSEP), but only if there was a significant improvement in RMSEP relative to a lower number of components because including any more components results in over-fitting so that model predictive power decreases (ter Braak et al., 1993). We checked that the final transfer functions had a high R² for prediction, and a low maximum bias.

We propagated the calibration errors into the down-core reconstructions by assessing the calibration error associated with each pollen sample using the “v1” option in rioja (Juggins, 2017). The errors arising from uncertainties in the modern pollen dataset tend to be large, so this approach emphasizes temporal differences in reconstruction uncertainty (results based on combining the two sets of uncertainties are shown in Supplementary Information).

Correcting for Changing [CO₂] Concentration

Prentice et al. (2017) developed a procedure to correct for direct effect of [CO₂] on plant physiological processes. The procedure (Appendix 1) requires the specification of [CO₂] and mean annual temperature (MAT). We used the ice-core [CO₂] record (Bereiter et al., 2015), smoothed using a LOESS smoothing spline with a span of 0.1. We estimated MAT for each site from the reconstructed MTCO and GDD₀ at that site (Appendix 2). This calculation also allowed us to estimate the mean temperature of the warmest month (MTWA) (Appendix 2), which we then be used to generate a time series of temperature seasonality (MTWA−MTCO). We then converted the down-core reconstructions of √MI (and their maximum and minimum values) at El Cañizar de Villarquemado to MI, and applied the [CO₂] correction to these MI values.

RESULTS

CCA showed a strong correlation between species abundance and the three climate variables in the modern pollen dataset, with correlations of 0.83, 0.61, and 0.47 respectively (Table 1). The VIF scores for each bioclimatic variable were low (<6), and the CCA indicates that they each have an independent contribution to explaining the variation in abundance. This was confirmed by the ANOVA-like permutation test, which showed that the bioclimatic variables and the three variability axes are all significantly different from one another (Table 1). As a further check, partial CCA was carried out using each predictor in turn (with the other two predictors as covariates). This confirmed that each predictor has a highly significant effect that it is independent of the others (Table 2).

Table 1. Summary statistics for the first three axes of CCA in the modern pollen dataset. The analysis was based on 6458 sites, 195 taxa, and three bioclimatic variables: mean temperature of the coldest month (MTCO, °C), growing degree days above 0°C (GDD₀), and the moisture index (MI). Summary statistics for the ANOVA-like permutation test (999 permutations) are also shown.

Table 2. Summary statistics for the ANOVA-like permutation test (999 permutations) of partial CCA in the modern pollen dataset. Partial CCA analysis used each predictor in turn, with the other two predictors as covariates.

For TWA-PLS regression based on the full complement of 195 taxa, we used results from component 2 for MTCO and MI and component 1 for GDD₀ because these are the significant results (95% confidence level) with the lowest RMSE and highest R². The R² values are 0.70, 0.66, and 0.52 for MTCO, GDD₀, and √MI, respectively (Supplementary Table 2). Nevertheless, close examination of the down-core reconstructions revealed anomalous peaks in reconstructed MI, particularly at the end of MIS 5. These correspond to samples with unusually high values of either Poaceae or Polypodiales (Fig. 3), and where the sedimentary record indicates variability in environmental conditions and that the basin was occupied by wetlands. Under these conditions, the anomalously high peaks of Poaceae likely correspond to reeds, while the anomalously high peaks of Polypodiales probably represent inwash. Both Poaceae and Polypodiales were therefore removed from the final TWA-PLS model (Table 3), reducing the total number of taxa considered from 195 to 193. This alteration made no change to the number of components or the goodness-of-fit of the model (R² values are 0.70, 0.68, and 0.51 for MTCO, GDD₀, and √MI, respectively), but made the reconstructions of MI for the anomalous samples less extreme, more comparable to reconstructed values in adjacent samples, and therefore more plausible (Fig. 4). The removal of these taxa had no significant effect on the MTCO and GDD₀ reconstructions (Supplementary Fig. 1). We checked whether particularly high or low values in the temperature reconstructions were a result of anomalous characteristics of the pollen assemblages, specifically whether there was evidence of pollen degradation (as measured by the fraction of indeterminable grains) or the samples were characterized by low biodiversity (as measured by the N2 index) (Hill, 1973). There was no correlation between the abundance of indeterminable grains and anomalously high or low reconstruction values. However, depauperate samples tended to produce more extreme temperature values than adjacent more diverse samples (Supplementary Fig. 2). We therefore excluded samples with N2 <2 from the final reconstructions. Excluding samples with N2 <3 had little further effect on the reconstructions, but increased the patchiness of the reconstructions by removing a large number of samples.

Figure 3. Simplified stratigraphic and pollen diagram from El Cañizar de Villarquemado (a more complete diagram is given in González-Sampériz et al., 2020). The first column shows changes in inferred sedimentary environments, including the alternation between lacustrine and non-lacustrine conditions. The simplified pollen diagram shows the changing abundance of conifers, mesophyte trees, Mediterranean, and steppe taxa. We also show the changing abundance of aquatics and Polypodiales. The final column shows the biodiversity index (Hill's N2). The transitional phase at the beginning of MIS 5 corresponding to the Zeifen-Kattegat Oscillation is labelled Zeifen in the column showing the marine isotopic stages.

Figure 4. The impact of removing Poaceae and Polypodiales from the taxon set on reconstructions of moisture index (the ratio of annual precipitation to annual potential evapotranspiration, MI) during MIS 4 and MIS 5a. The reconstructed √MI has been re-expressed as MI, but no CO₂ correction has been applied. The red line (without P/P) shows the values of MI once Poaceae and Polypodiales are removed. Removing these two taxa reduces anomalous peaks, where they were particularly abundant, but has little impact on the reconstructions for the rest of the core.

Table 3. Leave-one-out cross-validated predictions of the tolerance-weighted averaging partial least squares (WA-PLS) regression models used for the climate reconstructions. The P values were derived from a randomization t-test. The final model is based on 193 taxa, omitting Poaceae and Polypodiales. Selected components in the final model are marked in bold.

The final reconstructions (Fig. 5, Supplementary Tables 3, 4) show an increase in both summer (GDD₀) and winter (MTCO) temperature between 130 and 127 ka. Both summer and winter temperature generally declined through MIS 5—although there are fluctuations, they do not correspond exactly with the chronological boundaries of sub-stages within MIS 5 (Fig. 5). Furthermore, the changes in summer and winter temperature are not in phase. Minimum winter temperatures occurred earlier than minimum summer temperature in MIS 5e, so that winter temperatures were already increasing while summer temperatures continued to decrease after ca. 120 ka. In contrast, during MIS 5a, warming in summer occurred around the same time as winter cooling, such that the temperature seasonality was significantly enhanced between ca. 80 and 70 ka. There was no pronounced cooling, either in summer or winter, at the transition into the glacial (Fig. 5). The record from both MIS 3 and MIS 2 is not continuous, and the available samples may show the response to millennial-scale events; thus, it is difficult to characterize the trends in temperature. However, MIS 3 appears to have been somewhat warmer than both MIS 2 and MIS 4 in summer (Supplementary Table 4). MIS 3 was not significantly warmer in winter than other intervals during the glacial, and there was a gradual decline in winter temperature from MIS 4 through MIS 3 to MIS 2 (Supplementary Tables 3, 4). MIS 1 was characterized by warming in both summer and winter, although the reconstructions show considerable variability superimposed on this trend.

Figure 5. Reconstructed mean temperature of the coldest month (MTCO, °C), growing degree days above a base level of 0°C (GDD₀), and CO₂ corrected moisture index (the ratio of annual precipitation to annual potential evapotranspiration, MI). The reconstructions are based on 193 taxa, after the removal of Poaceae and Polypodiales. Only samples with a Hill's N2 biodiversity index >2 are plotted. The marine isotope stages (MIS) and substages are shown by vertical dotted lines and labelled; we also show the transition interval between MIS 6 and MIS 5e. Red dots indicate the modern climate at the site. The error bars on the reconstructions, shown by shading around the mean value, are sample-specific errors (v1 in rioja).

The implied increase in temperature seasonality during MIS 5e, 5c, and 5a compared to during MIS 5d and 5b corresponds to the differences in seasonality in insolation between these stages (Fig. 6), primarily driven by summer insolation. Higher seasonality during the early part of MIS 1 also corresponds to increased seasonality in insolation compared to the present day (Fig. 6). Changes in seasonal differences in insolation across the glacial were muted. However, temperature seasonality was larger during MIS 3 and MIS 2 than during the interglacial intervals. Intervals of increased temperature seasonality during MIS 3 and MIS 2 cannot be explained by changes in the seasonality of insolation.

Figure 6. The correlation of the temperature seasonality and insolation. The black line in the top plot is the difference between reconstructed mean temperature of the coldest month and the mean temperature of the warmest month calculated based on MTCO and GDD₀ (Appendix 2); we have applied a 1000 year time step with a span of 0.1 to draw attention to the major trends. The orange line in the top plot is the difference between July and January insolation in W m⁻² at 40.49°N (the latitude of El Cañizar de Villarquemado). The bottom panel shows mid-monthly insolation anomalies (compared to present) in W m⁻² at 40.49°N.

Conditions became progressively more humid from MIS 5e through to MIS 5c, while conditions were generally humid, but variable, during MIS 5a (Fig. 5, Supplementary Tables 3, 4). MIS 4 was also relatively humid, while MIS 3 was the most arid phase reconstructed during the glacial. However, the difference in reconstructed MI between MIS 4 or MIS 2 and MIS 3 is not large. This reflects the fact that [CO₂] decreased throughout the glacial, so that the impact of the CO₂ correction becomes larger from MIS 4 through MIS 3 and into MIS 2 (Fig. 7). The influence of changing [CO₂] is most marked in dry climates (Fig. 8), which is why this effect has such an important influence on our reconstructions of MI. The highest values of MI are reconstructed during the deglaciation; there was a pronounced increase in aridity during the Holocene. The reconstructed changes in MI are anti-correlated with changes in GDD₀ (r = −0.68), with decreased MI during intervals of summer warming. This suggests that the changes in MI may have been driven more by changes in evapotranspiration than by changes in precipitation.

Figure 7. Changes in [CO₂] and their impact on reconstructed moisture index (the ratio of annual precipitation to annual potential evapotranspiration, MI). The upper panel shows the changes in [CO₂] (from Bereiter et al., 2015), with a LOESS smoothing with a span of 0.1, and the lower panel shows reconstructed MI with (orange line) and without (blue line) taking account of the [CO₂] correction.

Figure 8. Scatter plot showing the impact of the [CO₂] correction on the reconstructed moisture index (MI). The colored dots represent the implied change in the reconstructions, grouped according to level of the actual [CO₂] at that time (in ppm).

The direct impact of low [CO₂] on tree growth accounts for the relatively low abundance of mesophytic and Mediterranean trees during the glacial. Reconstructed MI is somewhat higher than today. Nevertheless, intervals of increasing MI during the glacial—as in the first part of MIS 2—are marked by a decrease in steppe taxa and an increase in conifers. The increase in Mediterranean taxa during the later Holocene, at a time when reconstructed MI indicates progressive drying, reflects the reconstructed strong insolation-driven winter warming through the Holocene because Mediterranean trees occur at areas with intermediate levels of water availability, but are sensitive to winter temperature (Harrison et al., 2010).

Several abrupt changes are shown in the reconstructions, particularly during MIS 5a and in the glacial period. Some of these (Supplementary Fig. 3) occur at times that correspond to D-O events, as identified from Greenland ice core records (Dansgaard et al., 1984; Steffensen et al., 2008), including D-O 20 (72.28–70.28 cal ka) and D-O 19 (76.4–74 cal ka) in MIS 5a, and D-O 9 (40.11–39.81 cal ka) and D-O 8 (38.17–36.57 cal ka) in MIS 3 (for dating see Wolff et al., 2010). Heinrich Stadial 2 (26.45–24.25 cal ka) (Sanchez Goñi and Harrison, 2010) also corresponds to an interval of year-round cooling in our reconstructions. Gaps in the pollen record, and poor dating resolution in some parts of the record, preclude identification of other D-O and Heinrich events. The limitations of the age model mean that even the correspondence between abrupt changes registered in El Cañizar de Villarquemado and the Greenland D-O events is uncertain. Nevertheless, where abrupt D-O like events are registered, they were characterized by increased seasonality, which may contribute to the high seasonality recorded during some periods during the glacial (Fig. 6).

DISCUSSION AND CONCLUSION

The El Cañizar de Villarquemado record is characterized by rapid warming in winter temperature of >5°C and an increase in growing-season warmth of >4000 degree days over a period of ca. 2–3 kyr during the transition from MIS 6 to MIS 5e. Although there are fluctuations, there is an overall decline in both summer and winter temperatures through MIS 5. However, there is a long interval with poor pollen preservation during MIS 5b, which limits our ability to gain a complete picture of the evolution of climate during this interval. Temperature reconstructions for MIS 4, 3, and 2 are not significantly lower than the end of MIS 5, but this might be because the coldest intervals occurred during the intervals of low pollen preservation in MIS 3 and 2. The Younger Dryas interval is marked by relatively cold summers (Supplementary Table 3). There was a gradual warming in both summer and winter through the Holocene. The broad-scale changes in moisture are coherent with changes in GDD₀, with warmer summer intervals characterized by drier conditions and colder summers by wetter conditions. The MI reconstructions apparently indicate that the whole of the past ca. 130 kyr was wetter than today, although it is likely that the intervals of poor pollen preservation (for which we were unable to make reconstructions) were at least as arid as today. Rapid millennial-scale changes in temperature and moisture are superimposed on these longer-term trends, though not all D-O events can be identified in the record.

Many features of the El Cañizar de Villaquemado record are also shown in other quantitative reconstructions from the Mediterranean region. Both the Monticchio (Allen and Huntley, 2009) and the Lake Ohrid (Sinopoli et al., 2019) records show rapid warming at the transition from MIS 6 to MIS 5e. This warming occurs over a longer period in the Monticchio (ca. 5 kyr) and Ohrid (ca. 7 kyr) records. The apparent rapidity of this transition at El Cañizar de Villarquemado could reflect age-model uncertainties and/or regional idiosyncrasies in the vegetation composition. Differences between the sub-stages of MIS 5 are more pronounced in the Monticchio and Ohrid records than at El Cañizar de Villarquemado. Comparison of the modern analogue and WA-PLS reconstructions at Lake Ohrid shows that modern-analogue reconstructions (and, by implication, the response-surface approach used at Monticchio) tend to produce stronger fluctuations, which might contribute to explaining the more muted variability at El Cañizar de Villarquemado. However, the pronounced cold, dry interval registered in Monticchio and Ohrid during MIS 5b corresponds to an interval of low pollen preservation in Villarquemado. This, and the fact that the El Cañizar de Villarquemado site lies at the warmer and drier end of the climate gradient across the Mediterranean, could contribute to the apparent differences between the sites.

The coldest interval at Monticchio during MIS 2 is also represented by an hiatus in Villarquemado. The Younger Dryas was characterized by cooler summers and wetter conditions at Villarquemado, but only a small decrease in winter temperature. The wetter conditions and the muted winter temperature response are consistent with the record from Monticchio (Allen and Huntley, 2009). However, the Holocene record from Monticchio is very different, showing summer cooling and relatively stable moisture levels after 10 ka. Thus, while there are some similarities between the available quantitative records from the Mediterranean, they each show distinctive features. These might reflect differences in the size and geomorphic setting of the three lakes, and hence the pollen source area (Sugita, 1984; Prentice, 1985; Prentice et al., 1988) and, by implication, the geographic region represented by the climate reconstructions, but it is more likely that they reflect the pattern of climate changes across the region, as also indicated by qualitative information from other paleoclimate records (e.g., Harrison et al., 2002; Roberts et al., 2004, 2008; Lechleitner et al., 2018).

The temperature record at El Cañizar de Villarquemado shows intervals of enhanced seasonality during MIS 5 and MIS 1, largely, but not entirely, driven by changes in GDD₀. We have shown that there is a good correlation between these intervals of enhanced seasonality and orbitally forced changes in summer insolation during interglacial intervals. Temperature seasonality during the glacial was higher than during the interglacials and cannot be explained by insolation, although the within-glacial trends in seasonality do follow the changes in insolation. Izumi et al. (2013) showed that increased temperature seasonality at the last glacial maximum (LGM) compared to present, particularly in the northern extratropics, was a response to year-round forcing. They subsequently argued (Izumi et al., 2014) that this was a consequence of changes in land-ocean and high-low latitude temperature gradients rather than an independent temperature response to the large-scale forcing. Rehfeld et al. (2018), using a network of terrestrial and marine records, showed a fourfold decrease globally in temperature seasonality from the LGM into the Holocene. While this global decrease is larger than indicated at El Cañizar de Villarquemado, Rehfeld et al. (2018) pointed out that the change is larger at higher northern latitudes and negligible in the tropics. Thus, the higher seasonality during the glacial than during interglacials reconstructed at El Cañizar de Villarquemado appears to be a robust feature. Enhanced seasonality at El Cañizar de Villarquemado is also associated with intervals of abrupt warming that correlate with D-O events (e.g., D-O 9). However, enhanced seasonality during the D-O events appears to have been driven primarily by changes in winter temperature. Changes in MI on both orbital and millennial timescales are generally anti-correlated with changes in GDD₀ presumably because increased summer temperature and/or increased length of the growing season led to increased evapotranspiration and hence reduced MI.

There are gaps in the palynological record from El Cañizar de Villarquemado because of intervals of poor pollen preservation during MIS 5b, MIS 3, and MIS 2 (Fig. 3). The sedimentological record suggests that these were arid intervals, characterized by alluvial fan rather than lacustrine deposition, and oxidation of the sediments. A speleothem record from El Pindal (Moreno et al., 2010) also shows hiatuses in formation during MIS 2, consistent with our interpretation that the depositional hiatus at El Cañizar de Villarquemado is indicative of aridity. Similar situations have been identified in other palynological sequences from the region during arid intervals (Valero-Garcés et al., 2000, 2004; González-Sampériz et al., 2004, 2005; Vegas-Villarubía et al., 2013). Hyper-arid periods are a problem for pollen preservation and, while further work may improve the pollen record at El Cañizar de Villarquemado, it is unlikely that we will be able to obtain a quantitative record of climate during such intervals. Nevertheless, El Cañizar de Villarquemado provides a relatively complete and well-dated record from continental Iberia (González-Sampériz et al., 2010, 2020; Moreno et al., 2012; Valero-Garcés et al., 2019) and it is important to document changes in the drier part of the circum-Mediterranean region.

All reconstruction methods that use modern pollen-climate relationships are sensitive to the choice of training datasets, vegetation diversity, and the potential absence of analogue assemblages (Gavin et al., 2003; Jackson and Williams, 2004; Bartlein et al., 2011). While the training dataset that we have used provides a good coverage of Europe, Morocco, and the Middle East, it lacks modern pollen from much of northern Africa. However, it contains more than 6000 samples and includes samples from very cold and very warm environments to allow reconstruction of climate both much warmer and much colder than today (Turner et al., 2020). Analysis of the GAMs for individual taxa also shows that they have ecologically plausible relationships with climate variables (Wei et al., 2020). Furthermore, uncertainty estimates on the reconstructions, calculated using a method that propagates the calibration errors from the modern training dataset into the down-core reconstructions by assessing the calibration error associated with the particular assemblage in an individual pollen sample, are generally rather small. We have shown that pollen preservation issues, as indicated by intervals when the number of indeterminable grains was higher than average, do not affect our reconstructions. However, our analyses show that intervals of very low biodiversity are often characterized by more extreme values than other intervals. We have taken this into account by screening the down-core samples and excluding samples with very low diversity.

Much of the discussion about the lack of modern analogues has focused on interpretation of assemblages of species that are not found together today (Overpeck et al., 1985; Jackson and Williams, 2004; Williams and Shuman, 2008). However, one important non-analogue situation that is ignored in all previous statistical reconstructions is the impact of [CO₂] different from today on plant assemblages, with important effects on apparent climatic moisture inferred from pollen data (Prentice et al., 2017; Cleator et al., 2019). Taking the impact of [CO₂] into account in our reconstructions reduces the variability of MI during glacial intervals. Prentice et al. (2017) showed that this correction produced a reconciliation of apparently contradictory interpretations of pollen and geomorphic data for hydroclimatic changes in southeastern Australia at the LGM. High lake levels coincident with steppic vegetation are a feature of several glacial records from the Mediterranean region, including the Padul record (Camuera et al., 2018, 2019) and sequences from northeastern Iberia, such as the playa lakes of the Central Ebro Basin (Valero-Garcés et al., 2000, 2004; González-Sampériz et al., 2004, 2005). Camuera et al. (2019) argued that this reflects low evapotranspiration because of the low temperatures. However, much larger changes in temperature than implied by the El Cañizar de Villarquemado record would be required (e.g., Prentice et al., 1992). A simpler explanation is that the abundance of xerophytes and the reduction in trees reflect the ecophysiology impact of low [CO₂].

The El Cañizar de Villarquemado reconstructions provide a detailed picture of the response of western Mediterranean climate and vegetation to changes in external forcing in this sensitive region over a long period of time. It would be useful to generate quantitative reconstructions from other long records because preliminary comparisons with Lago di Monticchio and Lake Ohrid indicate some complexity in the response of climate to changes in external forcing within the circum-Mediterranean region. However, given that the largest impact of glacial-interglacial changes in atmospheric circulation is likely to be on precipitation and plant-available moisture, it will be important to take account of the impact of changing [CO₂] in these reconstructions.