Microsatellite Marker Analysis Reveals the Complex Phylogeographic History of Rhododendron ferrugineum (Ericaceae) in the Pyrenees

Genetic variation within plant species is determined by a number of factors such as reproductive mode, breeding system, life history traits and climatic events. In alpine regions, plants experience heterogenic abiotic conditions that influence the population's genetic structure. The aim of this study was to investigate the genetic structure and phylogeographic history of the subalpine shrub Rhododendron ferrugineum across the Pyrenees and the links between the populations in the Pyrenees, the Alps and Jura Mountains. We used 27 microsatellite markers to genotype 645 samples from 29 Pyrenean populations, three from the Alps and one from the Jura Mountains. These data were used to estimate population genetics statistics such as allelic richness, observed heterozygosity, expected heterozygosity, fixation index, inbreeding coefficient and number of migrants. Genetic diversity was found to be higher in the Alps than in the Pyrenees suggesting colonization waves from the Alps to the Pyrenees. Two separate genetic lineages were found in both the Alps and Pyrenees, with a substructure of five genetic clusters in the Pyrenees where a loss of genetic diversity was noted. The strong differentiation among clusters is maintained by low gene flow across populations. Moreover, some populations showed higher genetic diversity than others and presented rare alleles that may indicate the presence of alpine refugia. Two lineages of R. ferrugineum have colonized the Pyrenees from the Alps. Then, during glaciation events R. ferrugineum survived in the Pyrenees in different refugia such as lowland refugia at the eastern part of the chain and nunataks at high elevations leading to a clustered genetic pattern.


Introduction
The maintenance of genetic diversity within populations is essential for the conservation of species since it can enable adaptation to future environmental changes [1], [2]. With high genetic diversity, it is more likely that some individuals will possess alleles that could allow the population to adapt to a new habitat [3]. High genetic diversity may be maintained via gene flow among individuals and among populations. In recent years, many studies have shown strong impacts of global change on species' genetic diversity [4]. Because climate warming is greater at high than at low elevations [5] alpine biodiversity hotspots that are home to many endemic and rare species may be particularly threatened by climate change [6]. Thus, to predict how these biodiversity hotspots will deal with climate changes, it is important to understand the genetic patterns of the species and their changes over time [7].
Genetic structure and genetic divergence within and among plant populations result from a number of contemporary and historical factors acting at various temporal and spatial scales. Alpine habitats are typically characterized by small spatial changes in environmental variables such as topography, temperature, snow cover, soil moisture or bedrock [8] favorable to phenotypic differentiation and local genetic variability over short geographical distances. In addition to abiotic factors, life history traits can also influence genetic features of populations. For instance, selfing perennials and clonal species generally exhibit lower levels of genetic diversity and higher levels of population differentiation than outcrossing perennials [9]. In contrast, mixed mating systems maintain a high level of genetic diversity and may be advantageous for plants in variable alpine environments [6]. The combined effects of both life history traits and small-scale abiotic heterogeneity may result in local genetic adaptations and high genetic differentiation of populations, notably where gene exchanges between populations are restricted as is often the case in mountain habitats [10]. Therefore, plant-community fragmentation pattern, steep climatic gradients leading to asynchronous flowering phenology or disjoint pollinator communities [7], and topographical barriers may hamper gene flow (through seed and pollen movements) and may increase differentiation between populations.
At large scales, species genetic patterns are often shaped by past climate-driven range dynamics [11]. In temperate regions, species ranges have changed dramatically over the Pleistocene with its successions of glaciations and warming periods. In mountains, periodic advances and retreats of glaciers have, on several occasions, forced species to descend to lowland and peripheral refugia and/or to survive in isolated ice-free areas above glaciers (nunataks). Genetic isolation, bottlenecks and potential allopatric evolution of populations in refugia have been followed by repeated founding events during population expansions, which have led to loss of genetic diversity and heterozygosity [11]. Additionally, populations located at eroding margins of species distribution ranges following climate warming could have persisted in suitable isolated lowland habitats [12]. It has been suggested [13][14] that such small populations isolated for so long could display low within-population genetic diversity (as a result of genetic drift) but have disproportionally high levels of between-population differentiation and regional genetic diversity and distinctiveness [15].
The genetic pattern and phylogeography of many species [16], [17] have been investigated in various mountainous regions, but to date, only a very limited number of studies have focused on the Pyrenees [18][19][20]. The Pyrenees are a mountain range stretching across the isthmus that lies between the Iberian Peninsula and the rest of the European continent. Acting as a biogeographical barrier during postglacial expansions, they are a limit for numerous Iberian and northern species. However, lower lands at the eastern and western extremities of the mountain range have been corridors for migration of several plant and animal species [21,22]. Due to oceanic influences to the west and Mediterranean influences to the east of the chain, the climate varies greatly along a longitudinal gradient in the Pyrenees (Météo France data). This bioclimatic pattern potentially allowed some species to survive glaciations in various refugia along the chain. Studies have suggested that species have survived glaciations [23][24][25] in two main refugia -one in south central and another in the eastern Pyrenees -thanks to mild Mediterranean and oceanic climatic influences. However both geological and biological data [26] demonstrate that during glacial periods the ice sheet did not completely cover the highlands and rarely descended below 1000 m a.s.l except for big glaciers pushed down into large valleys and to the plains [23]. These nunataks and unglaciated peripheral area could have provided refugia for some species [18] and determined several recolonizing routes sculpturing a specific and complex genetic structure in the Pyrenees.
Here we investigate the genetic structure and phylogeographical history of the alpine shrub Rhododendron ferrugineum (Ericaceae) in the Pyrenees. This species is insect pollinated and reproduces both sexually (mixed mating system) and vegetatively through layering. It is widely distributed at high elevation, on non-calcareous bedrock. The aim of the present study is to understand to what extent the phylogeographical history of the species has shaped its present genetic patterns. The Asian origin of the genus Rhododendron [27] implies an east/west migration from Asia to Europe. In France, this would lead to the Alps being colonized by R. ferrugineum before the Pyrenees, with the Pyrenean populations most likely derived from populations from the Alps. After the migration to the Pyrenees, the species had to survive successive glaciation events. Different hypotheses of glacial refugia and postglacial colonization and their consequences on the genetic patterns of R. ferrugineum may be proposed: 1. The species could have survived glaciation episodes in one or both of the two main refugia at the eastern and south central part of the Pyrenees (HYP 1; [24]). Thus, the highest genetic diversity should be found in one or both of these areas and the genetic pattern would be shaped by founder effects i.e. gradual reduction of genomic variability along the colonizing routes from main refugia. This would have occurred because colonization usually involves only a fraction of the genetic diversity present in refugial areas [28]. 2. The species could also have survived in nunataks along the mountain chain (HYP 2). Hypothetically, high elevation nunatak populations would have experienced long periods of isolation, inbreeding and genetic drift [18]. This pattern would lead to a patchy distribution of genetic clusters with low levels of genetic diversity and high population differentiation; 3. A pattern with less differentiated clusters than in the previous hypothesis may be observed if the species had survived in large lowland refugia in which extensive gene flow had occurred. The recolonization of highland areas from lowland refugia could have left small, isolated and marginal rear-edge populations like those currently observed at low altitude along the chain (HYP 3). Indeed, the strong abiotic mountain heterogeneity may have allowed remnant populations to match suitable ecological conditions at intermediate or low altitudes. Moreover, the long life span of the clones [29,30] could have permitted their survival even with low seedling recruitment. Such rear-edge populations should exhibit more rare alleles and higher among-population genetic diversity than populations recently established from high altitude source populations [15]. These three hypotheses are not exclusive and it is conceivable that the genetic pattern observed may result from the different mechanisms at work successively or simultaneously within the bioclimatic regions of the Pyrenees inducing diverse phylogenetic patterns along the chain.

Species studied
Sampling of the species of interest was conducted in accordance to national and international guidelines (authorization from the ONF for the populations BetL, BetH, LecL, PraL, PuiH, from the Parc National des Pyrénées for the population TroH, from the Conseil général des Pyrénées-Orientales for population BouH, and no permissions were needed for the other sites).
Rhododendron ferrugineum L. (Ericaceae) is an evergreen shrub with a mean height of 70 cm that dominates subalpine landscapes in the Alps and Pyrenees between 1500 to 2200 m a.s.l, on north-towest facing slopes, sometimes reaching 90-100% of the vegetation cover. It is also found in the Jura Mountains, in the Northern Apennines and in the Apuan Alps [31]. In the Pyrenees, infrequent smaller and isolated populations can be found lower down, between 900-1500 m a.s.l. A single population of some scattered individuals remains at low altitude (1500 m a.s.l.) in the French Jura Mountains. Rhododendron ferrugineum is a clonal, selfcompatible [32] mass-flowering species (producing until 3000 flowers per m 2 ; [33]) pollinated by honeybees and bumblebees [32,34].

Study sites and sample collection
We sampled 29 populations across the Pyrenees, the population in the Jura Mountains and three populations in the Alps ( Figure 1, Table 1) ranging from 1070 m to 2200 m a.s.l. The populations were divided into three groups according to their elevation: low altitude (hereafter L; from 1000 to 1450 m a.s.l), intermediate altitude (hereafter I; from 1450 to 1900 m a.s.l.), high altitude (hereafter H; from 1900 to 2200 m a.s.l.). In the Pyrenees, populations were sampled along a longitudinal transect through three regions: Languedoc-Roussillon (eastern part), Midi-Pyrénées (central part) and Aquitaine (western part). Due to calcareous bedrock, populations in the western part of the Pyrenees are rare and we only sampled one. An estimation of population size is given in Table 1 [,1 ha, 1-10 ha, 10-20 ha, 20-50 ha, 50-100 ha, .100 ha). The mean geographical distance between populations from Pyrenees vs. Alps, Pyrenees vs Jura Mountains and Alps vs. Jura Mountains was 900 km (ranging from 515 to 1430 km), 695 km (from 590 to 830 km) and 335 km (from 165 to 640 km) respectively. The distances between populations across the Pyrenees ranged from 3 to 370 km with an average of 126 km. In each of the 33 populations (except for the very small JuraL and LapL populations), young leaves of twenty individuals were sampled (a total of 645 individuals sampled), separated by at least 5 m in order to avoid sampling clonal individuals, were dried in silica gel in the field and then stored at 220uC.

DNA extraction and microsatellite markers
For each individual sampled, 20 mg of leaf tissue was ground with two tungsten beads (diameter 3 mm) in a sterile 2 ml Eppendorf tube at 30 Hz for 90 s using TissueLyser II (Qiagen). Then, the quality and amount of DNA obtained from the DNeasy-Plant Mini Kit (Qiagen, Courtaboeuf, France) following the manufacturer's instructions was verified using a Nanodrop ND 1000 spectrophotometer (Peqlab, Erlangen, Germany).
We used nuclear microsatellites as suitable molecular markers to reconstruct regional phylogeographical patterns [18]. Chloroplasts (personal observation) or other types of genetic markers such as AFLP proved to be little informative in this species [35]. The 27 polymorphic microsatellite markers (Table S1) were previously developed [36,37] using pyrosequencing technologies (454 FLX Titanium, Roche Applied Science, Meylan, France). Forward primers were labeled with fluorochromes (Eurofins MGW Operon, Courtaboeuf, France) and were used in six multiplexed PCRs, optimized using Multiplex Manager 1.1 [38].
Amplifications were carried out in a 10 ml mix containing 10 ng of template DNA, 0.7 6 Qiagen Multiplex PCR Master Mix, 10 mM of each primer and RNase free water. Cycling conditions were 15 min at 95uC, 30 6 (30 s at 95uC, 90 s at 56uC, 45 s at 72uC) and 30 min at 60uC. Each set of reactions included a negative (water) and a positive (known genotype) control. Fragment lengths were read on an ABI 3730 sequencer (Applied Biosystems, Courtaboeuf, France) and scored with Genemapper version 4.0 software with the GeneScan -600 LIZ (Applied Biosystems, Courtaboeuf, France) as internal size standard.

Data analysis
Genetic diversity. All 645 individuals were tested for clonality with GenAlEx 6.5 software [39]. The presence of null alleles at each locus was examined using Micro-Checker software [40], and tests for departures from Hardy-Weinberg equilibrium (HWE) were tested by Fisher's exact tests using Genepop 3.4 [41]. We used GenAlEx 6.5 software to calculate for each locus and population the number of alleles, the allelic richness (Ar), the observed (H O ) and expected heterozygosity (H E ) and the inbreeding coefficient (F IS ). Population structure and differentiation. The pairwise F ST (level of genetic differentiation; [42]) among populations were calculated with Genodive v2.0 [43]. Then, the statistical significance of the pairwise F ST values was assessed by 1000 random permutations of individuals across populations with Genetix 4.05 [44].
To analyze the genetic structure among populations, we used Structure software (ver. 2.3.3; [45]) in order to determine the most likely number of homogenous clusters in the sample (K). The admixture model was applied with 50 runs for each K value from 1 to ''number of populations +1'', each run comprising a burn-in period of 20 000 iterations followed by 100 000 iterations. The optimal value of K was evaluated by considering the highest mean likelihood value of K, i.e., L (K), as well as the DK method [46]. These results were computed and visualized in the online interface Structure Harvester [47]. To evaluate the distribution of genetic variance within and among populations, and among groups when populations were nested into genetic clusters, altitude (low, intermediate and high altitude) or population size (according to the six population size groups in Table 1) we performed analyses of molecular variance (AMOVA; [48]) using Genodive v2.0. The significance of the results was tested using 999 random permutations of the data.
Mantel tests [49] were performed (''ncf'' package implemented in R, R Development Core Team, 2008) to test for isolation-bydistance (IBD; divergence due to drift and mutation proportionally increasing with geographical distance), investigating the correlation between matrices of genetic (F ST ) and geographical (km) distances. Mantel test were also used to determine wheter elevational separation was related to genetic differentiation. Calculation of the gene flow on the basis of rare alleles was conducted [50] and the number of migrants (N m ) was corrected for the number of samples.

Results
A total of 153 alleles at 27 microsatellite loci were detected across 645 individuals at all sites (Table S2). No clones were identified among the individuals sampled as revealed by the test of clonality.
Significant departures from HWE across loci were detected in 12 out of 33 (or 36%) sites. The RF87P2 marker was the only locus that showed no departures from HWE. Most of the loci showed significant departures from HWE in two to seven populations, except RF14P3, RF47P1 (11 populations each), RF56P1 (13 populations) and RF96P2 (31 populations). These instances likely reflect occasional departures from random mating rather than the presence of null alleles. Indeed, the probability of null alleles was low for all loci, according to Oosterhout's method [40], mean estimation of null alleles was 20.031 (ranging from 20.3675 to 0.1352). No differences in the results were observed when including the microsatellite marker RF96P2 or not, so we decided to keep it in the analyses.

Genetic diversity and allelic richness among populations
The number of alleles ranged from 2 to 15 (mean 6.25) alleles per locus. A total of 16 private alleles (restricted to a single population; Table 1) was detected (10.45%), mostly in AlpH (10) and westernmost (EspI, PieI) Pyrenean populations. These latter populations were small and isolated (,1 ha) except the large and central population EspI (.100 ha). Significant heterozygosity excesses (F IS coefficients ranging from 20.07 to 20.21) were detected in three eastern (LecL, LapL, PuiH) and four central populations of the Pyrenees (CamI, VicH, VicL, TecL). In contrast, significant heterozygosity deficiency was found in four populations with particularly high values of F IS coefficients in SomI (0.16) and EspI (0.31). Because the frequencies of null alleles did not appear to be significant, the highest F IS values could be due to inbreeding.

Genetic differentiation and structuration among populations
The F ST values (Table S3) were significantly different from zero for all population pairs, except for three pairwise comparisons: BetH/BetL and SouH/SouL that were geographically really close (less than 3 km) and CamI/EynH (35 km apart). The mean F ST value was relatively low (0.267) and most of the highest values (.0.33) were found in the Pyrenees between eastern and western populations, reaching 0.45 for LecL/AraH. The population from the Jura Mountains showed a highly variable level of F ST (from 0.25 to 0.42).
To infer the relationship between populations we used the clustering algorithm implemented in the program Structure. Comparing all 33 populations (K: 1 to 34), DK as a function of K reached a peak at K = 2 and then reached a second peak at K = 6 (Figure 2a). At K = 2, populations from the Alps, Jura and the eastern Pyrenees were assigned to the same cluster while populations from western and central Pyrenees were clearly grouped apart. Individuals from AlpI showed a similar proportion of membership to both genetic clusters. At K = 6, populations from the Alps and Jura were assigned to one cluster (CL1) and those of the Pyrenees were grouped in 5 clusters (CL2 to CL6) along a longitudinal gradient. Except between CL5 and CL4 or CL6, there was low admixture between the other Pyrenean clusters.
The analysis of molecular variance (AMOVA; Table 2) based on F ST showed that as expected for polymorphic loci such as microsatellites, the within-population variance component was rather high (.66%, P = 0.001). The among variation component (F ST ) ranged from 9% to 31% according to the grouping and was statistically significant in all combinations tested (P = 0.001). At a superior hierarchical level, the genetic clusters could explain a significant part of the total genetic variance (24.4%, P = 0.001) whereas it was not significant when grouping populations into altitudinal groups (P = 0.956) or population size (P = 0.222).
Allelic richness (Ar) was highest in CL1 (Alps and Jura, 5.2) and lowest in CL6 (western Pyrenees, 2.9) and ranged from 3.2 to 3.8 in the other clusters (Table 3). F ST values were all significantly different among clusters ( Table 4). The highest F ST values were found between CL2 (eastern Pyrenees) and CL5 or CL6 (western Pyrenees, 0.29 and 0.31 respectively). The F ST values between CL1 and the Pyrenean clusters were quite low (ranging from 0.16 to 0.20).
There was a positive relationship (Mantel test) between genetic distance and geographical distance matrices across populations located in the Pyrenees (R = 0.784, P = 0.001). The increase of isolation with geographical distance was consistent with the previous results (see Structure analyses). When including populations from the Alps and Jura, the Mantel test was still positive (R = 0.26) but not significant (P = 0.383), confirming the lack of strong differentiation among the populations from the Alps and Pyrenees. At the Pyrenean scale, the correlation among elevational distances and FST was slightly positive but not significant (R = 0.092, P = 0.06).  [51]). Although clonality and selfing are known to decrease genetic variability [52], the long life span of individuals [29,30] and the strong life-time inbreeding depression eliminating most selfed individuals (personal data) may contribute to maintaining relatively high genetic diversity even in small R. ferrugineum populations. Accordingly, most of the populations studied showed a F IS close to zero and a great proportion of the genetic variation was found within populations as observed in other R. ferrugineum populations [31] and other alpine species [19], [53], [18]. However, selfing has probably contributed to individual recruitment in populations harboring higher heterozygote deficiency (SomI, PieI, EspI, SpoI). Such high heterozygote deficiency has also been reported in the Apennines and the Alps [31]. Surprisingly, significant heterozygote excess    Alpine origin of the two Pyrenean R. ferrugineum lineages

Gene flow among populations
The Rhododendron genus is known to originate from Asia [27] and colonized Europe from east to west. Our hypothesis was that migration from the Alps to the Pyrenees occurred with a loss of genetic diversity along the colonization routes. In our study, we identified two lineages in the Alps: the first, found only in the northern population (AlpI), was related to western Pyrenean populations. This lineage would possibly be found in more populations with more extensive sampling in the Alps. The second lineage, related to the eastern and central Pyrenees, was found in the same population AlpI but also in south and central-eastern populations from the Alps and in the Jura Mountains. At K = 6, populations from the Jura Mountains and the Alps are found in the same cluster, that may indicate the occurrence of recent gene flow between the Alps and the Jura Mountains, possibly in low land refugia during glaciations. It would be interesting to complete our data by a large sampling of populations from the Alps and the Jura Mountains in order to detail the European phylogenetic pattern of the species. Compared to the Alpine populations, Pyrenean populations were characterized by a substantial loss of genetic diversity (Ar, H E , H O ) suggesting bottlenecks during the migration events [54], [16]. Cold, dry steppe/tundra occupied most of the land between the Alps and the Pyrenees during the last maximum glaciation [55] and could have allowed migration of Alpine species such as R. ferrugineum either from the southwestern Alps [56], [57] or from northern Alpine populations through the French Massif Central which is considered to be a key area in the colonization of the Pyrenees by Alpine species [58].

Genetic pattern and phylogeographical history of R. ferrugineum in the Pyrenees
We found that genetic diversity of R. ferrugineum was not correlated to altitude or population size; instead there was a sharp longitudinal split between the two genetic lineages. Other studies in alpine ecosystems similarly found no changes in genetic diversity depending on population size or along altitudinal gradients [59]. Along environmental gradients, several parameters such as overlap of generations or recruitment frequency may influence intrapopulation genetic variation and gene flow among populations may be sufficient to maintain relatively high genetic diversity even in the smallest populations.
The longitudinal genetic pattern seemed to result from the migration of two Alpine lineages and from the establishment of populations in refugia along the Pyrenean chain due to glaciation events. The strong differentiation between the eastern and western Pyrenean lineages and the low level of gene flow observed between the two suggested that very few admixtures occurred since the two lineages had become established. This may result from topographical barriers (high summits, glaciated areas), seed dispersion by gravity and short pollinator travels hampering gene flow among populations and genetic homogenization. An east-west genetic separation had also been observed in other Pyrenean species (i.e. Papaver alpinum, [57], Trifolium alpinum, [19]). The populations of the eastern lineage were little differentiated even though they were distant from one another (35 km between the two undifferentiated populations CamI and EynH), were of very different sizes, grew in contrasting environments, and experienced a low level of amongpopulation current gene flow. Very small and isolated populations (PraL, LecL, LapL) showed a relatively high level of intrapopulation heterozygosity (H E , H O ) and exhibited private alleles (PraL, LapL) despite relatively lower allelic richness likely due to the loss of rare alleles by genetic drift, common in such small populations [60]. These populations survive in rather unusual ecological conditions for this species, i.e. on north facing slopes in deep, shaded valleys and mainly beneath beech forest (Fagus sylvatica) or sometimes in the vicinity of Mediterranean Quercus ilex forests. These data suggested that they were more likely to be rear-edge remnant populations [12] left there after the recolonization of highland areas (HYP 3) rather than the result of a secondary colonization of lowland areas by high altitude individuals. Together, these findings and the absence of clear founder effects suggested that the eastern populations within CL2 have experienced recent and extensive gene flow. Such intense gene exchanges could have occurred in large lowland refugia during the cold periods of the Pleistocene and possibly later in genetically connected populations growing over high altitude plateaus, which are typical of eastern Pyrenees landscapes. These large high altitude populations could have then become fragmented possibly during the last 30 ky when climate became warmer and drier [26].
The four well separated and differentiated genetic clusters nested in the western Pyrenean lineage appeared as a sub-structure of the east/west configuration and presumably diverged more recently. Overall, we did not observe a gradient of genetic variability loss among these clusters as expected from a rangeexpansion wave [11]. Rather, we observed a consistent genetic heterozygosity and allelic richness and a significant isolation by distance (Mantel test) among them. Wright's isolation by distance model suggests that gene flow is reduced as geographical distance increases. This may explain the very low level of admixture between CL3/CL6. The patchy distribution of often homogenous and highly distinct clusters enclosing a similar level of within-and between-population genomic variability is compatible with the hypothesis of the survival of the species in numerous Pyrenean refugia as proposed for other Pyrenean species [19], [18]. Moreover, the level of differentiation, the sharp genetic boundaries between most clusters and the number of migrants per generation suggested a long isolation of populations and the action of genetic drift in several high altitude refugia (HYP 2) rather than the survival in large and genetically connected populations in low elevation unglaciated areas. The higher genetic diversity in SpoI likely accounts for such refugia as well as TouH (having one private allele) where genetic variability was greater than in other populations. However, the relatively higher level of admixture between CL5 and CL4 or CL6 suggested that they had recently exchanged genes. Because the admixture concerned mainly isolated and sometimes very distant (AraH and CL5 populations) populations at high and intermediate altitude, contemporary gene exchanges among them through pollen or seed dispersal was very unlikely. We suspected that, at least in this part of the Pyrenees, longitudinal gene exchanges could have occurred possibly in lowland refugia across Pyrenean piedmonts before the post-glacial recolonization of high altitude area through valleys (HYP 1). Such a recolonizing pattern has been observed in sessile oak [51] and suggested for Aster pyrenaeus [20]. Thus, the phylogeographical pattern in central and western Pyrenees could result both from population isolation and differentiation in nunataks and gene exchange in some lowland refugia. Among clusters, CL6 showed the lowest genetic variability (EspI, SomI and PieI) and genetic diversity gradually decreasing in this cluster along the sequence TecL, SomI and PieI. Thus, the genetic pattern of CL6 could be chiefly related to founder events and the loss of genetic diversity and heterozygosity [11] during westward expansion of the species.

Implications for conservation
The longitudinal genetic structuration found in R. ferrugineum and particularly the east/west separation between the two main genetic lineages has already been observed in other studies [57], [19]. This suggests that many other alpine species may have been subject to similar events as R. ferrugineum and, nowadays, could present a comparable genetic pattern. We state that these data now should be usefully considered in management projects for Pyrenean species particularly those having an extended east/west distribution. According to the European Economic Community laws (92/43/EEC habitat directives, 1992) R. ferrugineum heathlands must be locally preserved due to their high ecological and patrimonial value. We think that, henceforth, in aim to preserve genetic diversity of the species, the genetic pattern here revealed should be taken into consideration to choose sites that have to be protected. The eastern lowland isolated populations (LecL, PraL, LapL) grow in very specific habitat and present rare alleles. According to the ''rear edge population hypothesis'', they harbor all together the bulk of species' genetic diversity [12]. Scientific programs have to be developed to follow the population dynamic and the genetic diversity evolution of these very original populations and take an inventory of other currently unknown rear-edge populations. Moreover, studying population dynamics at the edge of a species' range is crucial to understand the response of species to environmental changes and their susceptibility to extinction [31]. Also the populations TouH (presenting a rare allele) and SpoI (having a high allelic richness) should be considered as populations of interest for conservation. However, other parameters should be taken into account such as surrounding plant community or pollinator diversity. For example, it is known that the number of bumblebees species, which are one of the main pollinators of R. ferrugineum, is high in the eastern Pyrenees, especially close to population EynH [61]. Consequently, it could be an interesting population to monitor.
In conclusion, the Pyrenees are usually considered to only be a latitudinal barrier to northward or southward species expansions. In agreement with previous findings [18,19], our study of the subalpine shrub R. ferrugineum highlights a much more complex role of the Pyrenees in species genetic pattern and biogeography. The species showed two genetically isolated lineages, which have probably colonized the Pyrenees from the Alps and Jura Mountains. Our findings suggest that the eastern lineage most likely survived in large lowland refugia and that low-altitude small and isolated rear-edge populations have survived until now after the upward recolonization of highland areas agreeing with our HYP 3. Extensive gene flow may have occurred in lowland refugia as well as among populations growing in high altitude plateaus. The sub-structure of the central and western lineage suggests that the species mostly survived glaciations in isolated nunataks (HYP 2) but that some lowland refugia could also have played a role locally (HYP 1). Our study provides a useful genetic reference framework to select R. ferrugineum heathlands in management projects in the Pyrenees and points out priority populations for protection.