Pedro Tarroso
email: ptarroso@cibio.up.pt

José Carrión

Miriam Dorado-Valiño

Paula Queiroz

Luisa Santos

Ana Valdeolmillos-Rodríguez

Paulo Célio Alves

José Carlos Brito

Rachid Cheddadi

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

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 [START_REF] Hewitt | The genetic legacy of the Quaternary ice ages[END_REF]. 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 [START_REF] Hewitt | The genetic legacy of the Quaternary ice ages[END_REF][START_REF] Taberlet | Quaternary refugia and persistence of biodiversity[END_REF][START_REF] Cheddadi | Revisiting tree-migration rates: Abies alba (Mill.), a case study[END_REF]. 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 [START_REF] Parmesan | A globally coherent fingerprint of climate change impacts across natural systems[END_REF][START_REF] Rebelo | Predicted impact of climate change on European bats in relation to their biogeographic patterns[END_REF], diversity depletion [START_REF] Araújo | Climate warming and the decline of amphibians and reptiles in Europe[END_REF][START_REF] Sinervo | Erosion of lizard diversity by climate change and altered thermal niches[END_REF], and, more dramatically, species extinction [START_REF] Thomas | Extinction risk from climate change[END_REF]. The biodiversity hotspots retain high levels of endemism and are considered as the best candidates for preserving species diversity for the future [START_REF] Myers | Biodiversity hotspots for conservation priorities[END_REF]. The Mediterranean Basin hotspot, in particular, played the role of refugium to diverse ecosystems over several hundreds of millennia [START_REF] Wijmstra | Palynology of the first 30 m of a 120 m deep section in northern Greece[END_REF][START_REF] Wijmstra | Palynology of the middle part (30-78 metres) of the 120 m deep section in northern Greece (Macedonia)[END_REF]Van der Wiel and Wijmstra, 1987a, b;[START_REF] Tzedakis | Buffered tree population changes in a quaternary refugium: evolutionary implications[END_REF]. Often, those areas where species have persisted during glacial times are referred to as glacial refugia [START_REF] Bennett | What do we mean by "refugia"?[END_REF]Carrión et al., 2010b;[START_REF] Hewitt | The genetic legacy of the Quaternary ice ages[END_REF][START_REF] Hu | Paleoecology meets genetics: deciphering past vegetational dynamics[END_REF][START_REF] Macdonald | Impacts of climate change on species, populations and communities: palaeobiogeographical insights and frontiers[END_REF][START_REF] Willis | Biodiversity baselines, thresholds and resilience: testing predictions and assumptions using palaeoecological data[END_REF], and the predicted high levels of diversity found at species level in these areas are corroborated at molecular level [START_REF] Hewitt | The genetic legacy of the Quaternary ice ages[END_REF][START_REF] Petit | Glacial refugia: hotspots but not melting pots of genetic diversity[END_REF]. 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. conservation effort dealing with the ongoing global climate change [START_REF] Anderson | Linking palaeoenvironmental data and models to understand the past and to predict the future[END_REF][START_REF] Willis | Biodiversity baselines, thresholds and resilience: testing predictions and assumptions using palaeoecological data[END_REF].

Species' glacial refugia have been generally defined based on species survival with a strong relationship with climate [START_REF] Hewitt | The genetic legacy of the Quaternary ice ages[END_REF][START_REF] Bennett | What do we mean by "refugia"?[END_REF][START_REF] Cheddadi | Spatial gradient of temperature and potential vegetation feedback across Europe during the late Quaternary[END_REF][START_REF] Médail | Glacial refugia influence plant diversity patterns in the Mediterranean Basin[END_REF]. Nevertheless, the term has been used recently with multiple definitions [START_REF] Bennett | What do we mean by "refugia"?[END_REF][START_REF] Ashcroft | Identifying refugia from climate change[END_REF]. 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 [START_REF] Bennett | What do we mean by "refugia"?[END_REF]. Palaeoenvironmental and molecular data have proven useful to locate species diversity and migration routes [START_REF] Petit | Glacial refugia: hotspots but not melting pots of genetic diversity[END_REF][START_REF] Cheddadi | Imprints of glacial refugia in the modern genetic diversity of Pinus sylvestris[END_REF][START_REF] Cheddadi | Revisiting tree-migration rates: Abies alba (Mill.), a case study[END_REF]. 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 [START_REF] Webb | Pollen responsesurface estimates of late-Quaternary changes in the moisture balance of the northeastern United States[END_REF][START_REF] Cheddadi | The climate of Europe 6000 years ago[END_REF][START_REF] Guiot | Palaeoclimatology: Back at the last interglacial[END_REF][START_REF] Davis | The temperature of Europe during the Holocene reconstructed from pollen data[END_REF][START_REF] Bartlein | Pollen-based continental climate reconstruction at 6 and 21 ka: a global synthesis[END_REF]. 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 [START_REF] Cheddadi | The climate of Europe 6000 years ago[END_REF][START_REF] Davis | The temperature of Europe during the Holocene reconstructed from pollen data[END_REF][START_REF] Roucoux | The response of NW Iberian vegetation to North Atlantic climate oscillations during the last 65kyr[END_REF][START_REF] Cheddadi | Spatial gradient of temperature and potential vegetation feedback across Europe during the late Quaternary[END_REF]Carrión et al., 2010b). During the Last Glacial Maximum (LGM), several species persisted in refugia located in the southern peninsulas [START_REF] Hewitt | The genetic legacy of the Quaternary ice ages[END_REF][START_REF] Tzedakis | Buffered tree population changes in a quaternary refugium: evolutionary implications[END_REF][START_REF] Petit | Glacial refugia: hotspots but not melting pots of genetic diversity[END_REF]Weiss and Ferrand, 2007;[START_REF] Bennett | What do we mean by "refugia"?[END_REF][START_REF] Hu | Paleoecology meets genetics: deciphering past vegetational dynamics[END_REF][START_REF] Médail | Glacial refugia influence plant diversity patterns in the Mediterranean Basin[END_REF][START_REF] Ohlemüller | Potential source and sink locations for climate-driven species range shifts in Europe since the Last Glacial Maximum[END_REF]. The Iberian Peninsula, with a milder climate than northern European latitudes [START_REF] Renssen | The two major warming phases of the last deglaciation at 14.7 and 11.5 ka cal BP in Europe: climate reconstructions and AGCM experiments[END_REF]Carrión et al., 2010b;[START_REF] Perez-Obiol | Mid-Holocene vegetation and climatic history of the Iberian Peninsula[END_REF], 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 [START_REF] Médail | Biodiversity hotspots in the Mediterranean Basin: setting global conservation priorities[END_REF][START_REF] Cox | The status and distribution of reptiles and amphibians of the Mediterranean Basin, IUCN, Gland[END_REF]. 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 [START_REF] Olson | Terrestrial Ecoregions of the World: A New Map of Life on Earth[END_REF]. This pronounced difference in climate patterns in the Iberian Peninsula also promotes differentiation of the biodiversity patterns [START_REF] Sillero | Biogeographical patterns derived from remote sensing variables: the amphibians and reptiles of the Iberian Peninsula[END_REF]. Additionally, the past vegetation and climate dynamics reveal a quite com- plex picture [START_REF] Roucoux | The response of NW Iberian vegetation to North Atlantic climate oscillations during the last 65kyr[END_REF][START_REF] Naughton | Present-day and past (last 25 000 years) marine pollen signal off western Iberia[END_REF][START_REF] Perez-Obiol | Mid-Holocene vegetation and climatic history of the Iberian Peninsula[END_REF]. 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 [START_REF] Olson | Terrestrial Ecoregions of the World: A New Map of Life on Earth[END_REF] 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[START_REF]Atlas Florae Europaeae. Distribution of vascular plants in Europe, The Committee for Mapping the Flora of Europe and Societas Biologica Fennica Vanamo[END_REF], 1976[START_REF]Atlas Florae Europaeae. Distribution of vascular plants in Europe, The Committee for Mapping the Flora of Europe and Societas Biologica Fennica Vanamo[END_REF][START_REF]Atlas Florae Europaeae. Distribution of vascular plants in Europe, The Committee for Mapping the Flora of Europe and Societas Biologica Fennica Vanamo[END_REF][START_REF]Atlas Florae Europaeae. Distribution of vascular plants in Europe, The Committee for Mapping the Flora of Europe and Societas Biologica Fennica Vanamo[END_REF][START_REF]Atlas Florae Europaeae. Distribution of vascular plants in Europe, The Committee for Mapping the Flora of Europe and Societas Biologica Fennica Vanamo[END_REF][START_REF] Peñalba Garmendia | Dynamique de végétation tardiglaciaire et holocène du centre-nord de l'Espagne d'après l'analyse pollinique[END_REF][START_REF]Atlas Florae Europaeae. Distribution of vascular plants in Europe, The Committee for Mapping the Flora of Europe and Societas Biologica Fennica Vanamo[END_REF], 1994;Jalas et al., 1996Jalas et al., , 1999;;[START_REF] Laurent | Refining vegetation simulation models: from plant functional types to bioclimatic affinity groups of plants[END_REF]. 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 [START_REF] Kühl | Probability Density Functions as Botanical-Climatological Transfer Functions for Climate Reconstruction[END_REF][START_REF] Cheddadi | Spatial gradient of temperature and potential vegetation feedback across Europe during the late Quaternary[END_REF][START_REF] Kühl | Climatic evolution during the Middle Pleistocene warm period of Bilshausen, Germany, compared to the Holocene[END_REF]. 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 [START_REF] Kühl | Probability Density Functions as Botanical-Climatological Transfer Functions for Climate Reconstruction[END_REF]. It has been observed that normal and log-normal (right skewed) distributions fitted to temperature and precipitation, respectively, tend to better represent the data [START_REF] Chevalier | CREST (Climate REconstruction SofTware): a probability density function (PDF)based quantitative climate reconstruction method[END_REF]. 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 [START_REF] Kühl | Probability Density Functions as Botanical-Climatological Transfer Functions for Climate Reconstruction[END_REF]. 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 [START_REF] Kühl | Probability Density Functions as Botanical-Climatological Transfer Functions for Climate Reconstruction[END_REF].

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 [START_REF] Kühl | Probability Density Functions as Botanical-Climatological Transfer Functions for Climate Reconstruction[END_REF][START_REF] Chevalier | CREST (Climate REconstruction SofTware): a probability density function (PDF)based quantitative climate reconstruction method[END_REF]. 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 [START_REF] Kühl | Probability Density Functions as Botanical-Climatological Transfer Functions for Climate Reconstruction[END_REF] as well as a weakness due to the exclusion of the quantitative data resulting from the pollen abundances [START_REF] Birks | Strengths and weaknesses of quantitative climate reconstructions based on Late-Quaternary biological proxies[END_REF]. Fluctuations in pollen abundances are related to multiple factors such as the species ecophysiology, differential pollen production, dispersal capacity, and other traits [START_REF] Hicks | When no pollen does not mean no trees[END_REF]. 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 [START_REF] Kühl | Probability Density Functions as Botanical-Climatological Transfer Functions for Climate Reconstruction[END_REF]. 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 [START_REF] Kühl | Probability Density Functions as Botanical-Climatological Transfer Functions for Climate Reconstruction[END_REF][START_REF] Kühl | Climatic evolution during the Middle Pleistocene warm period of Bilshausen, Germany, compared to the Holocene[END_REF]. Assuming a normal distribution, we extract the peak density value and the 95 % confidence interval from the density profile. 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 [START_REF] Davis | The temperature of Europe during the Holocene reconstructed from pollen data[END_REF][START_REF] Mauri | The climate of Europe during the Holocene: a gridded pollen-based reconwww[END_REF]. Although a specific year is selected, the time frame often includes ±500 years [START_REF] Mauri | The climate of Europe during the Holocene: a gridded pollen-based reconwww[END_REF], 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 [START_REF] Mauri | The climate of Europe during the Holocene: a gridded pollen-based reconwww[END_REF] 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 [START_REF] Jarvis | A Comparison among Strategies for Interpolating Maximum and Minimum Daily Air Temperatures, Part II: The Interaction between Number of Guiding Variables and the Type of Interpolation Method[END_REF] and it was used for the WorldClim variables [START_REF] Hijmans | Very high resolution interpolated climate surfaces for global land areas[END_REF].

To further summarise the spatial and temporal variability of the data we applied a functional principal component anal- ysis (fPCA). The fPCA extends the exploratory data analysis of the principal component analysis to functional data [START_REF] Bickel | Functional data analysis[END_REF], depicting both spatial and time patterns that are then summarised in a few components. [START_REF] Cheddadi | Spatial gradient of temperature and potential vegetation feedback across Europe during the late Quaternary[END_REF] 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 [START_REF] Pebesma | Multivariable geostatistics in S: the gstat package[END_REF], 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 Clim. Past, 12, 1137Past, 12, -1149Past, 12, , 2016 www.clim-past.net/12/1137/2016/ 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 [START_REF] Bartlein | Pollen-based continental climate reconstruction at 6 and 21 ka: a global synthesis[END_REF]. 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 [START_REF] Hewitt | The genetic legacy of the Quaternary ice ages[END_REF], 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 [START_REF] Berger | Long-term variations of caloric insolation resulting from the Earth's orbital elements[END_REF]. 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 [START_REF] Davis | The temperature of Europe during the Holocene reconstructed from pollen data[END_REF]. 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 [START_REF] Sykes | A Bioclimatic Model for the Potential Distributions of North European Tree Species Under Present and Future Climates[END_REF]. 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 [START_REF] Naughton | Present-day and past (last 25 000 years) marine pollen signal off western Iberia[END_REF][START_REF] Naughton | Climate variability across the last deglaciation in NW Iberia and its margin[END_REF] 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 [START_REF] Naughton | Climate variability across the last deglaciation in NW Iberia and its margin[END_REF] is in line with the slight increase in precipitation in all clusters between 14 and 13 ka (Fig. 3).

As described earlier in Europe [START_REF] Renssen | The two major warming phases of the last deglaciation at 14.7 and 11.5 ka cal BP in Europe: climate reconstructions and AGCM experiments[END_REF][START_REF] Heiri | Evidence for cooler European summers during periods of changing meltwater flux to the North Atlantic[END_REF], T jan shows wider changes in amplitude than T jul . The cold to warm transitions that occurred at ∼ 14.7 and 11.5 ka [START_REF] Renssen | The two major warming phases of the last deglaciation at 14.7 and 11.5 ka cal BP in Europe: climate reconstructions and AGCM experiments[END_REF][START_REF] Von Grafenstein | A Mid-European Decadal Isotope-Climate Record from 15, 500 to 5000 Years BP[END_REF] 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 [START_REF] Seppä | July mean temperature and annual precipitation trends during the Holocene in the Fennoscandian tree-line area: pollen-based climate reconstructions[END_REF]. Such trend is more obvious in northern Europe and the Alps, while we rather observe a cooling at lower latitudes [START_REF] Davis | The temperature of Europe during the Holocene reconstructed from pollen data[END_REF]. 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 [START_REF] Cheddadi | The climate of Europe 6000 years ago[END_REF].

Between 6 and 3 ka, areas with low precipitation expand in the Iberian Peninsula (Fig. 3), which allows the expansion of the Mediterranean taxa [START_REF] Naughton | Present-day and past (last 25 000 years) marine pollen signal off western Iberia[END_REF]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 [START_REF] Valiño | Climatic changes since the Lateglacial/Holocene transition in La Mancha Plain (South-central Iberian Peninsula, Spain) and their incidence on Las Tablas de Daimiel marshlands[END_REF].

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 [START_REF] Hewitt | The genetic legacy of the Quaternary ice ages[END_REF][START_REF] Naughton | Present-day and past (last 25 000 years) marine pollen signal off western Iberia[END_REF]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 [START_REF] Branco | Postglacial dispersal of the European rabbit (Oryctolagus cuniculus) on the Iberian peninsula reconstructed from nested clade and mismatch analyses of mitochondrial DNA genetic variation[END_REF][START_REF] Godinho | Genealogy of the nuclear beta-fibrinogen locus in a highly structured lizard species: comparison with mtDNA and evidence for intragenic recombination in the hybrid zone[END_REF]Weiss and Ferrand, 2007;[START_REF] Miraldo | Phylogeography and demographic history of Lacerta lepida in the Iberian Peninsula: multiple refugia, range expansions and secondary contact zones[END_REF]. 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 [START_REF] Ashcroft | Identifying refugia from climate change[END_REF]. However, as shown by large-scale landscape analysis (Carrión et al., 2010b, a) and climate reconstructions [START_REF] Davis | The temperature of Europe during the Holocene reconstructed from pollen data[END_REF][START_REF] Cheddadi | Spatial gradient of temperature and potential vegetation feedback across Europe during the late Quaternary[END_REF], 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 [START_REF] Hewitt | The genetic legacy of the Quaternary ice ages[END_REF]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 [START_REF] Sandel | The influence of Late Quaternary climate-change velocity on species endemism[END_REF], and areas of high velocity are often associated with species extinction [START_REF] Nogués-Bravo | Climate predictors of late quaternary extinctions[END_REF]. Our results indicate higher stability in the southern part of the Peninsula, similar to other studies based on climate data [START_REF] Ohlemüller | Potential source and sink locations for climate-driven species range shifts in Europe since the Last Glacial Maximum[END_REF]. 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 [START_REF] Mee | The ecological and evolutionary implications of microrefugia[END_REF]. Microrefugia are known to occur in the northern areas of the Iberian Peninsula (e.g. [START_REF] Fuentes-Utrilla | Extending glacial refugia for a European tree: genetic markers show that Iberian populations of white elm are native relicts and not introductions[END_REF], 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.

Figure 1 .

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

Figure 3 .

 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 .

 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 .

 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.

  .

	Name	Source	Longitude Latitude Altitude 14 C	Range Biome
	Albufera Alcudia	EPD; Burjachs et al. (1994)	3.125	39.792	11	3000-11 000 Mediterranean
	Algendar	EPD; Yll et al. (1997)	3.958	39.958	80	3000-9000 Mediterranean
	Antas	EPD; Pantaleón-Cano (1997)	-1.792	37.208	14	3000-9000 Mediterranean
	Barbaroxa	Queiroz (1999)	-8.792	38.042	38	3000-7000 Mediterranean
	Cala Galdana	EPD; Yll et al. (1995)	3.958	39.958	80	3000-8000 Mediterranean
	Cala'n Porter	EPD; Yll et al. (1997)	4.125	39.875	81	4000-9000 Mediterranean
	CC-17	Dorado Valiño et al. (2002)	-3.875	39.042	617	3000-12 000 Mediterranean
	Charco da Candieira	EPD; Van der Knaap and Van Leeuwen (1994)	-7.542	40.375	1221	3000-14 000 Mediterranean
	Gádor	Carrión et al. (2003)	-2.958	36.875	1413	3000-6000 Mediterranean
	Golfo	Queiroz (1999)	-9.125	38.542	53	3000-14 000 Mediterranean
	Guadiana	Fletcher et al. (2007)	-7.458	37.292	52	3000-13 000 Mediterranean
	Hoya del Castilho	EPD; Davis and Stevenson (2007)	-0.542	41.292	271	6000-10 000 Mediterranean
	Lago de Ajo	EPD; McKeever (1984)	-6.125	43.042	1744	3000-15 000 Temperate
	Laguna de la Roya	EPD; Allen et al. (1996)	-6.792	42.208	1780	3000-15 000 Mediterranean
	Lake Racou	EPD; Guiter et al. (2005)	2.042	42.542	1906	3000-12 000 Mediterranean
	Las Pardillas Lake	Sánchez-Goñi and Hannon (1999)	-3.042	42.042	23	3000-11 000 Mediterranean
	Navarres 1	EPD; Carrión and Dupré (1996)	-0.708	39.125	278	4000-15 000 Mediterranean
	Puerto de Belate	EPD; Peñalba Garmendia (1989)	-2.042	43.042	622	3000-8000 Temperate
	Puerto de Los Tornos	EPD; Peñalba Garmendia (1989)	-3.458	43.125	893	3000-9000 Temperate
	Quintanar de la Sierra EPD; Peñalba Garmendia (1989)	-3.042	42.042	1546	3000-15 000 Mediterranean
	Roquetas de Mar	EPD; Pantaleón-Cano (1997)	-2.625	36.792	94	3000-6000 Mediterranean
	Saldropo	EPD; Peñalba Garmendia (1989)	-2.708	43.042	645	3000-8000 Temperate
	Sanabria Marsh	EPD; Allen et al. (1996)	-6.708	42.125	1220	3000-14 000 Mediterranean
	San Rafael	EPD; Pantaleón-Cano (1997)	-2.625	36.792	94	3000-15 000 Mediterranean
	Santo André	Santos and Sánchez-Goñi (2003)	-8.792	38.042	38	3000-15 000 Mediterranean
	Siles	Carrión (2002)	-2.542	38.375	1246	3000-15 000 Mediterranean
	Padul	Pons and Reille (1988)	-3.708	37.042	1236	5000-15 000 Mediterranean
	Lourdes	Reille and Andrieu (1995)	-0.042	43.042	727	3000-15 000 Temperate
	Monge	Reille and Andrieu (1995)	-0.042	43.042	727	3000-14 000 Temperate
	Moura	Reille (1993)	-1.542	43.458	40	4000-12 000 Temperate
	Banyoles	EPD; Pèrez-Obiol and Julià (1994)	2.708	42.125	172	3000-15 000 Mediterranean

www.clim-past.net/12/1137/2016/ Clim. Past, 12, 1137-1149, 2016

Clim.Past, 12, 1137Past, 12, -1149Past, 12, , 2016 www.clim-past.net/12/1137/2016/

Data availability

The distribution data for Mediterranean species [START_REF] Legind | Global Biodiversity Information Facility[END_REF] were downloaded from GBIF portal (http://www.gbif. org). The raw fossil pollen data [START_REF] Leydet | European Pollen Database[END_REF] is accessible for most used sites at the European Pollen database (http://www.europeanpollendatabase.net). For the remaining sites, data is available by request to the author of the respective publication (see Table 1 for a description of the data sets and sources). The climate data [START_REF] Hijmans | WorldClim -Global Climate Data[END_REF]) is accessible at the WorldClim website (www.worldclim.org). The output reconstructions, anomalies and variance maps are available in NetCDF format as supplementary material.

The Supplement related to this article is available online at doi:10.5194/cp-12-1137-2016-supplement.

Acknowledgements. Pedro Tarroso was funded through 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