Habitats hold an evolutionary signal of past climatic refugia

Climatic refugia have often been associated with hotspots of richness and endemism, and identified on the basis of molecular or paleobotanical information. Here, we apply a phylogenetic analysis to 18,000 plant communities distributed across the Pyrenees, a south European mountain range, to identify climatic refugia from imprints of relictuality inferred from species’ evolutionary distinctiveness (ED). We produced a genus-level phylogenetic tree to calculate the standardized mean ED value of plant communities (cED). Then, we explored which habitats concentrate the plant communities with the highest cED and the interrelated effect of past (long-term climatic stability) and present (topographic and spatial position) factors. Results show strong differences of cED among habitats: forests ranked first, followed by some open habitats like high altitude wetlands. Climate stability and roughness positively influenced cED. A weak negative association resulted between the two diversity measurements (richness and endemism rate) and also with cED. We propose that forests acted as “mobile refugia” during the glacial-interglacial periods, supported by paleoenvironmental reconstructions revealing continuous presence at regional scale of key broadleaved trees at that time. Azonal habitats like the endemic-poor humid communities at high elevation would have also played an important role as more permanent microrefugia. Our approach identifies a variety of habitats and plant assemblages that have successfully withstood past climate change in different ways, and therefore would hold an important evolutionary potential to cope with current climate change. Given their potential role in preserving biodiversity, they should be integrated in future conservation agendas.


Introduction
The identification of past climatic refugia -areas that provided suitable climatic conditions for species to retreat or withstand during unfavorable climatic periods over geologic time-, has been one of the most interesting scientific challenges in the last decades (e.g., Bennet et al. 1991;Médail and Diadema 2009;Tzedakis et al. 2013). Identifying past climatic refugia, however, is far from being simple, because the use of classical molecular or paleobotanical methods (phylogeographic and genetic population structure analyses, fossils, macro-charcoal, pollen-stratigraphical data) is costly and cannot be extensively applied. Despite the difficulty, the identification of climatic refugia is of primary interest in the current global change scenario because they host lineages that have managed to persist over great spans of time in evolutionary terms. Deciphering the biological features of those surviving species, along with the environmental factors of the locations where they occur, is therefore crucial to foresee the vulnerability of current biodiversity in the face of climate change, and for prioritization in conservation planning (Gentili et al. 2015a;Keppel et al. 2015;Morelli et al. 2020).
Climatic refugia have been traditionally associated with local or regional climatic stability within large areas, as well as with hotspots of diversity and endemism, based on the assumption that stable climates enabled the survival of both ancient taxa and new ones from rapidly diversifying lineages (see review in Harrison and Noss 2017;Keppel et al. 2018). This process has been supported by some studies: small-ranged endemics tend to concentrate in areas of reduced climatic oscillations since the Pleistocene (e.g., Fjeldsa et al. 1999;Sandel et al. 2011), or in unusually cold areas spotted in a warmer matrix in Europe (Ohlemüller et al. 2008). High species richness, however, may be the result of different processes: low extinction rates preserving diversity from old radiations, recent diversification bursts resulting in young and rapid radiations, or a mixture of both processes. Moreover, past events might have blurred current imprints, as both warming and cooling periods operated repeatedly over geological times, and thus modern species either suffered strong geographic range fluctuations during Pleistocene waxing-and-waning cycles, or endured by adapting to new conditions (Smyčka et al. 2018). Therefore, it is not always easy to establish how past climatic events affected current molecular or richness patterns. Smyčka et al. (2018) have recently tried to unravel the main drivers of current patterns of plant endemism in the Alps using a phylogenetic and spatially explicit approach, and concluded that different kinds of evolutionary processes and climatic refugia co-occurred.
The concept of "refugium" evokes the persistence of taxa through major environmental changes at geological scales, and in this context phylogenetic methods constitute a powerful tool because they reach back in time to integrate information about evolutionary relationships, lineage ages and phylogenetic isolation. Although contemporary environmental heterogeneity plays a key role in shaping patterns of richness and endemism (Stein et al. 2014;Irl et al. 2015), historical factors have also left an imprint on them (Jezt et al. 2004). For instance, recent studies have shown that current phylogenetic endemism and diversity patterns in North American forests are influenced by the paleoclimate legacy of the glacial-interglacial climatic variability (Ma et al. 2016), and that tree assemblages with large phylogenetic age differences were associated to relatively high long-term climate stability (Feng et al. 2017). Likewise, phylogenetic studies have revealed that mature radiations and paleo-endemisms occur in climatically and geologically stable environments (Linder 2008), whereas neo-endemisms tend to concentrate in areas of high topographical relief in a Mediterranean hotspot (Molina-Venegas et al. 2017). In this context, we would expect an accumulation of evolutionary distinct lineages in climatic refugia because extinction rates are expected to be lower in scenarios of lower oscillations in temperature and moisture.
Here we propose to incorporate such evolutionary perspective to identify potential refugia by focusing on the phylogenetic singularity (sensu Pavoine et al. 2005) of local assemblages. So far, most attempts to identify refugia have focused on physical and permanent well-defined geographical areas where genetic lineages of analyzed species have persisted through Tertiary or Quaternary climatic fluctuations (i.e., Médail and Diadema 2009). Instead, in this study we will investigate imprints from current plant communities across habitats, whose species are the result of environmental filtering through time (Cadotte and Tucker 2017). The effect of major changes or fluctuating climatic events should have promoted both extinction and diversification, leaving an imprint in current configurations of plant assemblages with different amounts of old and new lineages. In order to identify which plant communities have accumulated the highest concentration of singular lineages through time we use the evolutionary distinctiveness metric (ED; Isaac et al. 2007), a score of evolutionary uniqueness that measures how isolated a species is within the Tree of Life. ED represents species' relative contributions to phylogenetic diversity and disentangles the value of relictual species from recently evolved species even though both types may be rare or endemic (Byrne 2003). Then, we use a community-based measure of ED (cED; see also Morelli et al, 2016a, b), as an index of refugial character to explore which habitats comparatively shelter a higher proportion of phylogenetically more distinct species.
We look for these refugial imprints across a south European range considered a biogeographic crossroad: the Pyrenees. This E-W oriented mountain chain was at the south border of the strongest glaciation activity, becoming a physical barrier for species to move easily between the Iberian Peninsula (a classical southern refugium), and the rest of Europe during the glacial-interglacial period. As a consequence of their complex topography, mountain ranges can be envisaged as environmental mosaics with a wide range of microclimatic conditions and patches of vegetation, offering the possibility for many different species and habitats to fulfill particular requirements and migrate altitudinally to track suitable conditions when the regional or continental climate changes (Rahbek et al. 2019). Such heterogeneity, typical of mountain landscapes, seems to have played a key role in promoting local or regional persistence of organisms not only during the deep past, but also in the recent past as demonstrated by Suggitt et al. (2018), who found lower extinction rates over the twentieth century for plants and insects in areas of higher topographic variation. The combination of an intense glaciation-deglaciation activity due to its southern latitudinal position, along with the environmental patchiness, makes the Pyrenean range an excellent area to extensively prospect the existence of refugia from imprints left in current plant communities, and test their occurrence with the support of available palynological evidences (González-Sampériz et al. 2005, 2017. By using the largest available database of plant inventories in the Pyrenean chain (~ 18,000 phytosociological relevés, ~ 360,000 plant records), we first explore which habitats concentrate species with high ED and exhibit the highest richness and proportion of Pyrenean endemisms, and then we test to what extent these three diversity variables are correlated with historical (climatic stability, probability of being under the ice during the "Last Glacial Maximum"-LGM), ecological and topographic factors (relative position in the mountain chain, terrain roughness, aspect, and altitude). Such information is used to test the following hypotheses: (1) long-term climatic stability and ruggedness promoted the persistence of current "relictual" plant assemblages (i.e., with a high proportion of species with high ED), and (2) climatic refugia allowed not only the persistence of old lineages but also the diversification of new ones, thereby resulting in positive correlations between 1 3 cED, richness and endemism rate. To validate our results with past evidence, we integrate Pyrenean palynological data from last Glacial and Lateglacial sequences (between ca. 30-10 ka cal BP). We finally discuss in which way our findings can be applied to improve conservation management in face of the current climate emergency.

Study area: the pyrenean flora
The Pyrenees spread over ca. 50,000 km 2 between the Cantabric and Mediterranean seas (430 × 150 km approximately), contain more than 150 peaks over 3000 m a.s.l., and constitute a natural border between Spain and France ( Fig. 1). They encompass three distinct vegetation regions (Eurosiberian, Mediterranean and Alpine), and include contrasting landscapes related to aspect, continentality, bedrock type, landform and ancient anthropic land use. The northern side and the western part are strongly influenced by the Atlantic Ocean and high annual rainfall, with a submontane belt supporting mixed deciduous forests, acidophilous oakwoods, beech forests, and fir forests, heaths and meadows. The southern side shows a more complex altitudinal zonation: the Mediterranean belt is substituted altitudinally by the submontane (marcescent oak woods with sclerophyllous oakwoods and xerophilous scrublands), montane (Scots pinewoods together with meso-xerophilous grassland and box scrub, and spotted beech and fir forests on moister slopes), subalpine (forests of Fig. 1 The Pyrenean range, located between the Iberian Peninsula and the rest of Europe. Grid of 10-km UTM cells used to geolocate plant inventories, altitudinal range according to colors, and ice cover extension during the Last Glacial Maximum mountain pine and grasslands maintained by summer cattle), and alpine formations (vegetation mosaics dominated by rocks) (Ninot et al. 2017).

Data gathering: plant communities, habitats and environmental variables
In 2019 we downloaded from the Iberian and Macaronesian vegetation database (SIVIM, Font et al. 2010) all the Pyrenean phytosociological relevés located above 400 m a.s.l. with a geospatial reference of UTM of 10 km side or finer precision (UTM: Universal Transverse Mercator; the most frequent geolocation unit adopted by botanists; see Fig. 1). They were collected over the last 60 years. The phytosociological method is based on full surveys of plants typically co-occurring within habitats, and by analyzing these natural assemblages we avoid the bias associated to the use of herbarium records at large scale (i.e. low completeness and not necessarily co-occurring species within grid cells; see also Pardo et al. 2017). Throughout this paper we will use the term inventories instead of relevés because we only used the list of plants and not their abundance. Inventories without information of altitude, aspect, or containing less than five plant records were discarded. In order to avoid taxonomic conflicts, taxonomic names were standardized at the species level according to the Atlas of the Pyrenean Flora (FLORAPYR; Gómez et al. 2017). The final data set was made of 17,923 inventories and 360,306 plant records of 2,537 plant species (about 74% of the Pyrenean flora) distributed all over the Pyrenees and encompassing all relevant environmental heterogeneity (Fig. 2).
Each community was then classified within one of the 14 broad and typical Pyrenean habitats according to the European Nature Information System habitat classification (EUNIS; https:// eunis. eea. europa. eu/; Table 1 and Fig. 2). Habitats with very high anthropogenic influence such as ruderal communities and meadows were discarded. To better accommodate the classification hierarchy of this system to the habitat diversity of the Pyrenees, we associated them either to the first or second EUNIS levels based on phytosociological alliance or expert knowledge (for example, F3 and F4 types were grouped). The number of inventories used in this study per EUNIS habitat ranged between 344 (inland surface waters) and 3264 (alpine and subalpine grasslands) (Table 1).
Species richness corresponded here to the number of plants in each inventory, and endemism rate to the proportion of taxa restricted to the Pyrenean range. The geolocation of each inventory was defined as the longitude and latitude of the centroid of the corresponding 10 km-UTM grid cell. Aspect and altitude were part of the characterization of the inventories stored in the original SIVIM database. Aspect was grouped in two categories depending on the dominant component: north-facing (north, northeast and east facing slopes) and south facing (south, southwest and west slopes). In addition, for each 10 km-UTM grid cell we computed past and present variables that, a priori, could be related to the persistence of old lineages: (1) the climatic stability (the opposite concept to the also referred as "temperature anomaly" in the literature), which was estimated as the difference between the mean temperature over the 1981-2005 period according to the CHELSA model and the mean temperature during the Last Glacial Maximum (LGM) according to Karger et al. (2017); (2) the probability of being covered by ice during the LGM, which refers to the percentage of ice sheet cover according to the Community Climate System Model NCAR-CCSM4 (a coupled climate model for simulating the earth's climate system, composed of four separate models simultaneously simulating the earth's atmosphere, ocean, land surface and sea-ice, and one central coupler component; Ehlers and Gibbard 2004); (3) the mean altitude; and (4) the mean Terrain Ruggedness Index (TRI) based on a 25 m resolution digital elevation model. Since 10 km-UTM grid cells might be too coarse to detect the effect of some variables, we repeated all the analysis for the subset of inventories that had been georeferenced at 1 km-UTM grid cells (N = 6,456). Topographic variables were computed using ArcGIS Version 10.5 (ESRI 2011).

Phylogenetic data
We built a genus-level phylogenetic tree of the Pyrenean flora following the pipeline proposed by Roquet et al. (2013). We downloaded from Genbank three chloroplastic regions (rbcL, matK and ndhF) for all available taxa, plus the ITS region for certain families (which we aligned separately by taxonomic clustering). We were able to retrieve sequences of these markers for all genera except Cytinus, Ptychotis and Xatardia, which represent 0.001% of our taxa dataset and none of them is a dominant species. All coding sequence clusters were aligned with MACSE (Ranwez et al. 2011) and non-coding ones with MAFFT (Katoh and Standlye 2013) and trimmed all alignments with TrimAl (Capella-Gutiérrez et al. 2009). We concatenated all alignments to obtain a supermatrix. We then conducted maximumlikelihood (ML) phylogenetic inference analyses with RAxML (Stamatakis 2014), applying the most appropriate partitioning scheme and substitution model obtained with Parti-tionFinder (Lanfear et al. 2012) and a supertree constraint at the family-level obtained with the online software Phylomatic v.3 (tree R20120829). Specifically, we performed 100 independent tree searches. The best ML tree (the one with the highest probability) was dated applying the penalized likelihood method in treePL (Smith and O'Meara 2012) and the following node calibrations: we fixed the node corresponding to the ancestor of eudicots at 125 Ma based on the earliest eudicot fossil (Hughes and McDougall 1990) and applied minimum age constraints to 15 nodes based on fossil information extracted from Smith et al. (2010). To deal with unknown within-genera structures, we simulated 10 scenarios of within-genera random branchings using a Yule process as implemented in the R library apTreeshape (Bortolussi et al. 2012). These 10 trees represent a distribution of possible hypotheses about evolutionary relationships in our dataset sensu Rangel et al. (2015).

Phylogenetic structure analysis
For each species of the Pyrenean flora, ED was computed as the sum of the branch lengths from the species tip to the root of the tree divided by the number of species subtended to each branch (Redding and Mooers 2006) using the equal.splits method from the evol.distinct function of the picante R package (Kembel et al. 2010). In order to have a communitylevel evolutionary distinctiveness measure (cED thereafter) we summed the ED values of the species present at each inventory (see also Morelli et al. (2016a, b)). Given that this measure is highly correlated with the species richness of each community, we compared the observed cED with a null model. The community matrix was randomly shuffled 999 times while keeping species richness at each inventory constant, and cED was computed for each inventory through every iteration. We then calculated the standard size effect of the observed cED (ses.cED) as follows: where cED obs is the observed cED at each inventory, cED rand is the mean value of the distribution of cED values obtained with the null model, and cED sd is the standard deviation of such distribution. Positive values of ses.cED indicate that the community is composed by species more evolutionarily distinct than expected by chance, and therefore we assume it has a "refugial" character. As we had ten different versions of our phylogeny, depicting different evolutionary hypotheses, we repeated the aforementioned process for each tree and averaged the results in order to include the different possibilities. In order to avoid the influence of extreme values in our data we discarded inventories with ses.cED values ses.cED = cED obs −cED rand ∕cED sd 1 3 that were 1.5 times the interquartile range lower or higher than the 25% or 75% quartiles respectively. Then, to decipher which habitats accumulated species of higher ED, we ranked all inventories based on their ses.cED values and then calculated the number and relative proportion of inventories belonging to each habitat that fell within the highest 10% values of ses.cED. To facilitate the reading, throughout this paper we will use cED as the standardized community ED (ses.cED). Metrics derived from regional phylogenies might suffer from biases due to incompleteness (Molina-Venegas and Lima 2021). However, the phylogeny used in this study includes more than 99% of genera of the Pyrenean flora and the resulting metrics were compared between habitats. In consequence, the results from our analysis are relative to each other, which avoids the possible negative effects of that bias.
Phylogenetic structure analyses were performed both with and without Gymnosperms and Pteridophytes because they may strongly influence phylogenetic structure indices, given their phylogenetic origin is much older than Angiosperms' age (Feng et al. 2017). We only present here the ones obtained for Angiosperms, but similar results were produced including the other two groups.

Statistical analysis
We evaluated differences between habitats in terms of cED, richness and endemism rate, as well as the effects of temperature stability, percentage of ice cover in UTM cells during the LGM, altitude, aspect and TRI, using three different generalized additive regression models (GAM; normal distribution for cED and richness; and binomial for endemism rate) with the mgcv package (Wood 2011) in R (version 4.0.4; R Core Team 2021). We also included a spline with the effect of the interaction between latitude and longitude of UTMs where inventories were localized, as we expected some spatial structure in our data. The 95% confidence intervals of estimated marginal means predicted by the model were used to test differences between habitats. For comparative purposes, we also produced maps with the estimated mean cED, richness and endemicity rate for 10-km UTMs from the same set of inventories and GAM models.
To test the hypothesized association between refugia, richness and abundance of endemics, we analyzed the correlation between cED, species richness and endemism rate using Pearson's correlation coefficient.

Results
Different habitats stand out in each of the three diversity measures analyzed in this study ( Fig. 3 and Table 1). Habitats hosting the largest number of species were woodland fringes and clearings and tall forb stands (E5), followed by broadleaved deciduous forests (G1), and dry grasslands (E1). Dry grasslands (E1), deciduous forests (G1), and alpine and subalpine grasslands (E4) showed the largest average values of species per inventory. Rocky habitats, i.e. screes (H2) and inland cliffs, rock pavements and outcrops (H3), hold the highest endemism rate. Broadleaved forests (G1, G2) and mires, bogs and fens (D) ranked first in terms of cED.
Models significantly explained cED, richness and endemism rate, with deviances ranging from 51.7% (cED) and 33.2% (richness). Between 28% (endemism rate) and 62% (richness) of such deviance was absorbed by the habitat variable. All included factors except ruggedness significantly contributed to explain cED in the 10-km UTM model, whereas richness was significantly affected by all factors and endemism rate was not affected by ruggedness nor by historical factors (Table 2). Overall, altitude always had a very low but significant effect on the three dependent variables: negative for cED and richness, and positive for endemism rate. Other factors also exerted contrasted effects: temperature stability and south exposure, for example, affected negatively cED and positively richness, and the opposite effect was exerted by the percentage of ice  Table 1). Right: Mean ± sd of richness per inventory (above), and estimated marginal mean ± sd (from the GAM) of standardized community evolutionary distinctiveness (cED) for each habitat type (below)

Table 2
Outputs of the Generalized Additive Models (GAM) for standardized community evolutionary distinctiveness (cED), richness, and endemism rate (only angiosperms included). For habitats, inland surface waters is the reference level. In all cases geolocation was significant (p value < 2*10 -16   1.155 *** cover in the past. South exposures also affected negatively the endemism rate. The particular geolocation of inventories also had a significant effect on the three diversity variables but resulted in different patterns (p < 0.0001; Fig. 4): the highest cED values corresponded to the perimeter of the study area (NE and S-central part), whereas the lowest cED values concentrated in the central (highest) part of the chain, the one of lowest endemicity and highest richness. To unveil the variables whose effect could be scaledependent, we re-run models with a subset of 6,456 inventories geolocalized at 1-km UTM. Results concerning the positive or negative effects of each factor were coherent in both scales (10-km and 1-km), but ruggedness emerged as significantly positive for cED and endemism rate, climate stability did no longer affect cED and richness, and south exposure did not affect richness either.
Finally, we explored the relationship between all three diversity measures: richness, endemism rate and cED. Although the strength of the relationships was weak, it was always negative: species richness with cED (r = − 0.02; p < 0.01) and endemism rate (r = − 0.26; p < 0.001), and also between the two last measures (r = − 0,10; p < 0.001).

Fig. 4 Left: Heatmaps
showing predicted values of richness, endemism rate and standardized evolutionary distinctiveness across longitude and latitude variables in the Pyrenean range, resulting from GAMs that include past and present factors. Right: Estimated spatial distribution of standardized community evolutionary distinctiveness (cED, above), richness (middle), and endemism rate (below) in the Pyrenean range (different Generalized Additive Models were set for each map; see text for further explanation). Predicted values follow the color scale 1 3

Discussion
Building on an eco-evolutionary approach, we identified putative climatic refugia in the Pyrenean range from imprints of relictuality in plant assemblages, based on the mean value of the Evolutionary Distinctiveness of species composition (cED). This approach successfully revealed clear differences in evolutionary originality among habitats across the Pyrenean range, and the important influence of climatic stability and ruggedness on cED. Our results also showed that different factors contributed in different ways to explain the three diversity measures analyzed (cED, richness, and endemism rate), and the existence of weak negative relationships between them.

Contemporary and historical factors associated to high cED
Our study demonstrated that habitats hold different amounts of cED in the Pyrenean range. Widespread mid-altitude broadleaved deciduous woodlands harbored not only high species richness but also the highest cED, followed by broadleaved evergreen woodlands at the lowest altitudes and the southern periphery. On the contrary, the very rich grasslands in terms of number of species (both dry, and (sub)alpine grasslands at high altitude), exhibited the lowest cED. Such differences among habitats, often associated to different altitudinal belts, would explain the negative effect of altitude on cED.
Climatic stability is usually associated with the existence of refugia (Fjeldå et al. 1999;Jezt et al. 2004;Harrison and Noos 2017). Our results do support such hypothesis, as relatively small differences in temperature between present day temperatures and the ones during the LGM had a significant positive effect on cED. Our model also captured the positive effect of roughness on cED too, but only at fine spatial scale. The resolution of 10-km UTM models resulted too coarse for a mountain area like the Pyrenean range, and clearly not appropriate to characterize sampling units (inventories) as small as 4-400 m 2 . The effect of climate stability, however, disappears at finer scale, when roughness emerges. More detailed information on ruggedness of sampled communities could probably offer a neater picture of its importance, as this variable has been shown to be transcendental for the probability of persistence of small organisms (Suggit et al. 2018).
Our results also indicated a slightly higher cED in communities at northern compared to southern expositions, which may be related to the higher frequency of cold-adapted (Quaternary) than warm-adapted (Tertiary and Mediterranean) plants in the Pyrenees. Nevertheless, the map of predicted cED (Fig. 4) depicts a bit more complex scenario than just S versus N sides, where the peripheral part stands out over the central part. Our pattern fits well the proposal of Médail and Diadema (2009), who in a comprehensive review situated 33 out of the 52 putative refugia in Mediterranean submontane and mountain margins, and reported numerous refugia both in northern and southern parts (their review did not include the whole Pyrenean range but only the southern part). The paleoenvironmental history of southern European mountains like the Pyrenees has been very complex and involved climatic events occurring since the Tertiary through the Quaternary (Carrión and Leroy 2010 and references therein). Their alpine flora retreated to southern or lower mountain areas during glacial periods, and survived at high elevations after ice retreated over the warmer interglacials, with the consequent upwards colonization by forests. Several cold-adapted plant lineages radiated during the glacial-interglacial cycles, which promoted diversification events for a number of Mediterranean lineages through isolation and allopatric differentiation (Birks 2016). For these reasons, the Pyrenean chain is expected to hold a heterogeneous array of climatic refugia for both warm and cold-adapted plants. In fact, available palynological data from last Glacial and Lateglacial pollen sequences from the Pyrenees show a mosaic of vegetation patches resulting in a heterogeneous landscape, with abundance of steppes during the most intensive cold-arid episodes in both non-glaciated mountains and lowlands, and continuous occurrence of woody taxa even at mid-and relatively high areas (Jalut et al. 1992;González-Sampériz et al. 2005, 2017.

Different kinds of refugia for different kinds of plant species
Our results highlight the role of all the three kinds of forests as past climatic refugia. Needless to say, canopy structure creates a microclimate that contributes to stabilize relative humidity in summer droughts, cancels frost in wintertime, and provides shelter from strong winds, which altogether buffers the impact of large-scale climatic changes (von Arx et al. 2013;de Frenne et al. 2013;Davis et al. 2019). Such local buffering may be the critical mechanism for understory plants to persist till present days even under harsh regional environmental conditions. In fact, recent changes in the abundance of understory plants have been proven to be more closely related to the rate of subcanopy microclimate than macroclimate change (Zellweger et al. 2020).
Palynological findings and studies of fossil macrocharcoal (Uzquiano 1992;Alcolea 2017;Aranbarri et al. 2016) support the proposal of Peñalba (1994) that northern Iberia was covered by a patched network of tree populations during the last glacial age González-Sampériz et al. 2010, 2017. The presence of deciduous forests, remnant fragments or even small clonal stands would have then allowed ground-layer communities to survive over periods of extreme temperatures. Besides, paleoenvironmental reconstructions from Lateglacial pollen records have shown relatively rapid responses of tree formations to deglaciation and climatic amelioration, as well as continuous presence (although at low proportions) of some taxa during glacial times. Pollen grains of the main species with high ED such as Tilia, Populus, Alnus, Fagus, and deciduous and evergreen Quercus, have been recorded in many palynological sequences of both southern and northern Pyrenees, including sites at different altitude and latitude (Table 3). Our analyses indicated a slightly higher cED in communities at northern compared to southern expositions, which fits the continuous presence and higher frequency of Tilia and Fagus pollen grains in northern palynological sequences of the Pyrenees than in southern ones (see Table 3, only a selection of records are included here). In addition, in mountain areas of continental central Europe resembling our study area, other indicators such as macroscopic charcoal have also demonstrated the existence of long-term, stable tree populations during the coldest intervals of the last glacial period, such as the LGM and Heinrich events (Willis et al. 2000;Stewart 2003;Birks and Willis 2008). Our approach, thus, highlights the role of forests as climatic refugia based on their proven existence during the glacial times. Forests allowed the persistence of some old lineages that were dragged along with tree dispersals to track the most suitable environmental conditions during climatic changes, or transported through water courses in the case of riparian forests (Huntley and Birks 1983). In sum, our findings fit well with palynological records and also with putative refugia identified by Médail and Diadema (2009) in the Mediterranean Basin, where several mesophilous trees in relatively low altitudinal areas were proposed, notably in locally moist and warm sites such as valley bottoms and wetlands. In those places, local topography permitted rapid altitudinal displacements as response to climatic changes. Whereas the refugium concept is Table 3 Pollen records of taxa with high ED in selected palynological sequences with full Glacial and/or Lateglacial chronologies from both northern and southern Pyrenees, organized by altitude Note that, except for Populus (which usually has poor taxonomic preservation and it is not frequent to find it in palynological sequences) all taxa have been recorded in both northern and southern Pyrenees (despite differences regarding altitude). References: Hyena coprolites of Gabasa cave (González-Sampériz et al. 2005  mostly thought in terms of a physical, spatial location, our approach widens the perspective to consider some vegetation formations as refugia themselves, thanks to their role in buffering climatic changes for understory plants. As a result, forests could be considered as mobile refugia through space. We also found that some kinds of open habitats hold an important concentration of communities with high cED: mires, bogs and fens (D), temperate and Mediterranean-montane scrub and heathlands (F3-F4), and inland cliffs, rock pavements and outcrops (H3). High altitude wetlands have been identified as present climatic refugia for cold-adapted organisms in alpine systems in the recent review of Brighenti et al. (2021), again on the basis of stronger decoupling from air warming than terrestrial environments due to greater heat capacity of water. Similar to high altitude wetlands, rocky habitats are patchy and have low species richness, but contrary to them they are very stressful for plant life and ranked first in endemism rate. Given their azonal character, Birks and Willis (2008) considered that rocky habitats sheltered different degrees of "relictness", which would translate into intermediate values of cED. Our study showed that they hold intermediate cED values in the case of inland cliffs, rock pavements and outcrops. The internal environmental heterogeneity of those rocky habitats distributed over a very long altitudinal gradient, as well as steep N and S aspects that make them ranging from hidden to very exposed to solar radiation deserves further analysis,which could unveil masked patterns now aggregated. As an example, using high-resolution infrared thermometry and sensors of temperature with data-loggers,  demonstrated the buffering effect of daily and seasonal temperatures in rocky crevices and cliffs where two relict Pyrenean plants have survived till present day: a Tertiary warm-adapted and a Quaternary cold-adapted. Several authors have stressed the importance of present-day cold rocky landforms as microclimate refugia (''cryptic refugia'', Birks and Willis 2008) in mountains all over the world, for species that were much more widespread during the LGM (Gentili et al. 2015a(Gentili et al. , b, 2020Brighenti et al. 2021). Gentili et al. (2020) recently demonstrated the widespread availability of isolated microhabitats for typical alpine plants, which could increase by deglaciation in the current warm period. The patchiness of these habitats may have fostered the speciation of endemic plants through allopatry, and the fewer competitive interactions found there may have favored their persistence, which could explain the high endemism proportion in rocky habitats.

Relictuality, richness and endemicity in the Pyrenean flora
Our analysis revealed that patterns of richness and cED do not match well; i.e., the presence of old, poor-branched lineages is not associated with hotspots of richness. Interestingly, not a single historical or contemporary factor exerted the same positive or negative effect for these two diversity metrics in our models. Our results therefore, do not support the hypothesis of refugia necessarily sheltering high richness (see review in Keppel et al. 2012). cED was also not found to be positively related to endemism rate, thus not supporting the "endemic hotspot-refugia" pattern found by other studies (Svenning and Skov 2007;Médail and Diadema 2009;Harrison and Noos 2017). These result are not new. In the Mediterranean Basin, for example, Médail and Diadema (2009) already placed some putative refugia outside diversity hotspots. And at global scale, Orme et al. (2005) found that hotspots of species richness did not overlap hotspots of endemism, suggesting that different mechanisms are responsible for the maintenance of different aspects of diversity. Looking for climatic refugia, Smyčka et al. (2018) found that the ED of species was only weakly associated with endemic status in the Alps, and concluded that different kinds of evolutionary processes and climatic refugia cooccurred. Clearly, the association between diversity patterns and refugia seems to be complex and both scale-and context-dependent, and more in-deep analysis would help to interpret the fingerprints of past events.

Conservation implications
Phylogenetic measures of diversity broaden our comprehension of the processes underlying current biodiversity patterns, going well beyond those focused on the visible tips of ecosystems (e.g., species richness and endemism rate). For that reason, they are becoming increasingly relevant to applied conservation, particularly for prioritizing species in conservation planning (e.g., Pavoine et al. 2005;Isaac et al. 2007), and for identifying areas that maximize the protection of the multiple facets of biodiversity (e.g., Pollock et al. 2017). Our approach identifies which plant assemblages have coped better with past climatic changes, and reveal that the loss of different kinds of habitats implies a differential loss of distinctiveness and evolutionary potential to face climatic changes. These results may have a straightforward application into the future because areas of high distinctiveness should be considered in the conservation agenda as a way of maximizing evolutionary potential and future ecosystem functionality in a sustainable way (Keppel et al. 2015). Combined with the pressure of threats, our analysis can help us to prioritize the ones that deserve urgent actions and their feasibility to be undertaken. Whereas rocky habitats are expected to be the least affected by ongoing global change because the low probability of being transformed, the highly distinctive alpine wetlands are very exposed to climate change because of its association to areas where glaciers melt, and other runoff waters. Forest management, on the other hand, is a traditional practice that can be easily used to preserve ED.
Identifying areas where contemporary climate change is milder, the so-called "slow lanes" for resident biodiversity (Morelli et al. 2020), should constitute the first and most sensible step of mitigation strategies (Keppel et al. 2015). Our study proposed the use of datasets of inventories to identify plant communities and habitats of high value in the context of climate change. It could be used to create environmentally heterogeneous networks to cope with different kinds of climatic uncertainty, allowing species of different thermal affinities to remain, move and/or transform as the climate does. By considering the most unique and successful ancient components of plant assemblages, our approach provides a new layer of information in the framework of the climate-change refugia management concept (Morelli et al. 2016b), contributing to a more integrative and sustainable preservation of biodiversity.