Holocene climates of the Iberian Peninsula: pollen-based reconstructions of changes in the west–east gradient of temperature and moisture

The Iberian Peninsula is characterized by a steep west–east moisture gradient at present, reflecting the dominance of maritime influences along the Atlantic coast and more Mediterranean-type climate further east. Holocene pollen records from the Peninsula suggest that this gradient was less steep during the mid-Holocene, possibly reflecting the impact of orbital changes on circulation and thus regional patterns in climate. Here, we use 7214 pollen samples from 117 sites covering part or all of the last 12 000 years to reconstruct changes in seasonal temperature and in moisture across the Iberian Peninsula quantitatively. We show that there is an increasing trend in winter temperature at a regional scale, consistent with known changes in winter insolation. However, summer temperatures do not show the decreasing trend through the Holocene that would be expected if they were a direct response to insolation forcing. We show that summer temperature is strongly correlated with plant-available moisture (α), as measured by the ratio of actual evapotranspiration to equilibrium evapotranspiration, which declines through the Holocene. The reconstructions also confirm that the west–east gradient in moisture was considerably less steep during the mid-Holocene than today, indicating that atmospheric circulation changes (possibly driven by orbital changes) have been important determinants of the Holocene climate of the region.

1 Introduction

The Iberian Peninsula is characterized by a steep west–east gradient in temperature and moisture today, reflecting the dominance of maritime influences along the Atlantic coast and more Mediterranean-type climate further east. Projections of future climate change suggest that the region will become both warmer and drier, but nevertheless show that this west–east differentiation is maintained (Andrade et al., 2021a). The changes in temperature are projected to be larger and the occurrence of extreme temperature episodes more frequent in the south-central and eastern parts of Iberia than in Atlantic coastal areas (Carvalho et al., 2021). Similar gradients are seen in future projections of precipitation change, with the largest reductions in precipitation in the south-central region (Andrade et al., 2021b). However, the stability of these west–east gradients during the Holocene has been questioned. In particular, the west–east gradient in moisture appears to have been less pronounced during the mid-Holocene (8–4 ka) when cooler summers and wetter conditions in the Atlantic zone (e.g. Martínez-Cortizas et al., 2009; Mauri et al., 2015) coincided with the maximum development of mesophytic vegetation further east and south (Aranbarri et al., 2014, 2015; Carrión et al., 2010, 2009; González-Sampériz et al., 2017).

However, much of the evidence for Holocene climates of the Iberian Peninsula is based on qualitative interpretations of vegetation changes, generally interpreted as reflecting changes in moisture availability (Morellón et al., 2018; Ramos-Román et al., 2018; Schröder et al., 2019). These records are extensive and they seem to indicate fairly complex spatial patterns of change. Kaufman et al. (2020) provides quantitative reconstructions of summer and winter temperature in their compilation of Holocene climate information, but there are only five terrestrial sites from the Iberian Peninsula. Iberia was also included in the quantitative pollen-based reconstructions of European climate through the Holocene in Mauri et al. (2015), which is an update of Davis et al. (2003). However, the geographical distribution of sites included is uneven and a large fraction of the records were from the Pyrenees and the Cantabrian mountains, with additional clustering of sites in coastal regions. Thus, the inferred patterns of climate over most of the central part of the peninsula are therefore largely extrapolated. Tarroso et al. (2016) have provided reconstructions of summer and winter temperature and mean annual precipitation since the Last Glacial Maximum for the Iberian Peninsula by using modern species distribution data to develop climate probability distribution functions (PDFs) and applying these to 31 fossil records. However, although they identified trends in precipitation during the Holocene, the temperature reconstructions do not seem to be reliable since they show no changes through time (9–3 ka), either for the Iberian Peninsula as a whole or for individual sub-regions, which is in contradiction to the other reconstructions. The current state of uncertainty about Holocene climate changes in Iberia is further exacerbated because quantitative reconstructions of summer temperature made at individual sites using chironomid data (Muñoz Sobrino et al., 2013; Tarrats et al., 2018) are not consistent with reconstructed summer temperatures based on pollen for the same sites.

We used the tolerance-weighted weighted averaging partial least squares regression with a sampling frequency correction (fxTWA-PLS) method introduced by Liu et al. (2020) as an improvement of the widely used weighted averaging partial least squares (WA-PLS; ter Braak and Juggins, 1993) method for reconstructing past climates from pollen assemblages. As presented in detail by Liu et al. (2020), this method is a more complete implementation of the theory underlying WA-PLS because it takes greater account of the climatic information provided by taxa with more limited climatic ranges and also applies a sampling frequency correction to reduce the impact of uneven sampling in the training data set. Liu et al. (2020) showed that fxTWA-PLS does indeed provide better reconstructions than WA-PLS. Here, we have further modified the algorithm implementing fxTWA-PLS, achieving an additional gain in performance. In the algorithm published by Liu et al. (2020), sampling frequencies were extracted from a histogram. In the modified algorithm, they are estimated using P-splines smoothing (Eilers and Marx, 2021), which makes the estimates almost independent of the chosen bin width (see Appendix A for details). In addition, the modified method applies the sampling frequency correction at two separate steps – the estimation of optima and tolerances and the regression step – a measure intended to produce more stable results. The modified method produces both improved R² values and reduced compression and maximum bias in reconstructed climate variables (see Table A1, Figs. A1 and A2). We will return to this point in the Discussion.

We have used this improved method to reconstruct Holocene climates across Iberia, and re-examined the trends in summer and winter temperature and plant-available moisture, using a new and relatively comprehensive compilation of pollen data (Shen et al., 2022) with age models based on the latest radiocarbon calibration curve (IntCal20; Reimer et al., 2020). We explicitly test whether there are significant differences in the west–east gradient of moisture and seasonal temperatures through time. We then analyse the relationships between the changes in the three climate variables and how trends in these variables are related to external climate forcing. These analyses allow us to investigate whether the west–east gradient in moisture was less steep during the mid-Holocene and explore what controls the patterns of climate change across the region.

2 Methods

Multiple techniques have been developed to make quantitative climate reconstructions from pollen (see reviews in Bartlein et al., 2011; Chevalier et al., 2020; Salonen et al., 2011). Modern analogue techniques (MAT; Overpeck et al., 1985) tend to produce rapid shifts in reconstructed values corresponding to changes in the selection of the specific analogue samples, although this tendency is less marked in the conceptually analogous response surface technique (Bartlein et al., 1986). Regression-based techniques, including weighted averaging methods such as weighted averaging partial least squares (WA-PLS; ter Braak and Juggins, 1993), do not produce step changes in the reconstructions but suffer from the tendency to compress the reconstructions towards the central part of the sampled climate range. However, this tendency can be substantially reduced by accounting for the sampling frequency (fx) and the climate tolerance of the pollen taxa present in the training data set (fxTWA-PLS; Liu et al., 2020). Machine-learning and Bayesian approaches have also been applied to derive climate reconstructions from pollen assemblages (Peyron et al., 1998; Salonen et al., 2019). However, comparison of fxTWA-PLS with the Bayesian model BUMPER (Holden et al., 2017) shows that fxTWA-PLS performs better in capturing the climate of the modern training data set from Europe (Liu et al., 2020).

Although fxTWA-PLS has clear advantages over other quantitative reconstructions techniques, there is still a slight tendency towards compression. We have therefore made a further modification to the approach as described in Liu et al. (2020). In the original version of fxTWA-PLS, the fx correction is applied as a weight with the form of $\mathrm{1}/{\text{fx}}^{\mathrm{2}}$ in the regression (step 7 in Table 1 in Liu et al., 2020). Here (see Appendix A), we make a further modification of fxTWA-PLS by (a) applying the fx correction separately in both the taxon calculation and the regression (step 2 and step 7 in Table 1 in Liu et al., 2020) as a weight with the form of $\mathrm{1}/\text{fx}$ and (b) applying P-splines smoothing (Eilers and Marx, 2021) in order to reduce the dependence of the fx estimation on bin width. The modified version further reduces the biases at the extremes of the sampled climate range.

There are no generally accepted rules as to the choice of variables for palaeoclimate reconstruction. No systematic comparison of these choices has been made. However, it is widely understood that plant taxon distributions reflect distinct, largely independent controls by summer temperature, winter temperature and moisture availability (see e.g. Harrison et al., 2010). Therefore, in common with many other studies (Cheddadi et al., 1997; Jiang et al., 2010; Peyron et al., 1998; Wei et al., 2021; Zhang et al., 2007), we have chosen bioclimatic variables that reflect these independent controls, with mean temperature of the coldest month (MTCO) to represent winter temperature, mean temperature of the warmest month (MTWA) to represent summer temperature and α, an estimate of the ratio of actual evapotranspiration to equilibrium evapotranspiration, to represent plant-available moisture. We choose not to use mean annual air temperature (MAAT) because it is a composite of summer and winter conditions; and we prefer to use an index of effective moisture availability (our estimate of α being one such index) to mean annual precipitation (MAP), whose significance for plant function depends strongly on potential evaporation (a function of temperature and net radiation). Our calculation of α takes account of this dependence. Growing degree days above a baseline of 0 ^∘C (GDD₀) would be a possible alternative to MTWA as an expression of summer conditions but is most relevant as a predictor of “cold limits” of trees in cool climates, whereas MTWA better reflects the high-temperature stress on plants in Mediterranean-type climates.

We used the modified version of fxTWA-PLS to reconstruct these three climate variables. The individual and joint effects of MTCO, MTWA and α were tested explicitly using canonical correspondence analysis (CCA). The modified version further reduces the biases at the extremes of the sampled climate range, while retaining the desirable properties of WA-PLS in terms of robustness to spatial autocorrelation (fxTWA-PLS; Liu et al., 2020).

The modern pollen training data set was derived from the SPECIAL Modern Pollen Data Set (SMPDS; Harrison, 2019). The SMPDS consists of relative abundance records from 6458 terrestrial sites from Europe, northern Africa, the Middle East and northern Eurasia (Fig. S1 in the Supplement) assembled from multiple different published sources. The pollen records were taxonomically standardized, and filtered (as recommended by Chevalier et al., 2020) to remove obligate aquatics, insectivorous species, introduced species and taxa that only occur in cultivation (see Table S1 in the Supplement for the list). Taxa (mainly herbaceous) with only sporadic occurrences were amalgamated to higher taxonomic levels (genus, sub-family or family) after ensuring consistency with their distribution in climate space. As a result of these amalgamations, the SMPDS contains data on 247 pollen taxa. For our analysis, we use the 195 taxa that occur at more than 10 sites.

Modern climate data at each of the sites in the training data set were obtained from Harrison (2019). This data set contains climate reconstructions of MTCO, growing degree days above a baseline of 0 ^∘C (GDD₀) and a moisture index (MI), defined as the ratio of annual precipitation to annual potential evapotranspiration. The climate at each site was obtained using geographically weighted regression (GWR) of the CRU CL v2.0 gridded data set of modern (1961–1990) surface climate at 10 arcmin resolution (New et al., 2002) in order to (a) correct for elevation differences between each pollen site and the corresponding grid cell and (b) make the resulting climate independent of the resolution of the underlying data set. The geographically weighted regression used a fixed bandwidth kernel of 1.06^∘ (∼ 140 km) to optimize model diagnostics and reduce spatial clustering of residuals relative to other bandwidths. The climate of each pollen site was then estimated based on its longitude, latitude and elevation. The MTCO and GDD₀ were taken directly from the GWR and MI was calculated for each pollen site using a modified code from SPLASH v1.0 (Davis et al., 2017) based on daily values of precipitation, temperature and sunshine hours, again obtained using a mean-conserving interpolation of the monthly values of each. For this application, we used MTCO directly from the data set but calculated MTWA from MTCO and GDD₀ based on the relationship between MTCO, MTWA and GDD₀ given in Appendix 2 of Wei et al. (2021). We derived α from MI following Liu et al. (2020). The modern training data set provides records spanning a range of MTCO from −42.4 to 14.8 ^∘C, of MTWA from 4.2 to 33.5 ^∘C and of α from 0.04 to 1.25 (Figs. 1 and S1).

Table 1 Details of the fossil pollen sites used. The fossil pollen data from the Iberian Peninsula were compiled by Shen et al. (2022) and obtained from https://doi.org/10.17864/1947.000369.

The fossil pollen data from the Iberian Peninsula were compiled by Shen et al. (2022) and the data set was obtained from Harrison et al. (2022). The taxonomy used by Shen et al. (2022) is consistent with that employed in the SMPDS. Shen et al. (2022) provides consistent age models for all the records based on the IntCal20 calibration curve (Reimer et al., 2020) and the BACON Bayesian age-modelling tool (Blaauw et al., 2021; Blaauw and Christeny, 2011) using the supervised modelling approach implemented in the ageR package (Villegas-Diaz et al., 2021). We excluded individual pollen samples with large uncertainties (standard error larger than 100 years) on the ages attributed in the new age model. As a result, the climate reconstructions are based on a fossil data set of 7384 pollen samples from 117 records covering part or all of the last 12 000 years (Fig. 2), with 42 individual records provided by the original authors, 73 records obtained from the European Pollen Database (EPD, http://www.europeanpollendatabase.net, last access: 10 March 2023) and 2 records from PANGAEA (https://www.pangaea.de/, last access: 10 March 2023). Details of the records are given in Table 1. The average temporal resolution of these records is 101 years. We then excluded a few samples where the reconstructed values of α exceed the natural limit of 0 and 1.26. Finally, 7214 samples from 117 records are used for the analyses of the climate reconstructions. Summer and winter insolation are also calculated using the PAST software based on the age and latitude of each sample (Hammer et al., 2001).

Variance inflation factor (VIF) scores were calculated for both the modern climates and the climates reconstructed from fossil pollen records in order to avoid multicollinearity problems and thus guarantee that the climate variables (MTCO, MTWA, α) used here represent independent features of the pollen records.

In addition to examining the reconstructions for individual sites, we constructed composite curves for the Iberian Peninsula as a whole. The composite curves provide a way of comparing the relationship between trends in the reconstructed climate changes and insolation changes. The curves were constructed after binning the site-based reconstructions using ± 500-year bins. We used 1000 bootstrap resamplings of the reconstructed climate values in each ± 500-year bin to avoid the influence of a single value or a single site on the mean climate value in this bin, and used the standard deviation of the 1000 values to represent the uncertainty of the mean climate value. We constructed linear regression plots to examine the longitudinal and elevational patterns in the reconstructed climate variables and assessed the significance of differences in these trends through time compared to the most recent bin (0.5 ka ± 500 years) based on p values, with the customary threshold of 0.05. We then compared the climate trends with changes in summer and winter insolation.

Table 2 Leave-out cross-validation (with geographically and climatically close sites removed) fitness of the modified version of fxTWA-PLS for the mean temperature of the coldest month (MTCO), mean temperature of the warmest month (MTWA) and plant-available moisture (α), with P-splines smoothed fx estimation, using bins of 0.02, 0.02 and 0.002, respectively; n is the number of components used; Avg.bias is the average bias, while Max.bias is the maximum absolute bias and Min.bias is the minimum absolute bias; RMSEP is the root-mean-square error of prediction; and ΔRMSEP is the per cent change rate of RMSEP, which is (RMSEP_n − RMSEP_n−1)$/$RMSEP_n−1 converted into percentage; when n=1, RMSEP₀ is the RMSEP of the null model. p assesses whether using the current number of components is significantly different from using one component less, which is used to choose the last significant number of components (indicated in bold) to avoid over-fitting. The degree of overall compression is assessed by linear regression of the cross-validated reconstructions onto the climate variable; b₁ and b₁.se are the slope and the standard error of the slope, respectively. The closer the slope (b₁) is to 1, the less the overall compression is.

Table 3 Canonical correspondence analysis (CCA) result of modern and fossil-reconstructed MTCO, MTWA and α. The summary statistics for the ANOVA-like permutation test (999 permutations) are also shown. VIF is the variance inflation factor, Df is the number of degrees of freedom, χ² is the constrained eigenvalue (or the sum of constrained eigenvalues for the whole model), F is significance and Pr (> F) is the probability. The CCA plots can be found in Fig. S11.

3 Results

The modified version of fxTWA-PLS reproduces the modern climate reasonably well (Table 2). The performance is best for MTCO (R² 0.75, RMSEP 4.70, slope 0.91) but is also good for α (R² 0.68, RMSEP 0.16, slope 0.78) and MTWA (R² 0.57, RMSEP 3.47, slope 0.71). The correlations between pollen records and each of the three bioclimate variables, as assessed by CCA, were strong for both modern climate data and fossil reconstructions (Table 3). The variance inflation factor (VIF) scores are all less than 6, so there are no multicollinearity problems (Table 3) (Allison, 1994), making it possible to reconstruct all three climate variables independently based on pollen data. Furthermore, the taxa that contribute most strongly to reconstructing colder/warmer or wetter/drier climates show predictable patterns consistent with their known ecological preferences (Table S2).

Figure 1 Climate space represented by mean temperature of the coldest month (MTCO), mean temperature of the warmest month (MTWA) and plant-available moisture as represented by α, an estimate of the ratio of actual evapotranspiration to equilibrium evapotranspiration. The grey points show climate values for a rectangular area (21^∘ W–150^∘ E, 29–82^∘ N) enclosing the SMPDS, derived from the Climate Research Unit CRU CL 2.0 database (New et al., 2002). The black points show climate values of the SMPDS. The red points show climate values of the Iberian Peninsula region in the SMPDS.

Figure 2 Map showing the location of the 117 fossil sites in the Iberian Peninsula used for climate reconstructions. Sites lower than 1000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ are shown as squares and sites higher than 1000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ are shown as triangles. The base maps show modern (a) mean temperature of the coldest month (MTCO), (b) mean temperature of the warmest month (MTWA) and (c) plant-available moisture as represented by α, an estimate of the ratio of actual evapotranspiration to equilibrium evapotranspiration.

Figure 3 Reconstructed anomalies in climate at individual sites through time. The sites are grouped into high- (> 1000 m) and low- (< 1000 m) elevation sites and organized from west to east. Grey cells indicate periods or longitudes with no data. The individual plots show the anomalies in reconstructed (a, d) mean temperature of the coldest month (MTCO), (b, e) mean temperature of the warmest month (MTWA) and (c, f) plant-available moisture as represented by α, an estimate of the ratio of actual evapotranspiration to equilibrium evapotranspiration. The anomalies are expressed as deviations of the mean value in each bin (± 500 years) from the most recent bin (0.5 ka ± 500 years) at each site.

Figure 4 Reconstructed composite changes (anomalies to 0.5 ka) in (a) mean temperature of the coldest month (MTCO), (b) mean temperature of the warmest month (MTWA) and (c) plant-available moisture as represented by α, through the Holocene compared to changes in (d) winter and (e) summer insolation for the latitude of the Iberian Peninsula, using ± 500 years as the bin. The black lines show mean values across sites, with vertical line segments showing the standard deviations of mean values using 1000 bootstrap cycles of site resampling.

Table 4 Assessment of the significance of anomalies to 0.5 ka through time with longitude and elevation. The slope is obtained by linear regression of the anomaly onto the longitude or elevation; p is the significance of the slope (bold parts: p < 0.05); x₀ is the point where the anomaly is 0 in the linear equation, which indicates longitude or elevation where the anomaly changes sign.

Winters were generally colder during the early to mid-Holocene than present, as shown by the coherent patterns of reconstructed anomalies at individual sites (Fig. 3a and d). Here, “present” means the most recent pollen bin (0.5 ka ± 500 years). The composite curve also shows a general increase in winter temperatures through time (Fig. 4a), consistent with the trend in winter insolation (Fig. 4d). The composite curve shows that it was ca. 4 ^∘C cooler than today at 11.5 ka, but temperatures were only ca. 0.5 ^∘C cooler than present after 2.5 ka. Winter temperatures today increase from north to south and are also affected by elevation; these patterns are still present in the Holocene reconstructions, but there is no spatial differentiation between western and eastern Iberia in the anomalies (Table 4, Fig. S2). The similarity of the changes compared to present geographically is consistent with the idea that the changes in winter temperature are driven by changes in winter insolation.

Summers were somewhat hotter than present in the west and cooler than present in the east during the early and mid-Holocene, as shown by the reconstructed anomalies at individual sites (Fig. 3b and e). This west–east difference could not arise if the changes in summer temperatures were a direct reflection of the insolation forcing (Fig. 4e). Indeed, the composite curve shows relatively little change in MTWA (Fig. 4b), confirming that there is no direct relationship to insolation forcing (Fig. 4e).

Figure 5 Changes in the west–east gradient of plant-available moisture as represented by anomalies in α relative to 0.5 ka at individual sites through the Holocene. The red lines show the regression lines. The shades indicate the 95 % confidence intervals of the regression lines

There is a strong west–east gradient in α at present (Fig. 2), with wetter conditions in the west and drier conditions in the east. However, the reconstructed anomalies at individual sites (Fig. 3c and f) suggest that the west was drier and the east was wetter than present in the mid-Holocene, resulting in a flatter west–east gradient. The west–east gradient is significantly different from the present between 9.5–3.5 ka (Fig. 5, Table 4), implying stronger moisture advection into the continental interior during the mid-Holocene. The change in gradient is seen in both high- and low-elevation sites (Fig. S3). There is also significant change in α with elevation between 9.5–4.5 ka (Table 4, Fig. S4).

Figure 6 The relationship between mean temperature of the warmest month (MTWA) and plant-available moisture as represented by α (a) in the modern climate data set and (b) in the Holocene reconstructions.

Summer temperatures are strongly correlated with changes in α, both in terms of spatial correlations in the modern data set at a European scale and in terms of spatial and temporal correlations in the fossil data set from the Iberian Peninsula (Fig. 6). The patterns of reconstructed anomalies in MTWA and α at individual sites are also coherent (Fig. 3b, c, e, and f), showing drier conditions and hotter summers in the west and wetter conditions with cooler summers in the east during the early to mid-Holocene. The west–east gradient in MTWA was significantly different from the present between 9.5 and 3.5 ka except 8.5 ka (Table 4, Fig. S5), roughly the interval when the gradient in α was also significantly different from the present. Again, the change in the west–east gradient is registered at both high- and low-elevation sites (Fig. S6). However, there is no significant change in MTWA with elevation except at 8.5 and 7.5 ka (Table 4, Fig. S7).

4 Discussion

The modified version of fxTWA-PLS (fxTWA-PLS2) (Tables 2 and A1) shows improved performance compared to the previous version (fxTWA-PLS1). Cross-validation R² values are higher for MTCO and MTWA but almost unchanged for α. The maximum bias is decreased for all three variables, especially for MTCO. The compression problem is also reduced for MTCO (b₁ increases from 0.82 to 0.91) and MTWA (b₁ increases from 0.69 to 0.71), while it remains roughly the same for α. The overall performance statistics thus show substantial improvements for MTCO and MTWA, while they show little change for α. However, Fig. A1 shows that “non-physical” reconstructions beyond the natural limits of α (0–1.26) are greatly reduced, especially at the lower limit. There are also fewer outliers in Figs. A1 and A2 for all three variables. Thus overall, the modified version further reduces the reconstruction biases, especially at the extremes of the sampled climate range. This improvement probably occurs because of the separate application of $\mathrm{1}/\text{fx}$ correction during both the calculation of optima and tolerances of taxa and during the regression step – instead of applying an overall weight of $\mathrm{1}/{\text{fx}}^{\mathrm{2}}$ at the regression step, which can result in some extreme values (with low sampling frequency) being weighed too strongly and appearing as outliers.

Figure 7 Comparison between reconstructed composite changes in climate anomalies. The first column represents this paper, the second column represents Mauri et al. (2015), the third column represents Kaufman et al. (2020) and the fourth column represents Tarroso et al. (2016). The composite curves from this paper and Kaufman et al. (2020) are calculated from individual reconstructions, using anomalies to 0.5 ka and a bin of ± 500 years (time slices are 0.5, 1.5, …, 11.5 ka). The composite curves from Mauri et al. (2015) are converted directly from the gridded time slices which are provided with anomalies to 0.1 ka and a bin of ± 500 years (time slices are 1, 2, …, 12 ka). The composite curves from Tarroso et al. (2016) are also converted directly from the gridded time slices provided, with anomalies to 0.5 ka and a bin of ± 500 years (time slices are 3, 4, …, 12 ka). Note that Tarroso et al. (2016) applied a smoothing to the data such that the plots in their paper do not show the excursion in MTWA at 8 ka. In all of the plots, the black lines show mean values across sites, with vertical line bars showing the standard deviation of mean values using 1000 bootstrap cycles of site/grid resampling.

Figure 8 Simulated mean values of mean temperature of the coldest month (MTCO), mean temperature of the warmest month (MTWA) and mean daily precipitation in the Iberian Peninsula between 8 and 0 ka, smoothed using 100-year bins. Here, BP means before 1950 AD. The black lines represent the Max Planck Institute (MPI) Earth system model simulations, the red lines represent the Alfred Wagner Institute (AWI) Earth system model simulations, the blue lines represent the Institut Pierre-Simon Laplace Climate Model (IPSL-CM5) TR5AS simulations and the gold lines represent the Institut Pierre-Simon Laplace Climate Model (IPSL-CM6) TR6AV simulations. The four simulations were forced by evolving orbital parameters and greenhouse gas concentrations. The four models have different spatial resolutions, with the finest resolution being 1.875^∘ × 1.875^∘ (AWI, MPI) and the coarsest resolution being 1.875^∘ × 3.75^∘ (IPSL-CM5, TR5AS).

The fxTWA-PLS2 reconstructions show that there was a gradual increase in MTCO over the Holocene, both for most of the individual sites represented in the data set (Fig. 3) and for Iberia as a whole (Fig. 4). Colder winters in southern Europe during the mid-Holocene (6 ka) are a feature of many earlier reconstructions (e.g. Cheddadi et al., 1997; Wu et al., 2007). A general warming trend over the Holocene is seen in gridded reconstructions of winter season (December, January, February) temperatures as reconstructed using the modern analogue approach by Mauri et al. (2015), although there is somewhat less millennial-scale variability in these reconstructions (Fig. 7). Nevertheless, their reconstructions show a cooling of 3 ^∘C in the early Holocene, comparable in magnitude to the ca. 4 ^∘C cooling at 11.5 ka reconstructed here. They show a gradual warming trend through the Holocene but the differences from the present are very small (ca. 0.5 ^∘C) after 2 ka, again consistent with our reconstructions of MTCO. Quantitative reconstructions of winter temperature for the five terrestrial sites from the Iberian Peninsula in the Kaufman et al. (2020) compilation all show a general trend of winter warming over the Holocene, but the magnitude of the change at some of the individual sites is much larger (ca. 10 ^∘C) and there is no assessment of the uncertainty on these reconstructions. The composite curve of Kaufman et al. (2020) shows an increasing trend in MTCO through the Holocene, although with large uncertainties (Fig. 7). In contrast to the consistency of the increasing trend in MTCO during the Holocene between our reconstructions and those of Mauri et al. (2015) and Kaufman et al. (2020), there is no discernible trend in MTCO during the Holocene reconstruction of Tarroso et al. (2016). Indeed, there is no significant change in their MTCO values after ca. 9 ka, either for the peninsula as a whole (Fig. 7) or for any of the four sub-regions they considered. Our reconstructed trend in winter temperature is consistent with the changes in insolation forcing at this latitude during the Holocene; it is also consistent with transient climate model simulations (Braconnot et al., 2019; Carré et al., 2021; Dallmeyer et al., 2020; Parker et al., 2021) of the winter temperature response to changing insolation forcing over the late Holocene in this region (Figs. 8 and S8). Thus, we suggest that changes in winter temperatures are a direct consequence of insolation forcing.

We have shown that there is no overall trend in MTWA during the Holocene (Fig. 4). According to our reconstructions, summer temperatures fluctuated between ca. 0.5 ^∘C above or below modern temperature. The lack of coherent trend in MTWA is consistent with the gridded reconstructions of summer (June, July, August) temperatures in the Mauri et al. (2015) data set and also with the five terrestrial sites from Iberia included in the Kaufman et al. (2020) data set. However, the patterns shown in the three data sets are very different from one another. Mauri et al. (2015) suggest the early Holocene was colder than today, and although temperatures similar to today were reached at 9 ka, most of the Holocene was characterized by cooler summers. Kaufman et al. (2020), however, showed warmer conditions during the early Holocene than at present, although they also show cooler conditions during the later Holocene. The differences between the three data sets could reflect differences in the reconstruction methods, or differences in the number of records used and in the geographic sampling. However, given the fact that all three data sets show similar trends in winter temperature, the lack of coherency between the data sets for MTWA points to there not being a strong, regionally coherent signal of summer temperature changes during the Holocene. Tarroso et al. (2016) also showed no significant changes in MTWA after ca. 9 ka (Fig. 7).

The chironomid record from Laguna de la Roya covers the late glacial and terminates at 10.5 ka (Muñoz Sobrino et al., 2013). The reconstructed July temperature during the early Holocene is ca. 12–13 ^∘C, which is considerably cooler than today at this site. However, the authors caution that these samples have poor analogues and the record should be interpreted with caution. Chironomid-based reconstructions of July temperature at Basa de la Mora (Tarrats et al., 2018), a high-elevation site in the Pyrenees, indicate temperatures within ± 0.5 ^∘C of the modern climate during the early to mid-Holocene (10–6 ka), similar to our regional composite reconstructions. However, they show persistently cooler conditions than present by ca. 1.5 ^∘C between 4.5 and 2 ka, not seen in our reconstructions. Furthermore, a direct comparison of our reconstructions of MTWA at Basa de la Mora (Fig. S9) to the chironomid-based reconstructions highlights that the two records show very different trajectories, since the pollen-based reconstruction of this site shows a consistent warming trend throughout the Holocene. Although Tarrats et al. (2018) argue that discrepancies between their temperature reconstructions and pollen-based reconstructions reflect the fact that the vegetation of Iberia, including the mountain areas, is largely driven by moisture changes and perhaps is not a good indicator of temperature, we have shown that there is sufficient information in the pollen records to reconstruct temperature and moisture independently (Tables 3 and S2). Thus, the cause of the differences between the pollen-based and chironomid-based reconstructions at Basa de la Mora is presumably related to methodology. In particular, the chironomid reconstructions use a training data set that does not include samples from the Pyrenees, or indeed the Mediterranean more generally, and may therefore not provide good analogues for Holocene changes at this site.

The lack of a clear trend in MTWA in our reconstructions (Fig. 4b) is not consistent with insolation forcing (Fig. 4e), which shows a declining trend during the Holocene, nor is it consistent with simulated changes in MTWA in transient climate model simulations of the summer temperature response to changing insolation forcing over the Holocene in this region (Fig. 8). The change in moisture gradient during the mid-Holocene, however, suggests an alternative explanation whereby changes in summer temperature are a response to land-surface feedbacks associated with changes in moisture (Fig. 6). Specifically, the observed increased advection of moisture into eastern Iberia would have created wetter conditions there, which in turn would permit increased evapotranspiration, implying less allocation of available net radiation to sensible heating, and resulting in cooler air temperatures. Our reconstructions show that the west–east moisture gradient in mid-Holocene (Fig. 5) was significantly flatter than the steep moisture gradient today (Fig. 2), implying a significant increase in moisture advection into the continental interior during this period. Mauri et al. (2015) also showed that summers were generally wetter than at present in the east but drier than at present in the west in the early to mid-Holocene, supporting the idea of a flatter west–east gradient.

We have shown that stronger moisture advection is not a feature of transient climate model simulations of the Holocene, which may explain why these simulations do not show a strong modification of the insolation-driven changes in summer temperature (Fig. 8). Although the amplitude differs, all of the models show a general decline in summer temperature. The failure of the current generation of climate models to simulate the observed strengthening of moisture transport into Europe and Eurasia during the mid-Holocene has been noted for previous versions of these models (e.g. Bartlein et al., 2017; Mauri et al., 2014) and also shown in Fig. S8. Mauri et al. (2014), for example, showed that climate models participating in the last phase of the Coupled Model Intercomparison Project (CMIP5/PMIP3) were unable to reproduce reconstructed climate patterns over Europe at 6000 years BP and indicated that this resulted from over-sensitivity to changes in insolation forcing and the failure to simulate increased moisture transport into the continent. Bartlein et al. (2017) showed that the CMIP5/PMIP3 models simulated warmer and drier conditions in mid-continental Eurasia at 6000 years BP, inconsistent with palaeoenvironmental reconstructions from the region, as a result of the simulated reduction in the zonal temperature gradient, which resulted in weaker westerly flow and reduced moisture fluxes into the mid-continent. They also pointed out the strong feedback between drier conditions and summer temperatures. The drying of the mid-continent is also a strong feature of the mid-Holocene simulations made with the current generation of CMIP6/PMIP4 models (Brierley et al., 2020). The persistence of these data–model mismatches highlights the need for better modelling of land-surface feedbacks on atmospheric circulation and moisture.

There are comparatively few pollen-based reconstructions of moisture changes during the Holocene from Iberia. Records from Padul show increased mean annual and winter precipitation during the early and mid-Holocene (Camuera et al., 2022; García-Alix et al., 2021). Reconstructions of mean annual and winter precipitation (Camuera et al., 2022) and the ratio of annual precipitation to annual potential evapotranspiration (Wei et al., 2021) also show wetter conditions at this time at El Cañizar de Villarquemado. Both of these sites lie in the eastern part of the Iberian Peninsula, so these reconstructions are consistent with our interpretation of wetter conditions in this region during the interval between 9.5 and 3.5 ka. Ilvonen et al. (2022) provide pollen-based reconstructions of mean annual, summer and winter precipitation from eight sites in Iberia, using WA-PLS and a Bayesian modelling approach. Although they focus on the contrasting pattern of hydroclimate evolution between northern and southern Iberia, the three easternmost sites (San Rafael, Navarres and Qintanar de la Sierra) show much wetter conditions during the early to mid-Holocene. With the exception of the record from Monte Areo, the records from further west are relatively complacent and indeed two sites (Zalamar, El Maillo) show decreased precipitation between 8 and 4 ka. Thus, these records are consistent with our interpretation that the west–east gradient of moisture was reduced between 9.5 and 3.5 ka (Fig. 5).

Speleothem oxygen-isotope data from the Iberian Peninsula provide support for our pollen-based reconstructions of changes in the west–east gradient of moisture through the Holocene. The speleothem records show a progressive increase in temperature from the Younger Dryas onwards, although the trend is less marked in the west than in the east (Baldini et al., 2019). This warming trend is consistent with our reconstructions of changes in MTCO through the Holocene. Speleothem records also show distinctly different patterns in moisture availability, with sites in western Iberia indicating wetter environments during the early Holocene and a transition to drier conditions from ca. 7.5 cal ka BP to the present (Stoll et al., 2013; Thatcher et al., 2020), while eastern sites record wetter conditions persisting from 9 to 4 cal ka (Walczak et al., 2015). This finding would support the weaker west to east moisture gradient shown by our results.

Pollen data are widely used for the quantitative reconstruction of past climates (see discussion in Bartlein et al., 2011), but reconstructions of moisture indices are also affected by changes in water-use efficiency caused by the impact of changing atmospheric CO₂ levels on plant physiology (Farquhar, 1997; Gerhart and Ward, 2010; Prentice et al., 2017; Prentice and Harrison, 2009). This has been shown to be important on glacial–interglacial timescales, when intervals of lower-than-present CO₂ result in vegetation appearing to reflect drier conditions than were experienced in reality (Prentice et al., 2011, 2017; Wei et al., 2021). We do not account for this CO₂ effect in our reconstructions of α because the change in CO₂ over the Holocene was only 40 ppm. This change, relative to modern levels, only has a small impact on the reconstructions (Prentice et al., 2022) and is sufficiently small to be within the reconstruction uncertainties. Furthermore, accounting for changes in CO₂ would not affect the reconstructed west–east gradient through time.

A more serious issue for our reconstructions may be the extent to which the vegetation cover of Iberia was substantially modified by human activities during the Holocene. Archaeological evidence shows that the introduction of agriculture during the Neolithic transition occurred ca. 7.6 ka in some southern and eastern areas of the Iberian Peninsula but spread slowly, and farming first occurred only around 6 ka in the northwest (Drake et al., 2017; Fyfe et al., 2019; Zapata et al., 2004). Anthropogenic changes in land use have been detected at a number of sites, based on pollen evidence of increases in weeds, the presence of cereals (e.g. Abel-Schaad and López-Sáez, 2013; Cortés Sánchez et al., 2012; López-Merino et al., 2010; Mighall et al., 2006; Peña-Chocarro et al., 2005) or the presence of fungal spores associated with animal faeces which has been used to identify the presence of domesticated animals (e.g. López-Sáez and López-Merino, 2007; Revelles et al., 2018). The presence of cereals is the most reliable source of data on human activities, but most cereals only release pollen during threshing and thus are not found in abundance in pollen diagrams from natural (as opposed to archaeological) sites (Trondman et al., 2015). Indeed, it is only after ca. 1 ka that the number of sites which record cereal pollen exceeds the number of sites at which cereals are not represented (Githumbi et al., 2022). Thus, while anthropogenic activities may have been important at the local scale and particularly in the later Holocene (e.g. Connor et al., 2019; Fyfe et al., 2019; Githumbi et al., 2022), most of the sites used for our reconstructions are not associated with archaeological evidence of agriculture or substantial landscape modification. Furthermore, the consistency of the reconstructed changes in climate across sites provides support for these being largely a reflection of regional climate changes rather than human activities.

We have used a modified version of fxTWA-PLS to reconstruct Holocene climates of the Iberian Peninsula because this modification reduced the compression bias in MTCO and MTWA, and specifically reduces the maximum bias in MTCO, MTWA and α. Although this modified approach produces better overall reconstructions (Appendix A), its use does not change the reconstructed trends in these variables through time (Fig. S10). Thus, the finding that winter temperatures are a direct reflection of insolation forcing while summer temperatures are influenced by land-surface feedbacks and changes in atmospheric circulation is robust to the version of fxTWA-PLS used. However, while we use a much larger data set than previous reconstructions, the distribution of pollen sites is uneven and the northern part of the Peninsula is better sampled than the southwest, which could lead to some uncertainties in the interpretation of changes in the west–east gradient of moisture. It would therefore be useful to specifically target the southwestern part of the Iberian Peninsula for new data collection. Alternatively, it would be useful to apply the approach used here to the whole of Eurasia, given that the failure of state-of-the-art climate models to advect moisture into the continental interior appears to be a feature of the whole region (Bartlein et al., 2017) and not the peninsula alone.

5 Conclusion

We have developed an improved version of fxTWA-PLS which further reduces compression bias and provides robust climate reconstructions. We have used this technique with a large pollen data set representing 117 sites across the Iberian Peninsula to make quantitative reconstructions of summer and winter temperature and an index of plant-available moisture through the Holocene. We show that there was a gradual increase in winter temperature through the Holocene and that this trend broadly follows the changes in orbital forcing. Summer temperatures, however, do not follow the changes in orbital forcing but appear to be influenced by land-surface feedbacks associated with changes in moisture. We show that the west–east gradient in moisture was considerably less pronounced during the mid-Holocene, implying a significant increase in moisture advection into the continental interior resulting from changes in circulation. Our reconstructions of temperature changes are broadly consistent with previous reconstructions but are more solidly based because of the increased site coverage. Our reconstructions of changes in the west–east gradient of moisture during the early part of the Holocene are also consistent with previous reconstructions, although this change is not simulated by state-of-the-art climate models, implying that there are still issues with resolving the associated land-surface feedbacks in these models. Our work provides an improved foundation for documenting and understanding the Holocene palaeoclimates of Iberia.

Appendix A: Theoretical basis

The previous version of fxTWA-PLS (fxTWA-PLS1)

The estimated optimum (${\widehat{u}}_k$) and unbiased tolerance (${\widehat{t}}_k$) of each taxon are calculated from the modern training data set as follows:

$$\begin{array}{}\text{(A1)}& {\displaystyle }{\widehat{u}}_k={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n{y}_{ik}{x}_i}{{\sum }_{i=\mathrm{1}}^n{y}_{ik}}},\text{(A2)}& {\displaystyle }{\widehat{t}}_k=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n{y}_{ik}{\left(x_i-{\widehat{u}}_k\right)}^{\mathrm{2}}}{(\mathrm{1}-\mathrm{1}/N_{\mathrm{2}k}){\sum }_{i=\mathrm{1}}^n{y}_{ik}}}},\end{array}$$

where

$$\begin{array}{}\text{(A3)}& N_{\mathrm{2}k}={\displaystyle \frac{\mathrm{1}}{{\sum }_{i=\mathrm{1}}^n{\left(\frac{y_{ik}}{{\sum }_{i^{\prime }=\mathrm{1}}^n{y}_{i^{\prime }k}}\right)}^{\mathrm{2}}}},\end{array}$$

where n is the total number of sites; y_ik is the observed abundance of the kth taxon at the ith site; x_i is the observed climate value at the ith site and N_2k is the effective number of occurrences for the kth taxon.

The fx correction is applied as weight in the form of $\mathrm{1}/{\text{fx}}^{\mathrm{2}}$ at regression at step 7 in Table 1 in Liu et al. (2020). The regression step uses robust linear model fitting by the R code:

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}\text{rlm}(x_i\sim {\text{comp}}_{\mathrm{1}}& +{\text{comp}}_{\mathrm{2}}+\mathrm{\dots }\\ & +{\text{comp}}_{\text{pls}},\text{weights}=\mathrm{1}/{\text{fx}}^{\mathrm{2}}).\end{array}\end{array}$$

The modified version of fxTWA-PLS (fxTWA-PLS2)

The distribution of y_ik is influenced by the distribution of the climate variable, so we need to apply the fx correction when calculating optimum and tolerance for each taxon as follows:

$$\begin{array}{}\text{(A5)}& {\displaystyle }{\widehat{u}}_k={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n\frac{y_{ik}{x}_i}{f_{x_i}}}{{\sum }_{i=\mathrm{1}}^n\frac{y_{ik}}{f_{x_i}}}},\text{(A6)}& {\displaystyle }{\widehat{t}}_k=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n\frac{y_{ik}{\left(x_i-{\widehat{u}}_k\right)}^{\mathrm{2}}}{f_{x_i}}}{\left(\mathrm{1}-\frac{\mathrm{1}}{N_{\mathrm{2}k}}\right){\sum }_{i=\mathrm{1}}^n\frac{y_{ik}}{f_{x_i}}}}},\end{array}$$

where

$$\begin{array}{}\text{(A7)}& N_{\mathrm{2}k}={\displaystyle \frac{\mathrm{1}}{{\sum }_{i=\mathrm{1}}^n{\left(\frac{\frac{y_{ik}}{f_{x_i}}}{{\sum }_{i^{\prime }=\mathrm{1}}^n\frac{y_{i^{\prime }k}}{f_{x_{i^{\prime }}}}}\right)}^{\mathrm{2}}}}.\end{array}$$

The modified version of fxTWA-PLS applies fx correction separately at taxon calculation and regression (step 2 and step 7 in Table 1 in Liu et al., 2020), both using weight in the form of $\mathrm{1}/\text{fx}$. The regression step (step 7) then becomes

$$\begin{array}{}\text{(A8)}& \begin{array}{rl}\text{rlm}(x_i\sim {\text{comp}}_{\mathrm{1}}& +{\text{comp}}_{\mathrm{2}}+\mathrm{\dots }\\ & +{\text{comp}}_{pls},\text{weights}=\mathrm{1}/\text{fx}).\end{array}\end{array}$$

The previous version uses fx values extracted from histograms, and different bin widths may result in different training results. The modified version applies P-splines histogram smoothing (Eilers and Marx, 2021) with third-order difference penalty, which makes the fx values almost independent of the bin width. The optimal smoothing parameter of the P-splines penalty was determined by the HFS (Harville–Fellner–Schall) algorithm (Eilers and Marx, 2021) for the Poisson likelihood for the histogram counts.

Table A1 Leave-out cross-validation (with geographically and climatically close sites removed) fitness of the previous and modified version of fxTWA-PLS (fxTWA-PLS1 and fxTWA-PLS2, respectively), for the mean temperature of the coldest month (MTCO), mean temperature of the warmest month (MTWA) and plant-available moisture (α), using bins of 0.02, 0.02 and 0.002, respectively; n is the number of components used; Avg.bias is the average bias, while Max.bias is the maximum absolute bias and Min.bias is the minimum absolute bias; RMSEP is the root mean square error of prediction; and ΔRMSEP is the per cent change rate of RMSEP, which is (RMSEP_n − RMSEP_n−1)$/$RMSEP_n−1 converted into percentage; when n=1, RMSEP₀ is the RMSEP of null model. p assesses whether using the current number of components is significantly different from using one component less, which is used to choose the last significant number of components (indicated in bold) to avoid overfitting. The degree of overall compression is assessed by doing linear regression to the cross-validation result and the climate variable; b₁ and b₁.se are the slope and the standard error of the slope, respectively. The closer the slope (b₁) is to 1, the lower the overall compression is. The fx correction is set intrinsically in functions of the fxTWA-PLS package for both versions in this paper, instead of relying on an outside input in Liu et al. (2020), so the values of fxTWA-PLS1 might be slightly different from values in Table 3 of Liu et al. (2020), but it does not affect the conclusion.

Figure A1 Training result using the last significant number of components for mean temperature of the coldest month (MTCO), mean temperature of the warmest month (MTWA) and plant-available moisture (α). The left column (a, c and e) shows the previous version (fxTWA-PLS1) and the right column (b, d and f) shows the modified version of fxTWA-PLS (fxTWA-PLS2). The 1:1 line is shown in black; the linear regression line is shown in red to show the degree of overall compression. The horizontal dashed lines indicate the natural limit of α (0–1.26).

Figure A2 Residual using the last significant number of components for mean temperature of the coldest month (MTCO), mean temperature of the warmest month (MTWA) and plant-available moisture (α). The left column (a, c and e) shows the previous version (fxTWA-PLS1) and the right column (b, d and f) shows the modified version (fxTWA-PLS2) of fxTWA-PLS. The zero line is shown in black; the locally estimated scatterplot smoothing is shown in red to show the degree of local compression.