Genetic and morphological diversity in Armeria ( Plumbaginaceae ) is shaped by glacial cycles in Mediterranean refugia by

Fuertes Aguilar, J., Gutiérrez Larena, B. & Nieto Feliner, G. 2011. Genetic and morphological diversity in Armeria (Plumbaginaceae) is shaped by glacial cycles in Mediterranean refugia. Anales Jard. Bot. Madrid 68(2): 175-197. Little is known of the direct effects of Quaternary glaciationdeglaciation cycles in plants within southern European refugia. This study, centered in the Sierra Nevada (S Spain), used RAPD and morphometric data from 36 populations of Armeria (Plumbaginaceae) from five taxa belonging to three species that are endemic to that region: A. filicaulis subsp. nevadensis, A. fili caulis subsp. trevenqueana, A. filicaulis subsp. alfacarensis, A. splendens, and A. villosa subsp. bernisii. The results based on genetic analyses at the population level (AMOVA, genetic diversity, genetic distance) and genetic and morphological analyses at individual level (haplotype phenetic distance, PCO, morphometrics) indicate that: (1) genetic diversity decreases with altitude, probably as a result of the postglacial recolonization processes, except in some secondary contact zones between taxa; (2) gene flow among interspecific populations, most likely facilitated by contraction of vegetation belts, led to the formation of hybrid taxa; (3) genetic distances among populations provide a useful basis for studying scenarios with frequent interspecific gene-flow since it allows distinguishing eventual cases of introgression from hybridogenous taxa.


Introduction
Glaciations have shaped present-day plant species distributions in the Northern Hemisphere (Hewitt, 1996(Hewitt, , 2000) ) through four main direct effects on populations: isolation, expansion, migration and extinction (Comes & Kadereit, 1998).During glacial periods, species that were adapted to arctic-alpine conditions colonized large extensions of tundra and steppe south of the ice sheets and alpine enclaves, while those adapted to temperate conditions were confined to southern refugia (Ferris & al., 1999).In contrast, interglacial periods favoured the recolonization of newly open habitats in deglaciated areas, thereby inducing migratory movements towards northern areas that allowed the admixture of previously fragmented populations (Gabrielsen & al., 1997).In some of the southern European areas, including high-altitude massifs, a more complex pattern is added to these latitudinal refugia dynamics.Due to the great diversity of habitats along high-elevation mountains, altitudinal shifts in plants induced by climatic changes took place at a much more reduced scale than latitudinal changes (Hewitt, 1999).As a consequence of overlapping latitudinal and altitudinal shifts associated with contraction-expansion cycles, southern refugia became extraordinary arenas for plant evolution that fostered a variety of outcomes including fragmentation, competition, hybridization and speciation (Brochmann & al., 1998;Steen & al., 2000).Most of the phylogeographic and genetic research on European refugia has focused on extinction and recolonization processes in the arctic-alpine flora of northern and central glaciated areas (Abbott & Comes, 2004;Ronikier & al., 2008), the identification of central and southern European refugia (Tribsch & Schönswetter, 2003;Schönswetter & al., 2005), and the pathways of recolonization from such refugia into deglaciated areas (Petit & al., 2003).However, little work has been done on how glaciationdeglaciation cycles affected the plants within the refugia of southern Europe (Kropf & al., 2006(Kropf & al., , 2008;;Petit & al., 2002;Nieto Feliner, 2011).
The Sierra Nevada, the highest massif in the Betic system (SE Iberian Peninsula), emerged during the Langhian (earliest middle Miocene) and estimations from sedimentary basins indicate that as early as in the Tortonian-Messinian transition (7.1 Mya), the massif had already reached 2000 m above sea-level (Braga & al., 2003).Since then, despite erosion forces, it has suffered a progressive uplift until its present-day highest elevation at the Mulhacén peak (3480 m).Climatic oscillations in the Sierra Nevada during the last 5 My (Tertiary-Quaternary) are well documented from paleosol evidence within adjacent basins (Günster & Skowronek, 2001) as well as from palynological records in caves and lacustrine sediments (Pons & Reille, 1988;Carrión & al., 2001a).Paleosol studies in periglacial areas reveal that at least three cold Quaternary episodes took place in this massif, identified as Riss, Würm, and Late Glacial, each one with a decreasing glacial activity (Gómez Ortiz & Salvador i Franch, 1996;Simón & al., 2000;Schulte, 2002).As a consequence of its orography, during Quaternary iceages the Sierra Nevada paradoxically contained the southernmost glaciers in Europe and at the same time the surrounding area was a refugial area for plants (Carrión, 2002;Carrión & al., 2003).
Previous studies based on trnL-trnF sequences (Gutiérrez Larena & al., 2002) revealed a shared pattern of plastid haplotypes among species suggesting an interspecific horizontal transfer in Sierra Nevada, with an initial colonization process being the establishment of Armeria splendens.The study of ITS se-quences (Fuertes Aguilar & Nieto Feliner, 2003;Nieto Feliner & al., 2004), revealed the prevalence of a ribotype (R3) that was almost restricted to the Sierra Nevada and which was present in all samples of the five taxa.However, in a few individuals from populations of A. filicaulis and A. villosa R3 co-occurred intra-genomically with other ribotypes, either R1 or R2, which were also detected outside the Sierra Nevada.Both studies discovered a species-independent geographically structured pattern of variation of sequence data that was also partially independent from morphological characters.Moreover, the geographical patterns of the nuclear and plastid sequences were also independent from each other.This finding, together with the emerging evidence of genetic exchange at a local scale between the different species (Nieto Feliner & al., 1998), underlined the necessity of a total genome marker study to detect local introgression and putative late generation hybrids.Due to the independence of morphological variation from DNA sequence variation (Fuertes Aguilar & al., 1999) in southern Iberian Armeria spp., a complementary morphometric analysis of the involved taxa was also considered relevant to help interpret the molecular data.
We chose the genetic fingerprint marker RAPD because it represents a random sample of the entire genome (Williams & al., 1990;Lorenz & al., 1994) and because it had provided sufficient variation in studies of related taxa of Armeria from other areas (Nieto Feliner & al., 2002).Most RAPD is nuclear and thus biparentally inherited in angiosperms (Weising & al., 1995;Rieseberg, 1996).This implies that ad-ditive patterns of such markers are expected in hybrids at the population level, and consequently they can be used to examine hybridization and introgression (e.g.Sale & al., 1996;Martin & Cruzan, 1999).RAPD phenotypes have been widely used to estimate the genetic structure of populations (Steward & Excoffier, 1996;Palacios & González Candelas, 1997;Jordano & Godoy, 2000;Segarra-Moragues & Catalán, 2003).Despite criticisms raised on the basis of reproducibility (Skroch & Nienhuis, 1995;Jones & al., 1997), a certain re-evaluation for this technique came from a review by Nybom (2004), who indicated that estimates of population genetic structure derived from markers such as RAPD, AFLP and ISSR produce very similar values and may be directly comparable.This view has been confirmed by several studies on plants occurring in the Mediterranean region (Morgan-Richards & Wolff, 1999;Golan-Goldhirsh & al., 2004).
Determining how morphological and genetic diversity patterns may have interacted within putative glacial refugia requires a study encompassing several techniques.Using such an approach, the aims of our study were to: (1) explore the genetic boundaries of the Armeria taxa inhabiting the Sierra Nevada and the possible occurrence of gene exchange among them, particularly in the hybridogenous A. filicaulis subsp.nevadensis; (2) assess the genetic structure and patterns of gene flow among populations from the different taxa; (3) compare morphometric and genetic variation across the taxa inhabiting the Sierra Nevada by verifying congruence with that previously obtained from sequence data; (4) interpret the data in the light of the glacial history of the massif to help clarify the biogeographical role of southern European massifs as Quaternary refugia.

RAPD sampling
The sampling for the molecular study included a total of 197 individuals from 36 populations across the whole altitudinal and latitudinal geographical distribution of the studied species (Fig. 1, Table 2).From the Sierra Nevada, we sampled populations of Armeria filicaulis subsp.nevadensis (4 pops.), A. filicaulis subsp.trevenqueana (1 pop.), A. splendens (3 pops.), and A. villosa subsp.bernisii (17 pops.).In addition, three populations of A. villosa subsp.bernisii from the surrounding massifs of Sierra de Gádor, Sierra de los Filabres and Sierra de Tejeda-Almijara were included.Five populations of A. filicaulis from the nearby Sierra de Tejeda-Almijara (subsp.filicaulis) and Sierra de Alfacar (subsp.alfacarensis), and three populations of A. villosa subsp.longiaristata from Sub-betic localities at Cazorla and Alfacar completed the sampling, For each population, samples of leaves from at least 5 individuals were collected and preserved in silica-gel until DNA isolation.One to five individuals per population J. Fuertes Aguilar & al. were collected for the morphometric study and as a voucher for identification.Additionally, 29 individuals from surrounding populations not included in RAPD analysis were included to complete the morphometric study (Table 2).Vouchers were deposited in the MA herbarium.

Primer selection and pilot study
In order to select those primers yielding an informative and reproducible banding pattern, a primer selection test was conducted from a set of 49 oligonucleotides used in previous studies in which withinspecies variability was observed (Swensen & al., 1995;Stewart & Excoffier, 1996;Jordano & Godoy, 2000;Nieto Feliner & al., 2002).Each primer was tested with nine samples obtained from three individuals belonging to three different taxa.Nine of the tested primers were selected based on number (30>n>5), resolution and intensity of bands.To assess the reproducibility of the PCR products, a survey was performed testing each of the nine primers with one individual per population.Two PCRs separated by at least 24 hours were performed and their respective banding pattern compared.As a result of this test of reproducibility, one primer was discarded and eight were chosen for the final study (OPC-06, OPC-15, OPC-16, OPD-12, OPJ-05, OPJ-09 from Operon Technologies, and UBC-111, UBC-489 from University of British Columbia catalogue).For further reliability, bands considered reproducible were restricted to a range between 100 and 1200 bp.

DNA isolation, PCR and RAPD analysis
DNA was isolated from leaf tissue using a modified method of Doyle & Doyle (1987) within 30 days after collection.0.2 g of silica-dried leaf tissue samples were homogenized and added to 0.20 ml of CTAB buffer (1.4 M NaCl, 100 mM Tris-HCl pH 8.0, 20 mM EDTA pH 8.0, 2% CTAB, 1% PVPP, 0.25% ß-mercaptoethanol).Samples were incubated at 65 °C for 60 min, and 0.5 ml of chloroform-isoamyl alcohol (24:1) was then added and the mixture was gently mixed for 10 min.After centrifugation at 15,000 rpm for 10 min, the aqueous layer was removed and the process repeated.The aqueous layer was subsequently precipitated with 2.5 volumes of -20 ºC ethanol.The DNA was spooled, air-dried, and re-suspended in 50 µl TE.DNA was quantified using a Gene Quant II UV spectrophotometer (Amersham Biosciences Europe GmbH, Cerdanyola, Barcelona, Spain).DNA samples were adjusted to a concentration of 35 ng/µl.RAPD reactions were performed in 50-µl aliquots containing 25 ng of template DNA, PCR buffer (50
A. villosa subsp.bernisii  Homology was assigned between bands from different lanes based on size and performed automatically, allowing a 1.5-2 % tolerance depending on primers.This procedure made feasible the comparison among different gels and the use of a uniform criterion for the large amount of processed samples.Variations in intensity between bands of equal molecular weight across samples were not considered.Those bands that could not be reproduced twice in the same individual during the pilot study were discarded.A final 1/0-matrix representing band presence/absence was then generated and exported to .xlsformat for further analysis.

Statistical and genetic analysis of RAPD
A multivariate analysis of the binary RAPD matrix by Principal Coordinates Analysis (PCO) was performed using NTSYS-pc 2.11a (Rohlf, 2002).To avoid bias produced by false negatives the Dice coefficient (S) (Dice, 1945), also known as the Sørensen coefficient, was used to construct the individual pairwise similarity matrix, as suggested by Wolfe & Liston (1998).The Jaccard coefficient was also calculated but results did not differ significantly (data not shown).A genetic distance version of the Dice coefficient was generated transforming the similarity matrix in NTSYS-pc using the formula D = (1-S) 1/2 , as pro- posed by Legendre & Legendre (1998).This distance matrix was then imported into PAUP*4.0b10(Swofford, 2002) to produce a Neighbour-Joining tree (Saitou & Nei, 1987) with support values of branches estimated through 1000 bootstrap analyses.The partitioning of the genetic variance within and among popu lations was evaluated through the analysis of molecular variance (AMOVA), as presented in the package ARLEQUIN 2.0 (Schneider & al., 2000).The analysis was performed at three hierarchical levels using two hypotheses: (1) considering the three species as separate groups (Armeria The distance matrix was then used to produce a tree using the Neighbour-Joining algorithm.Branch support was also assessed through 1000 bootstrap replicates in PAUP* 4.0b10.Gene diversity (Hpop) was estimated for each population from the frequency of RAPD bands using Nei's (1973) algorithm adopted for RAPD by Peever & Milgroom (1994).Genetic distances between populations using Nei's genetic distance (1972) with an estimation of support for each cluster, were also inferred using TFPGA (Miller, 1997).
To further identify the genetic structure underlying our dataset, a Bayesian non-hierarchical clustering method was applied using STRUCTURE version 2.3 (Hubisz & al., 2009), which calculates the likelihood for the individuals being assigned to a given number (K) of groups based on minimizing linkage between groups, and maximizing linkage within them.Ten replicates of each value of K (1 to 15) were run under the correlated allele model (Falush & al., 2007) with a MCMC and burn-in settings of 10 6 and 10 5 iterations, respectively.We followed two criteria to detect the optimal number of K groups: First, the method suggested in the original STRUCTURE paper, which consists of comparing mean log likelihoods from the different cluster assignments penalized by one-half of their variance (Pritchard & al., 2000); Second, the highest increment in DK values (Evanno & al., 2005) from the different cluster assignments.
The correlation analysis between genetic diversity and altitude was performed using a Pearson product-moment correlation with SPSS 12.0.1 statistical software package (SPSS, Chicago, IL)

Morphometric analysis
A total of 30 characters were measured (Table 3) for 212 individuals, including at least one individual from each population surveyed in the RAPD study.Characters were selected based on previous studies (Nieto Feliner & al., 1996, 2001), and included sixteen quantitative continuous, one quantitative discrete, seven qualitative discrete (2 or 3-state), and five ratios of quantitative continuous characters (Table 3).For quantitative continuous characters three measurements were taken in each individual and the mean value was calculated.Standardization of variables was performed before the analyses.To display relative positions of individual phenotypes in the morpho-space, a Principal Component Analysis (PCA) was performed as an ordination method.A minimum spanning tree (MST) based on Euclidean distance was superimposed on the scatter-plot to help identify the closest relatives among different swarms.All these analyses were carried out with NTSYS-pc 2.11a.

RAPD variation and genetic relationships among taxa
A total of 62 RAPD bands, ranging between 200 and 1200 bp were scored, with an average of 7.7 bands per primer.After processing, no replicated RAPD profiles were found, so the number of RAPD multiband phenotypes was equal to that of individuals (N = 197).Only two primers generated products that were useful to distinguish taxa.Primer OPC-06 produced a 950 bp band that discriminated Armeria villosa subsp.bernisii from the remaining taxa.Two more bands (760 and 360 bp) generated by OPC-15 grouped, exclusively, individuals of A. villosa together with A. filicaulis subsp.trevenqueana.No diagnostic markers suitable for the discrimination of any of the three species were identified.A marker was considered diagnostic if it was present at a frequency equal to or greater than 90% in all individuals of a given taxon and present at a frequency lower than or equal to 10% in all individuals of other taxa.We considered this relaxed criterion to be justified due to the frequent occurrence of hybridization within this genus (Fuertes Aguilar & Nieto Feliner, 2003).In the PCO analysis, the first three coordinates explained 20.8% of the variance (8.6%, 7.0% and 5.2% respectively).Phenotypes spread across the space defined by the first two axes, which together accounted for 15.6% of the total    variance (Fig. 2).No clear-cut groups matched any of the species.In the scatter-plot, three clusters could be distinguished on the basis of taxonomic ascription.One large cluster encompassed A. villosa samples (subsp.longiaristata and subsp.bernisii), a second included all A. filicaulis representatives except those of A. filicaulis subsp.nevadensis, and a third grouped representatives of A. splendens and A. filicaulis subsp.nevadensis.Individuals of A. villosa subsp.bernisii and A. filicaulis subsp.alfacarensis connected the first and second groups.

Morphological variation
Since morphological characters are the basis of the present taxonomy of the studied taxa, it is not surprising that morphology was in accordance with the taxonomic ascription.In the PCA analysis, the first three Diversity of Armeria in Sierra Nevada components explained 62.4% of total variance (43.3%, 11.6, and 7.6%, respectively).The projections of individuals onto the morpho-space created by the first three axes (Fig. 4) displayed four main groups.The first three (A, B, C) each corresponded to one of the species occurring in the Sierra Nevada: Armeria villosa, A. filicaulis and A. splendens, while the fourth (D) was composed of A. filicaulis subsp.nevadensis individuals, which were connected to the A. splendens and A. filicaulis groups in the MST.Individuals of A. villosa subsp.bernisii clustered with those of A. villosa subsp.longiaristata, thus supporting the coherence of the species.The position of A. villosa subsp.bernisii was intermediate between A. villosa subsp.longiaristata and A. filicaulis subsp.filicaulis.Three outlier individuals of A. filicaulis (295A, 296A, 297A) from the Alfacar population (Fil7) were placed closer to A. villosa subsp.longiaristata than to their remaining conspecific individuals, especially when PC1 and PC2 were considered.The plants from this locality and its surroundings have recently received taxonomic recognition as A. filicaulis subsp.alfacarensis  3).A minimum-spanning tree based on the Euclidean distance helps to identify connection between groups.Centroids identify clusters for major taxa in the morphospace.Symbols as in Fig. 1. (Gutiérrez Larena & al., 2004).The UPGMA-tree based on Euclidean distance (not shown) distinguished the same groups outlined in PCA.

Population structure and genetic variation
Under both AMOVAs, groups by species (3) and groups by taxa (7), variation within populations was exceedingly higher than variation among populations and among groups.The percentage of total variance explained by among-populations differences decreased (from 25.1% to 17.9%) and the among-group differences increased (from 8.75% to 16.77%) when groups were considered as taxa (Table 4).The values of genetic diversity are presented in Table 5.The highest values in genetic diversity in the Sierra Nevada were found in Armeria villosa subsp.bernisii from mid to low-elevation populations (Ber5 and Ber15 with 0.30).In contrast, the lowest values were found at the extremes of altitudinal range: at very low elevation in A. filicaulis (Fil3) from Puerto Blanquillo (900 m) and at high-elevation populations (2870 m) of A. splendens (Spl31) (0.16 and 0.15, respectively).Values of Hpop ranged from 0.17 to 0.21 in A. filicaulis subsp.nevadensis to 0.16-0.25 in the remaining A. filicaulis populations.In A. villosa, the values for subsp longiaristata (0.21-0.25) fell within the range detected for subsp.bernisii (0.17-0.29).The Bayesian analysis contributed to detect a genetic structure with clusters independent of their taxonomic ascription (Fig. 5).The similarity coefficients calculated from the STRUCTURE results including all individuals showed that 1 and 2 are the only values of K that gave an unambiguous distribution of the individuals.In addition, the highest likelihood was obtained for K = 2 and K = 6, and the DK values also J. Fuertes Aguilar & al. supported K = 6.The six identified groups matched the population clusters and subclusters identified by the ΦST -based NJ-tree of populations.Actually, the same geographic pattern and a moderate congruence with taxonomy could be observed in the neighbour-joining tree based on ΦST pairwise values among populations (Fig. 5).The first cluster (C1, Fig. 5) included populations of A. filicaulis (Fil1,Fil2,Fil3,Fil4,Filtre10) and  A . villosa subsp. bernisii (Ber11,Ber9) from the Tejeda-Almijara range and Western Sierra Nevada, plus Fil7, a population from the dolomites of Alfacar that was genetically linked to another dolomiticolous population (Ber8) from the Cerro Trevenque.The second cluster was divided into two subclusters, the first one (C2, Fig. 5) exclusively composed of A. villosa, and the second one (C3, Fig. 5) encompassing all the populations of A. splendens and A. filicaulis subsp.nevadensis plus two of A. villosa subsp.bernisii (Ber12 and Ber26).The remaining group of A. villosa subsp.bernisii (populations Ber14, Ber 15, and Ber16), all from the southern slopes of the Sierra Nevada, fell apart from the others, being closer to the lowland A. villosa subsp.longiaristata (C4, Fig. 5).

Genetic variation and population structure
It is generally accepted that the parameters describing the genetic variation, its partition among and within populations and the genetic divergence among populations, result from the combination of two types of factors: life history traits (e.g.breeding system, seed dispersal, life form, geographic range) of the taxon and the evolutionary history of its studied populations.Based on data in various studies using molecular markers, Nybom & Bartish (2000) and Nybom (2004) estimated the average values of withinpopulation diversity and among-population diversity for several life history traits, thus: Hpop = 0.25 and ΦST = 0.25 for long-lived perennials; Hpop = 0.20 and ΦST = 0.26 for endemics; Hpop = 0.28 and ΦST = 0.34 for narrowly distributed plants; and Hpop = 0.27 and ΦST = 0.27 for obligate outcrossers.In contrast to a study on continental populations of another species of Armeria, A. maritima, in which Baumbach & Hellwig (2003) reported low values of genetic variability per population (0.122-0.193) and high levels of differentiation (ΦST = 0.46) for this perennial outcrosser, the values reported here for part of our studied Armeria species are in accordance with the parameters indicated by Nybom & Bartish (2000) and Nybom (2004).The average values of genetic diversity observed in A. filicaulis (Hpop = 0.20, ΦST = 0.34) and A. villosa (Hpop = 0.24, ΦST = 0.37) match those of obligate outcrossers with comparable distribution areas.Indeed, the partitioning of the RAPD variation revealed the highest proportion of within-population variation (66.2 % and 65.4%), which is in accordance with the pattern in obligate outcrossers.An important implication from hierarchical analysis of variance is that differences among groups (either species or  subspecies) accounted for a small percentage of the total variance (Table 4), suggesting the importance of interspecific gene exchange.The lack of diagnostic markers represented by fixed private bands in totalgenome markers for the species and subspecies is good evidence that gene flow occurs across species barriers, as detected in Pistacia (Yüzbasioglu & al., 2008).However, some of our studied taxa showed deviations from the expected values for obligate outcrossers.Armeria splendens (Hpop = 0.17) and, to a lesser extent, A. filicaulis subsp.nevadensis (Hpop = 0.19) were closer to those of facultatively allogamous plants than to the values expected for obligate outcrossers with a self-incompatibility system, like most Armeria species (Baker, 1966).However, no evident signs of breakdown in the incompatibility system were observed in these species since the two morphs associated to the incompatibility system are found (Fuertes Aguilar & Nieto Feliner, unpubl.data).The possible occurrence of polyploids in A. filicaulis subsp.nevadensis and A. splendens, which could break the diallelic incompatibility system, was also discarded after a cytological survey of high-alti-tude populations of the two taxa (Fuertes Aguilar & Nieto Feliner, unpubl. data).Rather, we think that the low genetic diversity values observed in populations of A. splendens and A. filicaulis subsp.nevadensis, and especially in the case of the former, are not related to life history traits but to the recent history of their habitats.The small-sized scattered populations of A. splendens strictly depend on habitats surrounding lakes created by morraines after deglaciation.Present populations of A. splendens were therefore in the frontline of re-colonization of deglaciated environments after the last stadial, and might have suffered a loss of genetic variation due to founder effects and bottlenecks during migration (Gabrielsen & al., 1997;Piñeiro & al., 2011).In the case of A.filicaulis, the populations of subsp.ne vadensis, concentrated in N-NW slopes of the Sierra Nevada, are altitudinally and geographically isolated from the remaining populations of A. filicaulis, and therefore their inter-populational gene flow might be restricted and so favouring genetic drift and stronger genetic differentiation, as reported in other alpine taxa by

Genetic diversity along an altitudinal gradient
Correlation analyses between populational parameters and geographic distances have been used to study the history of postglacial migrations (Ga brielsen & al., 1997, Holderegger & al., 2002).In the present study, the linear distances are so small, the topography is so irregular, and the groups are so heterogeneous that any kind of corrected geographic distance would be inaccurate to describe true among-population relationships.Nevertheless, the distribution of the populations along an elevation gradient perpendicular to the direction of the glaciation/de glaciation front allows us to explore whether recolonization influenced any of the population parameters.The correlation analysis (Fig. 6) produced a significant negative correlation (r 2 = 0.410, d.f.= 18, p = 0.003) between elevation and population genetic diversity (Hpop) for Armeria villosa subsp.bernisii.The populations of A. splendens and A. filicaulis subsp.nevadensis analysed together (to compensate their low number) gave no significant result (r 2 = 0.497, d.f.= 7, p = 0.051), probably due to the low number of populations sampled, but they showed a trend towards a lower Hpop values as the altitude increased.To perform the analysis on A. villosa subsp.bernisii, a correction for elevation values was needed for populations on the southern slope of the Sierra Nevada.Climatic and vegetation conditions for Supra-Medi terranean and Oro-Mediterranean thermotypes on the southern slopes appear approximately 200 m higher than on the northern slopes (Lorite & Valle, 1999;Rivas-Martínez & al., 1986).Therefore, for these specific correlation analyses, we compensated the elevation values by decreasing the respective altitudes by 200 m for the nine populations of A. villosa.subsp.bernisii from the southern slope (Ber12 to Ber20, Table 5).The continuous decrease of Hpop for such populations as the elevation increases is clear: populations below 1850 m were in the 50% percentile, whereas the lowest values were found above that elevation.No other common factor groups these populations with low Hpop except for the population of A. villosa subsp.bernisii from Gador (Ber18), where the isolation from neighbouring populations (15 km to the closest one) and habitat degradation might be determinant.The fact that the lowest values for A. villosa.subsp.bernisii are similar to those exhibited by A. splendens and A. filicaulis subsp.nevadensis suggests that glaciation events could have influenced the genetic depauperation of those populations.However, the trend observed in A. villosa subsp.bernisii suddenly ends for the high-altitude populations (Ber13, Ber20, Ber21, Ber25), which depart from the common trend, showing values of Hpop 25-50% higher than those of populations immediately below in altitude.It is also noteworthy that the altitudinal belt between 1400 and 1900 m, where A. villosa.subsp.bernisii reaches its highest diversity values is the only gap in the elevational distribution of A. filicaulis.A possible explanation for this pattern is that A. villosa subsp.bernisii genetically assimilated previously existing populations of A. filicaulis subsp.nevadensis in this range.
We suggest that the pattern of genetic diversity along the altitudinal gradient in the Sierra Nevada is congruent with that presented in previous evidence based on plastid (Gutiérrez Larena & al., 2002) and ITS sequences (Nieto Feliner & al., 2002), and could be explained by two parallel waves of altitudinal colonization into deglaciated habitats, one by Armeria fili caulis subsp.nevadensis plus A. splendens and another by A. villosa subsp.bernisii and these are reflected by the two different regression slopes (Fig. 6).The higher than expected values in high-altitude populations of A. villosa.subsp.bernisii may be due either to the presence of populations of this taxon prior to glaciation (nunatak hypothesis) or, alternatively, to recent gene flow between A. filicaulis.subsp.nevadensis and A. splendens that created an additive effect in border populations produced via hybridization.Another possible indication of the genetic assimilation of other congeners' populations by A. villosa subsp.bernisii is the observed increase of genetic diversity in high-altitude populations (Ber12, Ber13, Ber26) (Fig. 6), two of which cluster with A. filicaulis subsp.nevadensis in the NJ-tree (Fig. 5).The process of assimilation of taxa by genetically "aggressive" species has been previously reported for A. villosa (Fuertes Aguilar & al., 1999) and it conforms to the compilospecies concept (Harlan & De Wet, 1963).A similar evolutionary scenario has been reported to occur during overlapping postglacial migration of different taxa in Arctic areas (Brochmann & al., 1996).

Genetic distance and morphological cohesion
The NJ-tree produced by analysing individual RAPD phenotypes shows a pattern analogous to that obtained with sequence data (Gutiérrez Larena & al., 2002;Nieto Feliner & al., 2004): few (6) recognizable groups that are taxonomically heterogeneous, mostly with a geographical and geological common ground, and with low BS (Fig. 3).The situation changes dramatically when the analysis is made on a populational basis, using estimates of ΦST as genetic distance (Fig. 5).The three main clusters detected are largely congruent with the morphological variation and taxonomy, except for those taxa with a putative hybrid origin and a few idiosyncratic populations, already detected through morphological and molecular incongruence (Fil7, Filtre10, Ber8, Ber9, Ber11).Cluster 1 gathers populations from the Tejeda-Almijara range and the nearby Western Sierra Nevada, thus suggesting connections through migration and genetic exchange between these areas (Gutiérrez Larena & al., 2006).The large cluster 2 encompasses almost exclusively populations of Armeria villosa, except for the recently described populations of A. filicaulis subsp.alfacarensis (Gutiérrez Larena & al., 2004).Cluster 3 groups A. splendens with populations of hybridogenous A. filicaulis subsp.nevadensis.This cluster also includes two populations of A. villosa subsp.bernisii (1950 m) that are sympatric to those of A. filicaulis subsp.nevadensis.This connection lends support to the hypothesis that A. villosa.subsp bernisii originated through genetic swamping of A. filicaulis populations with A. villosa invaders from outside the Sierra Nevada (Guitérrez-Larena & al., 2002).
An unexpected outcome obtained both from RAPD and morphological analyses is that populations that occur on two rather isolated crystalline dolomitic enclaves (Alfacar, Cerro Trevenque) are genetically closer than previously thought.Armeria filicaulis subsp.alfacarensis is connected in the MST to individuals of A. villosa subsp bernisii from Trevenque and all of them are placed in proximity to A. filicaulis subsp.trevenqueana in the Principal Component Analysis of morphometric data.The question arises whether there is a common origin for those dolomiticolous taxa, or whether the similarity is a result of introgression among them.A third alternative is a strong selective pressure for morphological and genetic characters adapted to such a demanding environment.Although there is a recent study supporting the relevance of substrate type as responsible for alpine plant distribution patterns the answer to these questions requires further research.

Glacial dynamics and paleoecological data
Paleoecological reconstructions for the Sierra Nevada environments during Pleistocene are still very incomplete.A significant difference in the effects of glaciations on Mediterranean massifs such as the Sierra Nevada compared to those at higher latitudes such as the Pyrenees and the Alps is that glaciers were smaller and much more adapted to relief in the former massif.As a consequence, water availability from melting ice was more geographically limited.This feature, combined with the reduced water available from precipitation during glacial periods, promoted the predominance of arid (steppic) ecosystems in areas only moderately distant from melting-ice water.The result was an enforcement of the island effect on some alpine plants.Another effect of climatic changes during glacial maxima was the contraction of bioclimatic belts.Given that the Sierra Nevada relief has remained more or less constant since the onset of Pleistocene glaciations, we can infer that conditions for the crioromediterranean bioclimatic belt descended to c. 1800 m on the north western slope of the massif (Gómez Ortiz & Salvador i Franch, 1996).This contraction caused a more marked ecological gradient between xeric and humid conditions along the slopes (Ferris & al., 1999).The sequence of changes has been recorded at the Padul site (950 m) for the Sierra Nevada, with the appearance of steppic vegetation (Artemisia, Ephedra, Chenopodiaceae) in basal Sierra Nevada zones during Last Glacial Maximum (Pons & Reille, 1988).There is also the possibility that high-altitude plants survived during the glaciations in nunataks, i.e. areas that remained ice-free, but the existence of nunatak survival in European mountain systems and its relevance for the interglacial recolonization process is still under discussion (Stehlik, 2000(Stehlik, , 2003;;Stehlik & al., 2002).The arrival of the genus Armeria in the southern Iberian Peninsula is difficult to date and paleopalynological data contribute little due to the low pollen production in the genus (Barrón & al., 2010).There is no evidence of Armeria fossil pollen recorded in the Padul sediments (Pons & Reille, 1988) but a sample from a hyena coprolite from Late glacial (12700 yr BP) c. 20 km north of the Sierra Nevada may belong to Armeria, possibly A. villosa (Carrión & al., 2001b;Carrión, unpubl. data.).This sample was originally assigned to Limonium, a genus with a pollen morphology indistinguishable from that of Armeria.However, Limonium is always associated with alkaline soils and never with the Quercus aff.pyrenaica forest that the paleoecological reconstructions indicate for this deposit.
Two previous studies based on plastid DNA and ITS markers proposed a scenario of colonization of the Sierra Nevada by successive migrations (Gutiérrez Larena & al., 2002;Nieto Feliner & al., 2004).On this view, Armeria splendens was the first species that colo -nized the Sierra Nevada. A. splendens belongs to a group of high-elevation species including A. alpina from the Alps and the Pyrenees, and A. bigerrensis and A. caespitosa from the Iberian Central System, with which it shares the same plastid haplotype (Gutiérrez Larena & al., 2002, Fuertes Aguilar & Nieto Feliner, unpubl. data).Of the three species present in the Sierra Nevada, A. splendens is also the only one that has exclusively the endemic ITS ribotype R3 (Nie to Feliner & al., 2004).In a subsequent phase, A. splendens would have originated A. filicaulis subsp.nevadensis by hybridization with A. filicaulis.This scenario is consistent with the finding that, based on RAPDs at the population level, A. splendens has its lowest genetic distance with A. filicaulis subsp.nevadensis populations (Nev27, Nev28, Nev29, and Nev30).The genetic diversity values of Armeria filicaulis subsp.nevadensis and A. splendens also decrease with increasing altitude, an effect observed in alpine plants during the colonization of glacier forelands along the successional series following ice retreat (Raffl & al., 2006).Another possible explanation for the low genetic diversity of A. splendens and A. filicaulis, as compared with A. villosa subsp.bernisii, is population size (Fig. 6), which typically is positively and significantly correlated with genetic diversity (Nybom, 2004).Armeria splendens and A. filicaulis have smaller population sizes (70-1000 individuals) than A. villosa subsp.bernisii, with populations of more than1000 individuals (Blanca & al., 2002).Other factors responsible for genetic depauperations and local extinction, such as genetic bottlenecks caused by brief periods of glacial activity, have also been recorded in the Late Glacial and Holocene for taxa in SE Spain (Carrión & al., 2001a;Carrión, 2002).
What kind of ecological conditions could have favoured the contact between Armeria splendens and A. filicaulis to produce the hybrid taxon A. filicaulis subsp.nevadensis?Assuming that both taxa have maintained their current ecological requirements, a descent in the altitudinal range of A. splendens so as to converge with that of A. filicaulis in W and NW Sierra Nevada may have occurred during the maximum contraction of vegetation belts in the Sierra Nevada.Such conditions may have only been met during the maximum glacial stages.The same kind of predominantly periglacial xeric conditions have been recorded in the Cañada de la Cruz (Sierra de Cazorla) series by Carrión & al. (2001a).Assuming niche conservatism, at the end of that period and during the present interglacial stage, A. splendens, and to a lesser extent, A. filicaulis subsp.nevadensis, could have extended their range upwards following their optima, until they reached their present distribution.In contrast to the dry and cold environmental conditions associated with the establishment of A. splendens and A. filicaulis in the Sierra Nevada, the climatic conditions for the wave of colonization in the massif by A. villosa are likely to have been considerably wetter and warmer.This species grows in the understory of Pinus and Quercus spp.forests, whose ranges expanded in Betic and Sub-betic areas during interglacial periods (Carrión & al., 2001a).The concomitant forest expansion across the massif probably facilitated the hybridization between A. villosa and A. filicaulis, resulting in the partial assimilation of the latter to create what is now known as A. villosa subsp.bernisii, and which is currently distributed over a large part of the forested slopes of the Sierra Nevada.This scenario is consistent with the clustering of two populations of the latter taxon (Ber12, Ber26) with those of A. splendens and A. filicaulis subsp.nevadensis in the NJ-tree based on genetic distances among populations (Fig. 5), as well as with the rare occurrence of the plastid haplotype L that is predominant in both A. villosa subsp.bernisii and A. filicaulis subsp.nevadensis (Gutiérrez Larena & al., 2002).A possible trace of the processes described above, which involved repeated genetic exchange between the Armeria congeners, is the predominance of haplotype I (from A. splendens) in populations of A. villosa subsp.bernisii and the occurrence of haplotype E, which is frequent in A. filicaulis but absent in populations of A. villosa.subsp.bernisii from the Sierra Nevada except at El Dornajo (Ber25) (Gutiérrez Larena & al., 2002).The expansion undergone by this latter taxon was probably the result of population fragmentation and long-distance dispersal, as reflected by the absence of a congruent geographic pattern among haplotype and ribotype distributions in the Sierra Nevada, and in the genetic distance tree (Fig. 5).A relevant part of this process of expansion has involved the migration of A. villosa.subsp.bernisii outside the Sierra Nevada and the ability to colonize the Tejeda-Almijara range (Gutiérrez Larena & al., 2006).
In summary, this work shows how a populational approach encompassing different species may help to unravel the evolutionary mechanisms (hybridization, genetic drift, gene flow) that were involved in the shaping of the present taxonomic and genetic diversity of the genus Armeria in the Sierra Nevada.In glacial refugia, such complex interactions between different taxa are often difficult to interpret using purely phylogeographic approaches.However, our present understanding of population genetics in arctic and alpine postglacial colonization, in conjunction with the infe-rence of ecological conditions during glacial/interglacial cycles, provides an explanation to how morphological and genetic diversity was moulded in southern European refugia.

Fig. 2 .
Fig. 2. Genetic variation among 197 studied individuals of Armeria spp. as displayed by a bi-dimensional scatter-plot using the two first coordinates of a principal coordinates analysis based on a Dice similarity matrix of RAPD's phenotypes.Symbols as in Fig. 1.

Fig. 3 .
Fig. 3. Neighbour-joining tree based on distance transformed from Dice similarity matrix from all analysed Armeria spp.individuals.Taxonomic ascription and individuals identified as in Table2.

Fig. 4 .
Fig. 4. Morphological variation among 212 studied Armeria spp.individuals, as displayed by tri-dimensional scatter-plot using, the three first components of a principal components analysis based on 30 morphological characters (Table3).A minimum-spanning tree based on the Euclidean distance helps to identify connection between groups.Centroids identify clusters for major taxa in the morphospace.Symbols as in Fig.1.

Fig. 5 .
Fig. 5. Left, Neighbour-joining tree of studied Armeria spp.populations using ΦST as an estimation of genetic distance showing four main clusters.Symbols as in Fig. 1.Right, genetic structuring of Armeria spp.individuals detected by genetic admixture analysis in STRUCTURE for K = 6.
Taxonomic synopsis of Armeria taxa involved in this study.

Table 2 .
List of populations and individuals sampled for RAPD and morphometric studies of Armeria spp.(Continuation).

Table 2 .
List of populations and individuals sampled for RAPD and morphometric studies of Armeria spp.(Continuation).

Table 2 .
List of populations and individuals sampled for RAPD and morphometric studies of Armeria spp.(Continuation).
A. villosa subsp.longiaristataJaén:La Matea, road to El Cerezo, 30SWH3515, 1280 m, roadside.262LA.villosa subsp.longiaristata Jaén: Cazorla, Sierra de Cazorla, road from la Nava de San Pedro to El Cerezo, a Missing codes under Population column correspond to individuals exclusively used in in the morphometric study.

Table 3 .
Morphometric characters used in the multivariate analysis of Armeria spp.individuals.

Table 5 .
Synthesis of geographical, genetic and morphological information of studied populations of Armeria spp. in southern Spain.Populations codes correspond to those indicated in Table1.