Spatial climate dynamics in the Iberian Peninsula since 15 000 Yr BP

Abstract. Climate changes in the Iberian Peninsula since the Last Glacial Maximum are associated with distributional shifts of major Mediterranean and European temperate species. The dynamic relationship between climate and species in the past may be retrieved from the fossil records available in the Iberian Peninsula. We have used an extensive set of pollen records to reconstruct spatial layers (1 kyr interval) of January minimum temperature, July maximum temperature, and annual precipitation over the time period between 15 and 3 ka. A functional principal component analysis was used to summarise the spatial evolution of climate in areas that share similar climate trends. When compared between them, the identified four areas show different climate trends over the studied period and are coherent with the existence of multiple refugial areas within the Iberian Peninsula.


Introduction
The distribution pattern of biodiversity today is the result of a dynamic process driven by geological events and climatic oscillations at a broad temporal scale (Hewitt, 2000).The climate change since the last glacial period was tracked by species through major range shifts, migrations, and/or ex-tinctions which may be analysed at the genetic level or from the fossil record (Hewitt, 2000;Taberlet and Cheddadi, 2002;Cheddadi et al., 2014).The relationship between climate and biodiversity will be maintained in the future with major consequences due to the current trend of climate warming related to anthropogenic activities, including range shifts (Parmesan and Yohe, 2003;Rebelo et al., 2010), diversity depletion (Araújo et al., 2006;Sinervo et al., 2010), and, more dramatically, species extinction (Thomas et al., 2004).The biodiversity hotspots retain high levels of endemism and are considered as the best candidates for preserving species diversity for the future (Myers et al., 2000).The Mediterranean Basin hotspot, in particular, played the role of refugium to diverse ecosystems over several hundreds of millennia (Wijmstra, 1969;Wijmstra and Smith, 1976;Van der Wiel and Wijmstra, 1987a, b;Tzedakis et al., 2002).Often, those areas where species have persisted during glacial times are referred to as glacial refugia (Bennett and Provan, 2008;Carrión et al., 2010b;Hewitt, 2000;Hu et al., 2009;MacDonald et al., 2008;Willis et al., 2010), and the predicted high levels of diversity found at species level in these areas are corroborated at molecular level (Hewitt, 2000;Petit et al., 2003).Understanding how the past processes impacted biodiversity patterns offers invaluable knowledge for the current species Published by Copernicus Publications on behalf of the European Geosciences Union.P. Tarroso et al.: Spatial climate dynamics in the Iberian Peninsula conservation effort dealing with the ongoing global climate change (Anderson et al., 2006;Willis et al., 2010).
Species' glacial refugia have been generally defined based on species survival with a strong relationship with climate (Hewitt, 2000;Bennett and Provan, 2008;Cheddadi and Bar-Hen, 2009;Médail and Diadema, 2009).Nevertheless, the term has been used recently with multiple definitions (Bennett and Provan, 2008;Ashcroft, 2010).The classic definition of refugia is related to the physiological limits of species that under an increasingly stressing environment experience distributional shifts to near suitable areas (Bennett and Provan, 2008).Palaeoenvironmental and molecular data have proven useful to locate species diversity and migration routes (Petit et al., 2003;Cheddadi et al., 2006Cheddadi et al., , 2014)).However, the locations and extension range of putative refugia still lack spatial consensus and quantification of the refugia's dynamic nature.Reconstructing past environments from proxy data will help understand climate dynamics and how it may have affected biodiversity patterns.The past climate changes, species distributions, and the interplay between them may be reconstructed from the fossil record.Fossil pollen records have proven to be an appropriate proxy for quantifying past climate variables (Webb et al., 1993;Cheddadi et al., 1997;Guiot, 1997;Davis et al., 2003;Bartlein et al., 2010).Using proxy data to derive a definition of refugia in terms of suitable climate in a spatial context may provide further insights into not only the persistence of species in the past but also the location of potential areas that may serve as future refugia for the species persistence.
Climate oscillations in Europe during the last 15 000 years exhibited latitudinal and longitudinal variations (Cheddadi et al., 1997;Davis et al., 2003;Roucoux et al., 2005;Cheddadi and Bar-Hen, 2009;Carrión et al., 2010b).During the Last Glacial Maximum (LGM), several species persisted in refugia located in the southern peninsulas (Hewitt, 2000;Tzedakis et al., 2002;Petit et al., 2003;Weiss and Ferrand, 2007;Bennett and Provan, 2008;Hu et al., 2009;Médail and Diadema, 2009;Ohlemüller et al., 2012).The Iberian Peninsula, with a milder climate than northern European latitudes (Renssen and Isarin, 2001;Carrión et al., 2010b;Perez-Obiol et al., 2011), served as a general refugium to several species that persisted in this area during the LGM.The current patterns of high biological diversity in the Iberian Peninsula derive partially from this favourable climate during harsh glacial conditions and highlight the importance of this area in the broader Mediterranean hotspot (Médail and Quézel, 1999;Cox et al., 2006).However, the Iberian Peninsula is not a geographically homogenous area.Currently Iberian Peninsula is divided into two main climate zones: the temperate at the northern portion of the peninsula and the mediterranean, occupying most of the central and southern part (Olson et al., 2001).This pronounced difference in climate patterns in the Iberian Peninsula also promotes differentiation of the biodiversity patterns (Sillero et al., 2009).Additionally, the past vegetation and climate dynamics reveal a quite com- plex picture (Roucoux et al., 2005;Naughton et al., 2007;Perez-Obiol et al., 2011).Thus, multiple areas were identified as small refugia which lead to the refugia-within-refugia "concept" (Weiss and Ferrand, 2007).Altogether it renders the Iberian Peninsula as a unique area for studying the late-Quaternary climate processes with a highly dynamic vegetation response (Carrión et al., 2010b) which is of a high importance for biodiversity conservation.
Our main objective in this study is to define areas within the Iberian Peninsula (Balearic Islands included) that share similar climate trends and which may have served as a potential refugium.We reconstructed three climate variables and quantified their changes over the past 15 000 years.We also summarised the geographical areas that have undergone similar climate changes and analysed their spatial dynamics between 15 000 and 3000 years.

Methods
The study area extends throughout the Iberian Peninsula and the Balearic Islands (Fig. 1).The method used to reconstruct past climate variables is based on the probability density functions (PDFs) of plant taxa identified in fossil pollen records, and it requires a georeferenced distribution of modern plant taxa and a database of modern climate variables.PDFs for each taxon were built relating the modern distributions in the climate space geographically.The raw fossil pollen data were gathered from author's contributions and from the European Pollen Database (www.europeanpollendatabase.net).Each selected site fits quality criteria regarding the number of radiometric dates (> 3 in each site) and a sampling resolution of at least 200 years.
Table 1.Origin and description of the data sources of fossil pollen used to reconstruct the climate in the Iberian Peninsula.Source is either the European Pollen Database (EPD) or author contribution.Longitude and latitudes correspond to the centroid of the nearest cell to the site and altitude as extracted from WorldClim data set, all at 5 spatial resolution.Each site has information about the number of 14 C dates available, the temporal range covered (see also Fig. 1 for the spatial distribution) and the respective biome following the classification of Olson et al. (2001) Using these criteria we selected a total of 31 records which cover different time spans between 15 000 and 3000 yr BP (Table 1; Fig. 1).Although having the LGM as a lower limit would have provided interesting data, the availability of sites for the spatial interpolation is very limited before the Holocene.Thus, we focused on the 15 000 years when there are still 10 sites available that are necessary for a reliable spatial data interpolation (Table 1; Fig. 1).For the climate reconstruction we assume that modern distributions are in equilibrium with climate over the species range.This is a reasonable assumption when considering the spatial resolution of this study.Using georeferenced full plants distributions reduces the bias resulting from local or isolated presence of species.These biases are also balanced by the inclusion of multiple taxa identified in each core for the climate reconstruction.

Data sources
The current distribution data for 246 taxa were obtained by georeferencing the Atlas of Florae Europaeae (Jalas and Suominen, 1972, 1973, 1976, 1979, 1980, 1983, 1986, 1989, 1991, 1994;Jalas et al., 1996Jalas et al., , 1999;;Laurent et al., 2004).We gathered additional occurrence data for the Mediterranean flora from the Global Biodiversity Information Facility data portal (http://www.gbif.org/).These data were checked and corrected by removing species' presences from botanical and herbaria collections and/or observations with lower spatial resolution than 30 (∼ 55 km).The final taxa list and the assignment of pollen taxa to and modern taxa distributions is given in Table S1.
The georeferenced geographical distributions were rescaled to the resolution of 30 (∼ 55 km).The global observed climate data  for January minimum temperature (T jan ), July maximum temperature (T jul ) and annual precipitation (P ann ) data were obtained from WorldClim database (Hijmans et al., 2005, www.worldclim.org) with 5 resolution (∼ 10 km).The climate data were downscaled to the same spatial resolution of the plant distribution data by aggregating the mean value to the resolution of 30 .All computing was performed using R (R Development Core Team, 2012) with the package rgdal (Keitt et al., 2012).

Reconstruction of past climate variables
The climate reconstruction method is based on the PDF of each taxon identified in a fossil dated pollen assemblage.Pollen taxa were assigned to georeferenced plant taxa (see Table S1).This approach was successfully used to reconstruct climate variables from fossil pollen data (Kühl et al., 2002;Cheddadi and Bar-Hen, 2009;Kühl and Gobet, 2010).Using the PDFs intersection of all taxa identified in a fossil sample we obtain the most likely climate value within which the fossil plant assemblage may occur (Kühl et al., 2002).It has been observed that normal and log-normal (right skewed) distributions fitted to temperature and precipitation, respectively, tend to better represent the data (Chevalier et al., 2014).To avoid sampling the climate spatial distribution instead of the species tolerance, we corrected for the potential bias by using binned climate within the species range as a weighting factor for each climate value (Kühl et al., 2002).The chosen bin size is 2 • C for temperature variables and 20 mm for precipitation.This procedure decreases the weight of the most frequent climate values and increases those, in the distribution of the species, that occur less frequently in the study area (Kühl et al., 2002).
The reconstructed climate using the PDF method results from combining the individual PDFs of the species identified in each pollen sample.The product of the PDFs provides the most likely climate value (Kühl et al., 2002;Chevalier et al., 2014).To identify a taxon as present in the sample, a threshold of three pollen grains was chosen.A minimum of five taxa present is required to reconstruct a climate value for each fossil sample.
Using presence data is seen as an advantage of the PDF method (Kühl et al., 2002) as well as a weakness due to the exclusion of the quantitative data resulting from the pollen abundances (Birks et al., 2010).Fluctuations in pollen abundances are related to multiple factors such as the species ecophysiology, differential pollen production, dispersal capacity, and other traits (Hicks, 2006).We have used the pollen proportions to weight the PDFs of the respective taxa.The minimum positive pollen proportion corresponds to the presence of the taxon, while the maximum defines its highest abundance within the fossil record.Using pollen proportions of a taxon within a time series instead of within a sample avoids the bias of differential pollen production and thus allows estimating the presence of a species relative to its maximum percentages in the whole record.
The pollen proportions were converted to alpha values, reducing the species climate tolerance towards the peak density values (Fig. 2).We assumed that the pollen proportion has an inverse relation to the proximity of near-optimal conditions.To avoid the selection of a unique climate value from the PDF when the maximum detection of a species occurs Variable Density 82 % 55 % 10 % Figure 2. Example of the influence of pollen proportion (pp) on the calculation of the density of taxa presence intersection.The shades of grey indicate the effect of different pp when the pollen adjustment value (pa) is set to 0.9 and arrows indicate the assumed presence range.The first case (dark grey) results from pp = 1.0, which represents the highest detectability and is assumed to be found near the core distribution area and, thus, near-optimum conditions.The presence is assumed in a narrow range around peak density with α = pp×pa 2 (corresponding to 10 % of the area).When pp = 0.5 (medium grey) the corresponding area is 55 % and with pp = 0.2 (light grey) is used the widest presence range (82 % of the PDF area).
-i.e. when its pollen proportion is found to be 1, we use a pollen adjustment value set to 0.9.This means that, at the maximum taxon detection, the PDF will be reduced to the area of the density corresponding to 10 % of the probability (Fig. 2).The maximum detection of a taxon indicates a nearoptimal climate niche and the adjustment value set to a value near but not equal to one allows some degree of uncertainty in the reconstruction.On the other hand, setting this value to zero will not allow any influence of the pollen proportion, resulting in a binary presence/absence reconstruction (Kühl et al., 2002).For each sample, the collection of the taxa tolerance intervals built this way are added resulting in a taxon profile, showing where in the climate space the frequency of the taxon is higher, taking into account the proximity to optimal conditions.The final climate reconstruction value is the product of the climate PDF with the taxa profile.The reconstructed value and associated uncertainty are usually extracted from the PDF as the mean and standard deviation (Kühl et al., 2002;Kühl and Gobet, 2010).Assuming a normal distribution, we extract the peak density value and the 95 % confidence interval from the density profile.The confidence interval range shows the uncertainty around the recon- structed value and is related to the standard deviation.The reconstructed values for each site were fitted with a smoothing spline to produce continuous time series, from which 1000year time slices were extracted.
In order to evaluate the robustness of the reconstruction method, we have compared modern reconstructed and observed climate data  from WorldClim database.For reconstructing climate from pollen data, we have used all samples available within the last 500 years.Climate values were averaged for all sites with more than one sample.The correlation between the two data sets was tested using a Pearson's correlation score.To provide the significance of the correlation value, a set of 999 replicates were performed where the observed climate variable was shuffled without repetition.Although this evaluation does not take into account neither the climate oscillations during the last 500 years nor the human disturbances, it still provides a broad evaluation of the reconstruction method because (1) it depicts per site the relationship between observed climate data with reconstructed values and (2) the slope direction of the regression and the related correlation signal indicate that the reconstruction is spatially coherent.A linear regression was used to estimate a baseline for calculating the anomalies at each site using the observed climate.The pre-industrial period around AD 1850 is commonly used as reference climatology to compute anomalies.This period is also often used as a baseline in climate models, facilitating data-model comparisons, and it is less biased with recent climate warming allowing past warming to be better depicted (Davis et al., 2003;Mauri et al., 2015).Although a specific year is selected, the time frame often includes ±500 years (Mauri et al., 2015), which is equivalent to the period we have used in our study.The regression allows a climate baseline to be built without artificially adding samples to compensate for differential number of samples available for recent periods in a 4-D (spatial plus time) interpolation (Mauri et al., 2015) and the linear equations provide all the information to generate the baseline with the observed climate data.

Spatial analysis of past climate
Thirteen climate grids, ranging from 15 000 to 3 000 calendar years BP (hereafter, "ka") with a 1000-year interval, were obtained for each reconstructed variable by spatial interpolation of the climate anomalies at each available site.The anomalies were computed for each site between the reconstructed climate and the modern reference climate.Anomalies were projected onto a 30 (∼ 55 km) resolution grid and interpolated onto a 5 (∼ 10 km) resolution grid using 3-D thin-plate smoothing splines with two spatial dimensions including altitude.This interpolation method generates accurate climate predictions (Jarvis and Stuart, 2001) and it was used for the WorldClim variables (Hijmans et al., 2005).
To further summarise the spatial and temporal variability of the data we applied a functional principal component anal-P.Tarroso et al.: Spatial climate dynamics in the Iberian Peninsula ysis (fPCA).The fPCA extends the exploratory data analysis of the principal component analysis to functional data (Bickel et al., 2005), depicting both spatial and time patterns that are then summarised in a few components.Cheddadi and Bar-Hen (2009) applied a fPCA to nearly the same timescale as the present study to depict January temperature patterns from European pollen data.Here we have broadened the approach to each climate time series available in each grid cell to produce gridded spatial components.The functional data were built by combining B-spline basis functions to fit the time series.We have retained the components that explain more than 90 % of the variance and rescaled the range from −1 to 1.We used hierarchical cluster analysis over the produced first components grids of each variable to identify areas in the Iberian Peninsula that share similar climate trends over the past 15 ka.Climate stability was computed for each variable as the mean absolute deviance from the current climate as available in WorldClim data set.
All analysis were performed using the R Project for Statistical Computing (R Development Core Team, 2012) with packages fields (Furrer et al., 2012), rgdal (Keitt et al., 2012), gstat (Pebesma, 2004), and fda (Ramsay et al., 2012).The climate reconstructions were performed with R scripts developed by the authors and available at request.

Results
The modern climate reconstructions (500 years) show a high degree of agreement with the observed climate data (RMSE Tjan = 5.01; RMSE Tjul = 3.85; RMSE Pann = 399.85;Fig. S1).These data sets show a positive linear trend and a significant positive correlation (p ≤ 0.006 for all variables), revealing that the reconstruction method predicts well the spatial distribution of climate.The standard error associated with the climate reconstruction is on average low but increases with age (Fig. S2).
The reconstructed three climate variables exhibit high spatial variability between 15 and 3 ka (Fig. 3, Fig. S3).The uncertainty associated with the spatial interpolations is usually low, suggesting a good sampling coverage, except in the northwestern area (Fig. S4).The Iberian Peninsula had extensive areas with extremely low T jan that gradually increased markedly after 10 ka.The pattern of T jul over the same time remains stable, with lower values before 12 ka.There is a decreasing trend of precipitation, especially after 10 ka (Fig. 4), which is marked mostly in the south-eastern part of its area (Fig. S3).The clustering of the first fPCA component for the three reconstructed variables are spatially structured (Fig. 5), and allow their overall trends to be summarised (Fig. 4).The first component of each variable explains more than 95 % of the variation (T jan : 95.5 %; T jul : 99.2 %; P ann : 99.5 %).Cluster C1 (27 % of the total area) is located mostly in northern and western Iberia and includes part of the north-Iberian mountain ranges but also low altitudinal coastal areas (average altitude is 679 ± 454 m).This is the wettest cluster with P ann ranging from 1054 to 1115mm, the coldest in July (21.7 < T jul < 24.2 • C) and with very low January minimum temperatures (−5.5 < T jan < 0.2 • C).The C2 cluster encompasses part of the Cantabrian mountain range and the central Iberian system (28 % of the total area with an average altitude of 859 ± 303 m).It occupies most of the northern plateau, where it has the lowest January temperature (−5.7 < T jan < −1.3 • C), whereas July (25.1 < T jul < 27.7 • C) is warmer than within C1.This shows high seasonal amplitude with low precipitation (537 < P ann < 621 mm), similar to C3 and C4.The dissimilarities between clusters C3 and C4 (25 and 20 % of the total area and average altitude of 613 ± 96 and 278 ± 232 m, respectively) concern mainly the temperature.Cluster C4 is warmer and wetter than C3.These are the warmest areas for both January (C3: −1.7 < T jan < 3.0; C4: 1.0 < T jan < 6.4 • C) and July (C3: 29.4 < Tjul < 33.4; C4: 27.2 < T jul < 30.3 • C) and with low annual precipitation (C3: 504 < P ann < 614; C4: 555 < P ann > 682 mm).The Balearic Islands are fully included in the C4 cluster (Fig. 5).
The mean absolute deviance from the current climate shows that the climate stability during the last 15 kyr was not spatially uniform (Fig. 6).T jan and P ann exhibited higher stability in the southern Iberia, although T jan has lower values of deviance (higher stability) towards the eastern coast and P ann towards the western coast.T jul exhibited lower de-viance at higher altitudes, particularly at the central system, northern mountains, and Pyrenees, but also in the southern Sierra Morena.

Discussion
Fossil pollen data provide a record of vegetation changes which constitutes a valuable proxy for reconstructing past climate changes, especially using large data sets (Bartlein et al., 2010).The method used here provides reliable climate reconstructions, despite the low number of sequences selected according to our quality criteria for spatial climate reconstruction, both in terms of sampling resolution and number of 14 C dates.The western part of the peninsula has a better data coverage which provides more robust spatial interpolations, particularly for the most recent to middle time periods analysed.Nevertheless, the spatial uncertainty related to the interpolation shows a uniform variance for all time periods (Fig. S4).The only exception is the north-western part of the study area, where the lack of data promotes higher uncertainty for the spatial interpolation.The residuals between observed climate and reconstructed climate were high, resulting also in a low coefficient of determination for the linear regression (Fig. S1).However, this is expectable since we were comparing observed climate data with reconstructed values of the last 500 years and averaging the climate variation in this period tend to increase the residuals.In addition, the anthropogenic impact on the ecosystems is likely also biasing the results.Nevertheless, a positive linear trend with a significant positive correlation was found between reconstructed climate and observed climate that allows us to produce a reference data set using this model and the observed climate.The results provided here reinforce the role of the Iberian Peninsula as a glacial refugium and holding enough climate variation (Fig. 5) to support a network of smaller refugial areas (Weiss and Ferrand, 2007).
Climate of the last 15 kyr was dynamic, with oscillations of temperature and precipitation occurring mostly at the southern part of the peninsula.Given the link between climate and species distributions (Hewitt, 2000), it is likely that these changes had an impact on the location, extent and evolution of the refugia and the recolonisation processes during the post-glacial period.Nonetheless, the reconstructed overall trend is a noticeable warming in winter temperatures after 15 ka, particularly between 12 and 9 ka (Fig. 4) that is likely due to the increase in the summer insolation in the northern hemisphere (Berger, 1978).This warming trend tends to reduce the area in the Iberian Peninsula with very low temperatures (Fig. 3).Although insolation peaks at 9 ka and decreases afterwards, it does not translate to a general cooling and in south-western Europe is seen an increase in insolation in both summer and winter (Davis et al., 2003).A striking pattern is the partitioning of the peninsula in spatially structured areas that shared similar climate trend over the late Quaternary (Fig. 5).The wettest and cold cluster C1 is predominantly located at the northern and north-western Iberia and occupies most of the current temperate climate zone.Although very similar to C2, it contrasts in the seasonal amplitude and precipitation amount.Interestingly, the pattern of current bioclimate zones in Iberian Peninsula is retrieved on the clusters scheme, suggesting the persistence of a transition area between very different climate zones, although the magnitude of the differences have changed in the past.
Our results show that January temperatures exhibited a general warming trend over the last 15 000 years which corresponds on average to an increase of ∼ 5.5 • C. The southern part of the peninsula is more resilient to change, particularly for T jan and P ann , whereas the northern part recorded major changes.This pattern is less obvious for July temperatures, where variations showed a smaller amplitude even though this variable is markedly different between clusters, thus contributing to the climate split of the study area (Fig. 4, Fig. S5).The minimum winter temperatures constrain the physiologic ability of plants to further development and, thus, are a major factor restricting distributions (Sykes et al., 1996).Higher summer insolation provides enough energy to plant growth and July temperature in the Mediterranean is a less limiting variable for growth than T jan which makes the reconstruction of summer months and its interpretation more complex.

The end of the Pleistocene
The 1000-year time interval provides enough resolution to analyse general patterns of climate evolution.However, abrupt climate events are not detectable.The end of the Oldest Dryas (OD; ending around 14.5 ka) is characterised in Iberia by a vegetation change compatible with cold and humid conditions (Naughton et al., 2007(Naughton et al., , 2015) ) and is followed by the Bölling-Allerød warm period (B-A; ending around 13 ka).Our results show a similar pattern with colder conditions between 15 ka and 13 ka and a higher humidity (Fig. 4), particularly evident in the central and southern clusters.Although all clusters show a similar trend, the C1 and C2 are colder than average.The general pattern in the Iberian Penin-sula is a contrast between a colder north and a warmer south but, nevertheless, an area dominated by low January temperatures (Fig. 3, Fig. S3, S5).The evolution of precipitation during the last 15 kyr in the Iberian Peninsula shows a very stable pattern: northern areas comprised in C1 had high precipitation values during the period analysed, while the south was wetter than today (Fig. S3, S5).The increase in the moisture availability during the B-A (Naughton et al., 2015) is in line with the slight increase in precipitation in all clusters between 14 and 13 ka (Fig. 3).
As described earlier in Europe (Renssen and Isarin, 2001;Heiri et al., 2004), T jan shows wider changes in amplitude than T jul .The cold to warm transitions that occurred at ∼ 14.7 and 11.5 ka (Renssen and Isarin, 2001;von Grafenstein et al., 2012) in Europe had a spatial impact that is noticeable in the reconstructed temperatures (Fig. 3, Fig. S3, S5).

The Holocene
The B-A warm stage is followed by the cold Younger Dryas (YD; between ∼ 12.9 and ∼ 11.7 ka), marking the beginning of the Holocene.This period records a warming trend for T jan , while T jan decreased abruptly (Fig. 4) with a reduction of the warmer areas between 13 and 12 ka (Fig. 3).
The Holocene warm period (approximately between ∼ 8.2 and 5.6 ka, depending on the location in Europe) is characterised by increasing summer temperatures (Seppä and Birks, 2001).Such trend is more obvious in northern Europe and the Alps, while we rather observe a cooling at lower latitudes (Davis et al., 2003).Our results point to a slight decrease in T jan and T jul around 7 ka, but the overall temperature pattern is rather stable.This is likely affected by the temporal resolution of this study, failing to clearly detect rapid events.P ann shows a slightly wetter climate at 7 ka (Fig. 3) which is consistent with earlier reconstruction for the southern European lowlands (Cheddadi et al., 1997).
Between 6 and 3 ka, areas with low precipitation expand in the Iberian Peninsula (Fig. 3), which allows the expansion of the Mediterranean taxa (Naughton et al., 2007;Carrión et al., 2010b, a).The increasing aridity trend in the south is balanced by the high precipitation values in the north (Fig. 4, Fig. S3, S5), contributing to the shaping of the current Iberia pattern of two contrasting bioclimatic regions: the north is temperate and wet while the south is a dry and warm.
The behaviour of the reconstructed variables at 5 ka is likely to be influenced by non-natural ecosystem changes due to human activities such as the forest degradation that began in lowlands and later in mountainous areas (Carrión et al., 2010a).These human impacts add confounding effects in the fossil pollen record and may lead to slightly biased temperature reconstructions after 5 ka.On the other hand, human impacts at larger scales, capable of leaving noticeable imprints on landscape, were likely to happen later (Carrión et al., 2010a) and, furthermore, there is evidence of a cooling and drier stage in the Iberian Peninsula after 5 ka (Dorado Valiño et al., 2002).

Climate role in Iberian refugia
The climate change since the LGM in the Iberian Peninsula had an impact on the persistence of temperate species, migrating pathways, and on the overall recolonisation processes during the post-glacial period within the peninsula (Hewitt, 2000;Naughton et al., 2007;Carrión et al., 2010b).During this period, climate favoured migrations and expansion processes that culminated in secondary contacts for several lineages previously isolated in patches of suitable habitat (Branco et al., 2002;Godinho et al., 2006;Weiss and Ferrand, 2007;Miraldo et al., 2011).Particularly the B-A warming phase and the warming stage after the YD, which we show here and which have highly affected the spatial organisation of the climate in the Iberian Peninsula, are likely favouring expansion processes of warmth-dependent organisms.Given the relationship between climate change and biodiversity patterns, the clustering scheme (Fig. 5) depicting areas with different climate evolution is consistent with the molecular evidence of a network of putative refugia within Iberia (Weiss and Ferrand, 2007).Refugia have been associated with climate and habitat stability, with both playing complementary roles (Ashcroft, 2010).However, as shown by large-scale landscape analysis (Carrión et al., 2010b, a) and climate reconstructions (Davis et al., 2003;Cheddadi and Bar-Hen, 2009), both have a strong dynamic nature in the Iberian Peninsula and likely promoted the formation of patches of suitable habitat during harsh conditions.The highly structured populations that many species exhibit in the Iberian Peninsula have contributed decisively to the idea of refugia diversity (Hewitt, 2000;Weiss and Ferrand, 2007).Overall, the information included in the multidimensional climate data allowed us to define areas characterised by a climate evolution during the late Quaternary with smaller amplitude of change (clusters C3 and C4).These areas showed higher stability of both T jan and P ann (Fig. 6).Cluster C4 coincides at a great extent with areas that offered more resilience to change between millennia (Fig. 5).Within these areas, temperature and precipitation were suitable to support the survival of temperate trees, likely acting as glacial refugia.On the other hand, the cold areas of the first and second cluster also associated with faster changes cluster likely diminished the suitability for the long term persistence of species.One might infer that the defined clusters are associated with potential isolation or dispersal events of species throughout the studied time span.Particularly, the fourth cluster (Fig. 5) includes areas that have already been described as glacial refugia for several animal and plant species (Weiss and Ferrand, 2007; see chapter 5 for a review of refugia in the Iberian Peninsula).In the area represented by this cluster, the reconstructed T jan indicate a mild climate with higher precipitation than currently, which is compatible with the persistence of species in these areas.The southern plateau, mostly comprised in the second cluster (Fig. 5), also recorded mild conditions which are often associated with southern refugia; however, the rapid T jul oscillations associated with a cold T jan and low precipitation may have prevented long-term persistence but are likely compatible with a recolonisation process.
The pattern of stability indicates a southern Iberia with less change, particularly T jan and P ann .High altitudes offer more resilience to change, particularly to July temperature, and lower areas may be swept rapidly with occurring changes (see Fig. S6).Our data suggest that, at the regional scale and with extensive time-series data, this relation is preserved.Areas of lower velocity of change, which are hence more stable, are associated with high levels of endemicity at global scales (Sandel et al., 2011), and areas of high velocity are often associated with species extinction (Nogués-Bravo et al., 2010).Our results indicate higher stability in the southern part of the Peninsula, similar to other studies based on climate data (Ohlemüller et al., 2012).However, our studied time frame extends to 15 ka, which does not cover the glacial maximum (∼ 21 ka).At that time period, a higher degree of fragmentation of the stability is expected due to colder conditions, and areas compatible with refugia would be also less contiguous.These could be seen as a macrorefugia, offering conditions for large population during glacial times (Mee and Moore, 2014).Microrefugia are known to occur in the northern areas of the Iberian Peninsula (e.g.Fuentes-Utrilla et al., 2014), but the spatial scale used and the number of pollen sites available render microrefugia undetectable in this study.

Conclusions
The reconstruction of past climates using biological data is an invaluable resource for the study of the dynamics of glacial refugial areas.Although there is a limited number of available sites and time range coverage, the spatial combination of fossil pollen data provides a continuous record with a climate signal that can be translated into spatially explicit analysis of climate dynamics.
The reconstructed climate variables for the post-glacial period show different patterns of evolution but are clearly marked by the lasting impact of climatic events.The Iberian Peninsula had areas that shared similar climate evolution during the late Quaternary.Some areas that we have suggested as potential refugia are consistent with those areas where genetic diversity was found to be high and which are often considered as refugial areas for several animal and plant species.
The analysis of these areas and the related climate provides new insights into the dynamics of refugia through time and space, which helps in gaining a better understanding of the evolution of biodiversity hotspots both at the species and the intraspecific levels.Linking past climate and diversity in the Iberian Peninsula is a major issue for conservation issues, especially under the expected future climate change.a PhD grant (SFRH/BD/42480/2007) and a post-doc grant (SFRH/BPD/93473/2013) and José Carlos Brito is supported by a contract (IF/00459/2013), both from Fundação para a Ciência e Tecnologia.José Carrión's contribution was funded by the project Paleoflora y Paleovegetación ibérica, Plan Nacional de I+D+i, Ref. CGL-2009-06988/BOS.Luisa Santos acknowledges the contribution of M. C. Freitas and C. Andrade (University of Lisbon), who provided the cores.The authors would like to acknowledge all contributors of the European Pollen Database and the Global Biodiversity Information Facility for making their data sets publicly available to the scientific community.We are very grateful to Basil Davis for his kind support and comments.We also thank William Fletcher and Maria Sanchez-Goñi for data contribution and comments, and also Penélope González-Sampériz for contributions.We are also grateful to Graciela Gil Romera, Ana Ortega, and an anonymous referee for the extensive reviews that greatly improved the quality of the manuscript.This is ISEM contribution no.2016-076.
Edited by: V. Rath

Figure 1 .
Figure 1.Study area with sample points.The black area inside each circle represents the ages available in each pollen sequence.

Figure 3 .
Figure 3. Distribution of the reconstructed climate variables in the Iberian Peninsula and Balearic Islands in the last 15 kyr.Colours show the proportion of area covered with each class of (a) minimum temperature of January, (b) maximum temperature of July, and (c) annual precipitation.

Figure 4 .
Figure 4. Minimum and maximum temperatures of January and July, respectively, and annual precipitation during the last 15 kyr.The solid line represents the average climate in the study area.The remaining lines are the average of each cluster found -C1: shortdash line; C2: dotted line; C3: dash-dot line; and C4: long-dash line.
Hierarchical cluster analysis of the functional PCA components of T jan , T jul , and P ann in the last 15 kyr found in the study area.The top dendrogram represents the size of the clusters of similar climate evolution and the relations between them.Numbers correspond to each identified cluster.

Figure 6 .
Figure 6.Average differences between millennia for each of the climate variables.Calculation of the differences are computed between a given age and the previous one.Isolines in each map indicate the average value of change. .