Skip to main content

The evolutionary history of Antirrhinumin the Pyrenees inferred from phylogeographic analyses



The origin and colonisation history after the Quaternary ice ages remain largely unresolved for many plant lineages, mainly owing to a lack of fine-scale studies. Here, we present a molecular phylogeny and a phylogeographic analysis of Antirrhinum, an important model system in plant biology, in the Pyrenees range. Our goal was to reconstruct the evolutionary and colonisation history of four taxa endemic to this region (A. majus subsp. majus, A. majus. subsp. striatum, A. molle, and A. sempervirens) by using a dense sampling strategy, with a total of 452 individuals from 99 populations whose collective distribution spans nearly the entirety of the Pyrenees and adjacent mountains.


Phylogenetic and phylogeographic analyses of the sequences of two plastid (trnS-trnG and trnK-matK) regions revealed the following: (i) historical relationship between the Pyrenees and Iberia (but not with the Alps); (ii) the long persistence of populations in the Pyrenees, at least since the Late Pleistocene; (iii) three different colonisation histories for populations from the Western, Central, and Eastern Pyrenees; (iv) the deep phylogeographic separation of the eastern and western populations; and (v) the colonisation of southern France from the Eastern Pyrenees.


The present study underlines the enormous influence of the glacial history of the mountain ranges on the current configuration of intra- and inter-specific genetic diversity in Antirrhinum, as well as the importance of periglacial areas for the survival of species during glacial periods of the Quaternary.


The three mountain systems (the Pyrenees, Alps, and Balkan mountains) located to the north of the major southern European peninsulas (Iberia, Italy, and Balkans) have played an important role in determining the current distributions of European biota, and they harbour exceptionnally high levels of biodiversity [13]. Their high elevation greatly influenced the climate-induced range shifts of many species during Quaternary climatic cycles [47], and local speciation promoted by isolation of populations in long-standing refugia and potential contact zones between postglacial recolonizing genetic lineages have been invoked to account for such diversity [8, 9].

The Pyrenees are especially renowned as an important plant diversity hotspot because of their high floristic richness (around 3500 species and subspecies of vascular plants) [10, 11], including considerable endemicity (about 4% of the plant species are endemic) and numerous plant associations [12, 13]. Yet we still know relatively little about the processes that led to species accumulation within this region. Their east–west arrangement clearly acted as a latitudinal barrier for postglacial dispersal predominantly out of the Iberian Peninsula. Moreover, they served as an arena in which postglacial lineage contacts occurred [9]. While the valleys in other European mountain ranges were only partially covered by ice even during the coldest periods of the ice age glaciations (e.g. Balkan Mountains and Carpathians), the Alps and the Pyrenees were almost entirely covered by large ice sheets [1417]. Such similarities in their glacial history may have resulted in analogous patterns of genetic structure of the plant species that inhabit the two east–west mountain ranges. Unfortunately, in contrast to the high number of studies that have focused on the evolutionary history of plant lineages in the Alps [see [18], few investigations have addressed the origin and colonisation history of Pyrenean plant species in a phylogenetic context at the appropriate spatial scale. Most of the previous phylogeographic studies that included the Pyrenean range have relied upon limited sampling strategies (but see [19]). Thus, the major patterns related to the origins and colonisation of Pyrenean plant lineages after the Quaternary ice ages remain largely unresolved.

The present study used Antirrhinum species from the Pyrenees to reconstruct the colonisation history of mountain plants in this region. Monophyly of this genus, an important model system in plant biology primarily distributed in the western Mediterranean basin, has been previously proved [20]. The Iberian Peninsula harbours the highest species diversity; however, a few Antirrhinum species are also found in other European areas, such as the Alps (A. latifolium) and Italy (A. siculum), or are distributed widely across the Mediterranean (A. tortuosum) despite the lack of any obvious long-distance dispersal syndrome. Historically, the number of Antirrhinum species has been the subject of taxonomic debate, resulting in different taxa circumscriptions. The difficulties in species delimitation, as well as the discrepancies between phylogenetic relationships and taxonomic classifications, have been interpreted as strong evidence for recent geographic speciation and extensive hybridisation since the Pliocene [2124].

Four Antirrhinum taxa from two sections are endemic to the Pyrenees and pre-Pyrenees range. The pre-Pyrenees range is formed by a complex system of foothill ranges which stretches from the western side (influenced by the Atlantic regime) to the Mediterranean coast on the eastern end of the Pyrenees (commonly called sub-Pyrenees). This range mostly runs parallel to the Pyrenees along the Iberian side, but it also includes the Corbières Range, towards the eastern end on the French side, where the Pyrenees’s slopes descend rather abruptly. The four endemic taxa occur in different habitats throughout these mountain ranges, ranging from lowland areas (800 m) to high mountain environments (>2500 m). A. molle and A. sempervirens (sect. Kickxiella) inhabit limestone cliffs from medium to high altitudes, whereas A. majus subsp. majus and A. majus ssp. striatum (sect. Antirrhinum) occur mostly in lowland areas of the Eurosiberian woodlands (see Additional file 1 for further clarification of the different taxonomic treatments). Previous phylogeographic studies found that four of the eight main Antirrhinum lineages are distributed in northeast Iberia [22]. Thus, the Pyrenees along with the pre-Pyrenees mountains are one of the areas with the highest numbers of genotypes and lineages. Recent studies have also been conducted to understand the cause and the maintenance of the parapatric distribution of two of these taxa (A. majus subsp. majus and A. majus subsp. striatum), as well as the ecological process that occurs in narrow contact zones [2527]. Thus, it seems reasonable to assume that the present-day distribution of Antirrhinum in the Pyrenees is the result of historical processes that have taken place across the mountain range and adjacent mountains during Quaternary times. We evaluate this possibility with molecular phylogenetic and phylogeographic analyses.

Based on previous phylogenetic and phylogeographic studies on plants (see Additional file 2), four hypotheses (not mutually exclusive) can be proposed for explaining the origin of Antirrhinum lineages that are found in the Pyrenees: (i) a predominantly southern (Iberian) origin; (ii) a predominantly northern (central European) origin; (iii) a long-standing isolation of Pyrenean populations; and (iv) recurrent secondary contacts between Iberian and central European lineages.

In this study, we first attempt to infer the historical connections between the Pyrenees and other geographic areas. We also investigate the genetic structure within the species and assess the patterns of genetic diversity within the Pyrenees region. By combining the results of these analyses, we explore potential processes that could explain the evolutionary history of Antirrhinum lineages throughout the Pyrenees.


Phylogenetic analysis

Detailed features of the two sequenced cpDNA regions are summarized in Table 1 (GenBank accession numbers are shown in Additional file 3). The total aligned length of the combined trnS-trnG/trnK-matK dataset of the 104 samples added to the 83 sequence matrix of Vargas et al. [22] was 1795 bp (excluding the outgroup). This revealed 82 variable sites of which 57 were parsimony-informative.

Table 1 Features of the two sequenced cpDNA regions

The 50% majority-rule consensus tree of the Bayesian analysis is shown in Figure 1. No conflicts are found between the Bayesian phylogenetic analyses for the trnK-matK and trnS-trnG separate matrices (Additional files 4 and 5). The strict consensus tree of the MP (maximum parsimony) analysis was basically congruent with the ML (maximum likelihood) analysis, although with a lower resolution and support values (Figure 1). One subclade recognized under Bayesian posterior probabilities (PP) contained all A. majus subsp. majus, A. majus subsp. striatum, A. sempervirens and A. molle accessions (PP =0.99) together with the accessions of species from southern Iberia (A. australe and A. cirrhigerum) and those species widespread across the Mediterranean (A. tortuosum and A. siculum). No analyses supported the monophyly of Pyrenean samples. The five samples of A. latifolium from the South-western Alps were clustered together (PP =1; ML-BS – maximum likelihood bootstrap – = 80%; MP-BP – maximum parsimony bootstrap percentages – = 63%) and clearly unrelated to Pyrenean samples. All three phylogenetic analyses showed high support values for clade 5 (PP =1.0; ML-BS = 95%; MP-BP = 94%), which comprised samples primarily distributed in the Eastern Pyrenees.

Figure 1

Phylogenetic analyses of Antirrhinum based on 190 trn S- trn G/ trn K- mat K sequences. The Bayesian tree is shown, in addition to values above branches for Bayesian posterior probabilities (PP) and parsimony bootstrap percentages (BP), values below branches indicate maximum likelihood bootstrap values (BS). A dash (-) above branches indicates disagreement between the maximum parsimony (MP) strict consensus tree and the Bayesian inference (BI) and maximum likelihood (ML) trees. Nodes with high support are indicated with an asterisk (*).

Estimation of divergence times

The values of standard deviation of the uncorrelated lognormal relaxed clock (0.31) and coefficient of variation (0.33) for rate heterogeneity within our trnS-trnG/trnK-matK dataset indicated low but significant rate heterogeneity among lineages (>0.1). Therefore the use of the uncorrelated clock was considered appropriate. Examination of the MCMC samples with Tracer 1.5 [28] confirmed adequate sample size, with ESS (effective sample size) values in the hundreds or thousands for both dating analyses. These also indicated adequate convergence and stationarity of the posterior probability distributions after discarding the burn-in. The topology of the maximum clade credibility (MCC) tree (Figure 2) was congruent with those of ML and MP analyses. The analysis revealed geographically structured lineages diverging since the Pliocene. A highly supported lineage (N1) containing samples of Northern Iberia (A. braun-blanquetii, A. graniticum, A. lopesianum, A. linkianum, A. litigiosum, A. pulverulentum and A. pertegasii), Mediterranean basin (A. siculum, A. australe, A. cirrhigerum, A. tortuosum) and the four Pyrenean species was estimated to split between 0.46 and 4.2 Ma – 95% Highest Posterior Density (HPD) – in the late Pliocene/early Pleistocene. A well-supported lineage primarily from the Eastern Pyrenees (N4) diverging between 0.13 and 1.7 Ma (100% HPD) was also observed. Additionally, we found highly supported sublineages of Pyrenean populations diverging in the Quaternary (see Figure 2).

Figure 2

Biogeographic reconstruction based on 190 concatenated trn S- trn G/ trn K- mat K sequences of Antirrhinum . The tree summarizes the geospatial Bayesian analysis. Pie charts represent posterior probability distributions of ancestral states at well-supported nodes. Colonisation routes supported by a BF > 3 are shown on the map (see Additional file 6). The inset shows the study mountain range (Pyrenees and neighboring areas). Colors in ellipses and branchs represent the areas of origin. Arrows specify directionality in the colonisation route, inferred from well-supported nodes of interest in the geospatial Bayesian analysis. Arrow width indicates relative support of migrations. Node bars represent the 95% highest posterior density intervals for the divergence time estimates. Numbers above branches are Bayesian posterior probabilities. Abbreviations used: L, Late Pleistocene; H, Holocene.

Ancestral area reconstructions

DPA analysis showed uncertainty in the geographic origin of the Pyrenean lineages; since several areas outside the Pyrenees showed similar probabilities (see Figure 2). Nevertheless, stronger support was found for the geographic origin of more recent subclades, leading to a robust reconstruction of dispersal events within and around the Pyrenees during the Pleistocene. Bayes factor (BF) tests (Additional file 6) for significant non-zero rates revealed a historically close connection between the Pyrenees and the Iberian Peninsula (BF = 21.93). In contrast, South-western Alps populations (where A. latifolium occurs) were not linked to the Pyrenees. The Central Pyrenees were the most likely ancestral location of western haplotypes (BF = 36.74). A colonisation episode involving a few populations may have also occurred from Central Pyrenees to the Southern French basin, as indicated by a well-supported connection between these two areas (BF = 10.45). However, the colonisation of the Southern French basin seems to have occurred mainly from the Eastern Pyrenees as indicated by a considerably higher Bayes factor value (BF = 404.17). Interestingly, a less supported link was also observed between Eastern and Central Pyrenees (BF = 4.13). These results suggest: (i) no close-relationship between Pyrenean lineages and those distributed in the South-western Alps; (ii) an origin of Pyrenean lineages during the Pleistocene; and (iii) recent dispersal episodes within the Pyrenees from two main spreading centers (Eastern and Central Pyrenees).

When the number of areas was reduced to four –Pyrenees and surrounding areas, South-western Alps, Iberian Peninsula (excluding the Pyrenees) and rest of Mediterranean areas– the inferred historical biogeographic scenarios from the S-DIVA (Statistical Dispersal-Vicariance Analysis) and DEC (Dispersal-Extinction-Cladogenesis) analyses were mostly congruent, but not very informative about the origin of the Pyrenean lineages (Additional files 7 and 8 respectively). However, the results agreed with a sister relationship between the South-western Alps lineage and an Iberian clade, indicating: (i) a historical connection between two distant areas and (ii) no recent relationship between the Alps and Pyrenean lineages.

Haplotype data analysis

Antirrhinum network

The statistical parsimony analysis of the Antirrhinum matrix found 11 haplotype lineages, of which seven included samples from NE Iberia. In addition, five haplotype lineages from the Pyrenees were formed by 18 haplotypes (haplotypes 2, 9-12 and 52 - 63, see Additional file 9) obtained from the 99 populations collected across the Pyrenees and adjacent areas. Some missing haplotypes (extinct or not found) and one loop were retrieved for the whole genus. Most of the haplotypes were grouped primarily together within lineages I and II of northern Iberia (Additional file 9). The five samples collected in the South-western Alps provided two new haplotypes (haplotypes 64 and 65) with a tip position, which formed part of lineage VI from southeastern Iberia.

Pyrenees network

The haplotype network analysis of the Pyrenees matrix distinguished 13 haplotypes (hereafter named H1–H13) distributed into five lineages (hereafter named lineages A–E). The H12 resulted from a 9 bp deletion. Neither loops nor missing haplotypes (extinct or not found) were observed (Figure 3c). No species had exclusive haplotypes, except for A. sempervirens that had four (H2, H3, H7 and H8) (Additional file 10).

Figure 3

Distribution of haplotypes and lineages based on 452 trn S- trn G sequences. (3a) Haplotypes distributed across the Pyrenees and delimited by a grid of 10x10 km quadrats in which haplotypes of all the populations at the same quadrat are grouped together (pie charts). (3b) Distribution of cpDNA lineages. (3c) Haplotype network of Pyrenean populations. Labels on the lines connecting haplotypes represent the plastid lineages. The dashed black line represents the approximate location of the genetic boundary indicated by SAMOVA. Yellow lines delimit the historical phytogeographic areas (Western, Central and Eastern Pyrenees). The blue shadow represents the approximate location of the pre-Pyrenees range.

Genetic structure and diversity

The geographical distribution of haplotypes and lineages is shown in Figure 3 (information per grids summarized in Additional file 10). The Eastern Pyrenees harbored only haplotypes of lineage E, which appeared to be primarily distributed in this area. The Western Pyrenees harbored exclusively the haplotype 4. The five lineages were distributed across the Central Pyrenees. The majority of 10×10 km quadrats (54) harbored only one haplotype, although some contained two (8 quadrats) and three (3 quadrats) haplotypes. Moreover, few quadrats had haplotypes of two (grid coordinates I14, I25, J21, J22, K19) and three (J17, F8) lineages, and all of them fell into the Central Pyrenees.

Haplotype frequencies and molecular diversity indices per region are summarized in Tables 2 and 3. Central Pyrenees contained the highest genetic diversity (Hd = 0.789, π = 0.004). In this area haplotype 5 occurred in the highest percentage (32.85%), followed by H11 (21.30%) and H4 (20.22%). Eastern Pyrenees showed a lower value of genetic diversity in comparison with Central Pyrenees (Hd = 0.529, π = 0.0012). Only haplotypes of lineage E were found in the Eastern Pyrenees, with H11 showing the highest percentage (56.56%). The lowest molecular diversity was found in Western region, with H4 as the only haplotype.

Table 2 Number and frequency of each haplotype (%)
Table 3 Genetic diversity in the three Pyrenean sections

FST values indicated that a phylogeographical pattern exists between the three regions (Table 4). The highest FST value (0.620) was found for genetic differentiation between Lineage E and the rest of lineages. AMOVA (analysis of molecular variance) analysis showed significant values for genetic differentiation among groups for the following two comparisons: (i) sampling locations partitioned according to Eastern, Western, and Central regions; and (ii) according to lineage E versus rest of samples. Nevertheless, variation values were higher among populations within groups (Table 5). The results of SAMOVA (spatial analysis of molecular variance) revealed two high FCT values for clusters not including groups of a single population (Additional file 11). Variance among groups reached a value over 80% for the partition into seven geographic groups corresponding mainly to different plastid lineages and sublineages, without a clear geographic structure. The division into K = 5 clusters had the second highest variance among groups (78.77%). Both divisions clearly differentiated lineage E in a separate group. Such split largely corresponded to the biogeographic boundary for Eastern Pyrenees, which can be visualized in the distribution map of haplotypes (Figure 3). With K ranging between 8 and 20, FCT values did not increase significantly and in most cases the newly defined groups comprised single populations. SAMOVA showed that a substantial portion of the cpDNA genetic variability was found among groups, whereas the differences among populations within groups accounted for less variation.The clustering BAPS (Bayesian Analysis of Population Structure) analysis resulted in a best partition of K = 8. These eight clusters primarily corresponded with the plastid lineages and sublineages (see Figure 1). In the admixture analysis, all individuals were unambiguously assigned to their respective group without any probability of being misplaced.

Table 4 Pairwise F ST -values
Table 5 Results of AMOVA

Demographic history

Neither Tajima’s D nor Fu’s FS was significantly different from zero for lineage E (Table 6). Nevertheless, the mismatch distribution analyses did not reject the sudden expansion model for this lineage. The unimodal peak and nonsignificant Harpending’s raggedness index (Figure 4) are indicative of recent demographic expansion [29]. For the group of the remaining lineages, Tajima’s D and Fu’s Fs were negative but not significant. Mismatch distribution analyses rejected the range expansion hypothesis for this group by a bimodal distribution and significant SSD (sum of square deviations) and RAG (Harpending’s raggedness index) values (Figure 4).

Table 6 Tajima’s D and Fu’s FS test
Figure 4

Results of mismatch distribution analyses. Pictures in right are results of mismatch distribution. Goodness of fit of the observed vs. the theoretical mismatch distributions under a sudden expansion model was tested using the sum of squared deviation (SSD) and raggedness index. Not significant P > 0.10.


Origins of Antirrhinumin the Pyrenees

Regarding the geographic origin of the Pyrenean lineages, an important historical connection between the northeast Iberian Peninsula and the Pyrenees is well supported (Figure 2, Additional file 6). It is notable that none of the analyses identified evolutionary relationships between the lineages from the Pyrenees and those from the South-western Alps, which is a connection found in other plant species with similar distributions e.g. Erodium, [30], Primula, [31]. The lack of a connection between these two areas is also supported by nrITS sequences (Liberal et al. unpublished, see additional file 12). This is a rather unexpected, but essential, result since the yellow-flowered populations that inhabit the Alps and the Pyrenees have recently been circumscribed into a single species (A. latifolium Mill.) by some taxonomists [32, 33]. The phylogenetic and dating analyses provide strong support for the divergence of Pyrenean sublineages in the Quaternary but poor support for sister-group relationships at the most basal nodes. The analyses presented here and those of Vargas et al. [22] support the first hypothesis herein posit that the Pyrenean lineages are the result of long-standing isolation and geographic differentiation from a widespread ancestor in the Pleistocene. However, additional nuclear data, such as single-copy nuclear markers, would be required to determine whether the unresolved basal relationships are attributable to hard or soft polytomies [34].

East–west geographical structure of the genetic diversity

The main genetic discontinuity identified by the analyses (SAMOVA, AMOVA, and FST) indicated the differentiation of the Lineage E from the remaining lineages. These two groups basically correspond to the division between the Eastern and Central-West Pyrenees (Figure 3). This cryptic transition area, which approximately follows the abrupt valley of the Segre River, is supported by the phylogeographic gap in the southeastern part of the Pyrenees that was already identified in previous studies of animals and plants e.g. [3538]. A similar phytogeographic boundary was previously proposed based on species distributions in the Pyrenees [39, 40]. The high level of differentiation between the Eastern Pyrenean lineage (Lineage E) and the remaining lineages, suggests the isolation of eastern populations over more than one glacial cycle.

Interestingly, the eastern and central areas may have experienced different historical processes during the last glaciations. According to the results of the mismatch distribution analyses, the observed genetic structure in the Eastern Pyrenees is congruent with a range expansion process, whereas this pattern was not observed in the Central Pyrenees. By contrast, the strong evidence for the genetic clustering of Central Pyrenean samples, which was detected by the BAPS and SAMOVA analyses, as well as the significant genetic differentiation among the different lineages that occur in the area, support long-standing genetic isolation in different local refugia without an extensive interglacial (or postglacial) range expansion as already inferred for Pinus uncinata [according to [19]. However, the homogeneously distributed genetic diversity indicates that moderate gene flow occurred among Central Pyrenean refugia during interglacial periods. Thus, two different historical processes were inferred: (i) limited extinction/expansion during glacial cycles in the central area, and (ii) population expansion (probably pre-dated by a demographic bottleneck based on the reduced number of lineages and haplotypes) in the eastern area, may explain the observed differences in the genetic structure and diversity of both regions. A similar interpretation was proposed for Ramonda myconi to explain genetic differences between populations from Central and Eastern Pyrenees [41] where significant evidence was presented to support the existence of an important refuge near the central area. Further support for a large refuge in Central Pyrenees was also found by palynological data [42, 43], the patterns of phylogeographic concordance among different species [44], and the high diversity of habitats, species richness, and endemic species [63 of the 160 endemic plant species of the Pyrenean flora are exclusively from the Central Pyrenees; [45]. In agreement with this and other studies e.g. [46], a scenario of multiple refugia in mountains adjacent to the Central Pyrenees (i.e., the southern mountains of the pre-Pyrenees) is gaining more support (see Figure 3).

Similar to the pattern found in Antirrhinum, Ramonda myconi also exhibits considerable differentiation but reduced genetic diversity in the eastern population [41]. Fossil pollen records and molecular studies suggest the presence of smaller refugia in the southeastern Pyrenees, along the Mediterranean coast [4749]. Recent studies indicate that climatic fluctuations during the last glacial period affected the Eastern Pyrenees (Mediterranean coast) more strongly than the Central Pyrenees, at least during the last 50,000 years. Overall, more pronounced oscillations between arid and humid periods prevailed in the easternmost regions [50, 51]. Thus, severe reductions in population sizes and the isolation of peripheral populations at the edge of the distribution might have led to genetic drift and local adaptations via selective forces see (see [5255]). Future comparative phylogeographic studies and coalescent demographic analyses could help to elucidate the historical processes responsible for the phylogeographic and floristic configuration of the Pyrenees.

Migration routes across the Pyrenees

The colonisation of the Western Pyrenees occurred from the Central Pyrenees. Indeed, the western populations have a low level of genetic diversity, which is shared by central populations. Such gradual genetic impoverishment observed from the Central to the Western Pyrenees is interpreted as a recent colonisation, which originated in the most central parts of the mountain range see [4, 5, 56, 57]. The results of the DPA analysis also indicate that the postglacial colonisation of the southern French basin occurred primarily from the Eastern Pyrenees (Figure 3). The geographic distribution of the haplotypes and the mismatch distribution analyses (Figure 4) suggest a gradual colonisation of southern France from northeastern Iberia, where the edge of the glacial distribution range may have retreated early see [58].

Similarities between the Pyrenees and the Alps

The causes of the separation of Antirrhinum lineages and plant species by the abrupt valley in the Pyrenees (Segre River) are unclear. Similar west–east vicariant lineages abutting along the Etsch Valley and the Brenner Pass in the Alps have been reported [18, 5961]. Indeed, this discontinuity accounts for the most important phytogeographic and biogeographic boundary in the Alps e.g. [6264]. Similar genetic and floristic discontinuities are found in the Alps and the Pyrenees, which suggest that plant extinctions due to Holocene glaciers and further recolonisation events were not distributed homogeneously across each of the two mountain ranges. The major effect of the glacial history of mountain valleys in the current configuration of intra- and inter-specific genetic diversity requires further investigation.

The present study also underlines the importance of periglacial areas for the survival of species during glacial periods of the Quaternary. The pre-Pyrenees range, is recognised as an important glacial refuge where many animal and plant species survived during cold episodes [30, 38, 65, 66]. Previous studies focused on the Alps have also recognised the importance of marginal and peripheral areas at lower altitudes with warmer climates as refugia for mountain species, where their present-day distribution were almost entirely glaciated reviewed in [18].


Our results do not support connections between the Antirrhinum populations of the Alps and those of the Pyrenees. Instead, the analyses agree with ancient connection with Iberia followed by persistence of Antirrhinum populations in the Pyrenees, at least since Late Pleistocene times. The three Pyrenean sections (Western, Eastern, and Central) have different colonisation histories. The Western Pyrenean populations are closely related to the central populations, which is here interpreted as evidence for a recent westward colonisation (see Figure 2). It is likely that the Central pre-Pyrenees played an important role in the maintenance of the genetic diversity of Antirrhinum by providing suitable conditions for the establishment of probably separate refugial populations during glacial periods. A significant influx of migrants from adjacent refugia (Central pre-Pyrenees) during the interglacial and postglacial periods would account for the genetic diversity found in the Central Pyrenees. However, the current genetic structure of Antirrhinum in the Eastern Pyrenees appears to have been shaped by range expansion, which was probably preceded by range contraction at the southernmost edge of this area, from where little gene flow may have occurred into the Central Pyrenees over several glacial cycles. According to this, the Eastern pre-Pyrenees also appear to have been a refugial area, although less important than the central pre-Pyrenees, from where Antirrhinum colonisation southern France.


Study species

Maternally inherited plastid markers often display stronger differentiation between populations than biparental markers due to their smaller effective size and lack of recombination [56, 6769]. Because of this, reconstruction of migration histories by means of seed dispersal have been traditionally based on organelle markers once maternal inheritance has been documented [70, 71]. To this end, artificial crossing experiments under greenhouse conditions were conducted between four Antirrhinum taxa (A. mollissimum, A. majus subsp. majus, A. controversum and A. charidemi) (Additional file 13).

Sampling strategy and DNA sequencing

We collected a total of 452 individuals from 99 populations of A. molle, A. sempervirens, A. majus subsp. majus and A. majus subsp. striatum from the Pyrenees and adjacent mountains. In addition, a total of five individuals from five populations of A. latifolium distributed across the Province/southern Maritime Alps (Additional file 3) were analyzed given that this species has historically been considered closely related to A. majus subsp. striatum [32, 33]. The number of populations sampled per species depended on their distribution range: A. molle (3), A. sempervirens (9), A. majus subsp. majus (61) and A. majus subsp. striatum (26) respectively. Particular effort was made to achieve a sampling as complete as possible for the entire mountain range (see Figure 3). All individuals were collected in the field and dried in silica gel. Total genomic DNA was extracted using Dneasy Plant Mini Kit (QIAGEN Inc., California).

Two sequencing strategies were adopted. First, trnK-matK and trnS-trnG intergenic spacer sequences were obtained for one individual from each population, and added on to the sequence data matrix from Vargas et al. [22]. This dataset (hereafter Antirrhinum matrix) was used to infer the colonisation history of Pyrenean lineages by phylogenetic and phylogeographic analyses. Second, trnS-trnG intergenic spacer sequences were obtained for the 452 individuals (1-6 individuals per population) from the Pyrenees. This sequence matrix (hereafter Pyrenees matrix) was used to investigate underlying spatial genetic structure and genetic diversity across the Pyrenees.

The trnK-matK and trnS-trnG regions were amplified as in Vargas et al. [22]. PCR products were outsourced for sequencing to a contract sequencing facility (Macrogen, Seoul, South Korea) on an ABI Prism® 3730xi DNA sequencer, using the same primer sets as for PCR. Resulting sequence data were assembled and edited using Geneious Pro v5 [72]. Both datasets were aligned with the MAFFT v.6.814b alignment tool [73], as implemented in Geneious Pro with default parameters. Further adjustments were made by visual inspection.

Phylogenetic analyses (Antirrhinum matrix)

Bayesian phylogenetic analyses were performed on the two separate matrices (trnS-trnG/trnK-matK) to examine plastid gene tree congruence. Gambelia speciosa and Misopates orontium were selected as the outgroup based on previous phylogenetic evidence [20]. In addition, phylogenetic analyses for 190 trnS-trnG/trnK-matK concatenated sequences were conducted using Bayesian inference (BI), as implemented in MrBayes v3.1.2 [74], maximum likelihood (ML), as implemented in PhyML 3.0 [75], and maximum parsimony (MP), as implemented in TNT 1.1 [76] (see Additional file 14 for details).

Estimation of divergence times (Antirrhinum matrix)

A relaxed molecular-clock approach implemented in BEAST v.1.6.2 [77] was used to estimate the divergence times between Antirrhinum Pyrenean lineages. This software estimates the phylogenetic tree and node ages simultaneously. Since no reliable fossils of Antirrhinum are known to date, we implemented a secondary basal calibration following [78] but with some modifications. The divergence time between the basalmost Antirrhinum lineages (crown node) was modelled as a normal distribution with mean = 4.5 Ma and standard deviation = 1.8, on the basis of a previous analysis of ndhF sequences of Antirrhineae that included two most distant Antirrhinum species (A. meonantum and A. majus subsp. majus, following results of phylogenetic analyses). The latter analysis incorporated a calibration of 74 Ma for the divergence time between Oleaceae and Antirrhineae [79], and minimum stem-age constraints for Lamiales families and tribes based on five fossils see [80] for details on fossils].

Ancestral area reconstructions (Antirrhinum matrix)

A discrete phylogeographic analysis (DPA) implemented in BEAST [81] was performed to assess the probability distribution of the geographic locations in each node. A total of 14 discrete areas were delimited: (i) the four Iberian quadrants (northeastern Iberia, NE; northwestern Iberia, NW; southeastern Iberia, SE; southwestern Iberia, SW), as divided by the geographical coordinates 40ºN/5ºW see [22]; (ii) Eastern, Central and Western Pyrenees, as the three recognized biogeographic regions within the Pyrenees (see below); (iii) the other two northern areas sampled nearby the Pyrenees (Southern French basin and South-western Alps); and (iv) the remaining five areas sampled across Mediterranean basin (northern Africa, Sicily, Sardinia, Italian Peninsula and Turkey). Statistical significance for the rates of the dispersal events was assessed via Bayes factor (BF) test as described by Lemey et al. [81]; see Additional file 14 for details].

Additional ancestral range reconstructions were conducted using the Bayesian time-calibrated molecular phylogeny with the aim to discriminate between northern and southern origin of Pyrenean lineages. For this purpose, only four areas were delimited (i) Iberian Peninsula, (ii) Pyrenees and adjacent areas, (iii) South-western Alps, and (iv) samples from the Mediterranean basin (excluding Iberia). Ancestors were allowed to be present in all of them. Distribution ranges of sequences (haplotypes) instead of species was used see [80]. Two alternative reconstruction methods were used: (a) Statistical Dispersal-Vicariance Analysis (S-DIVA) implemented in the program RASP 1.1 [82], and (b) dispersal-extinction-cladogenesis analysis (DEC) implemented in the software package Lagrange v2.0.1 [83]; see Additional file 14 for details].

Haplotype network analyses (Antirrhinum and Pyrenees matrices)

Haplotype networks using the Antirrhinum and Pyrenees matrices (hereafter Antirrhinum and Pyrenees networks, respectively) were constructed under the statistical parsimony framework in TCS 1.2.1 [84]. The maximum number of differences among haplotypes, as a result of single substitutions, was calculated with 95% confidence limits. Indels of varying lengths resulting from the expansion/contraction of polyA and polyT were excluded due to their high levels of homoplasy [85]. One indel of 9 bp was found in several samples of the population Villefranche-de-Conflent of A. majus. subsp. striatum for the trnS-trnG plastid region (Additional file 3) and was re-coded as a single character and treated as a fifth character state in. For the Pyrenees network, one sequence of A. braun-blanquetii (EU673480) was used as the outgroup, since this sequence appears to retain the ancestral haplotype in the Antirrhinum network analysis.

Genetic diversity and geographic structure (Pyrenees matrix)

An analysis of genetic diversity was carried out across the three recognized biogeographic regions in which the Pyrenees range is divided (Eastern, Central and Western Pyrenees) (see Figure 3a). The boundaries of these three biogeographic areas, although with slight differences, have been traditionally established by both geologists [86, 87] and phytogeographers [39, 40, 8891] on the basis of geologic, climatic and floristic data. The boundaries of these areas are approximately defined by the Gallego River, between Central and Western Pyrenees, and the Segre River, between Central and Eastern Pyrenees. Haplotype frequencies and molecular diversity indices for each biogeographic area were calculated using DnaSP v5 [92]. In addition, to identify potential hotspots of genetic diversity across the Pyrenees, individuals were geographically grouped by means of a 10×10 km grid. Charts representing haplotype frequencies were constructed for each grid cell, which was named by a generic letter–number code (Figure 3).

To identify genetic subdivisions among the Eastern, Central and Western Pyrenees, we performed an analysis of molecular variance (AMOVA) [93]. Pairwise FST statistics were also calculated to estimate genetic distances. Both analyses were performed by using ARLEQUIN [94]. Additionally, an AMOVA was performed in order to assess the partitioning of variance between the Lineage E (primarily distributed in the eastern part of the Pyrenees) and the rest of lineages (see below). Furthermore, to evaluate the optimal grouping of the sampled sites without a priori assumptions, a Bayesian analysis of population structure, implemented in the BAPS software, version 5.3 [95, 96], and a spatial analysis of molecular variance (SAMOVA) implemented in the software package SAMOVA 1.0 [97] were also performed (see Additional file 13 for details on analyses).

Demographic history (Pyrenees matrix)

Past range expansions were assessed for the two partitioned population groups (lineage E and remaining populations, see below) using two methods. First, we tested for range expansion using Tajima’s D [98] and Fu’s Fs statistics [99] which, assuming neutrality, are expected to be significantly negative under population expansion. Significance for each test was assessed by generating null distributions from 10000 simulations. Second, we also conducted a mismatch analysis of cpDNA sequence differences, comparing the observed frequencies of pairwise differences of haplotypes with those expected under a sudden expansion (pure demographic expansion) model [29, 100]. Under expansion, the distribution of the observed differences is expected to be unimodal, whereas, under population contraction or genetic subdivision, a multimodal distribution is expected. Statistical significance was determined by 10000 bootstrap replicates. Goodness of fit was assessed by the sum of square deviations (SSD) and the Harpending’s raggedness index [101] between the observed and the expected mismatch. These analyses were conducted using ARLEQUIN 3.5. [102].


  1. 1.

    Varga ZS, Schmitt T: Types of oreal and oreotundral disjunctions in the western Palearctic. Biol J Linn Soc. 2008, 93 (2): 415-430.

    Google Scholar 

  2. 2.

    Barthlott W, Mutke J, Rafiqpoor D, Kier G, Kreft H: Global Centers Of Vascular Plant Diversity. Plant Taxonomy to Evolutionary Biology Nova Acta Leopoldina NF 92, Nr 342. 2005, 61-83.

    Google Scholar 

  3. 3.

    EEA: Biogeographical Regions, Europe 2011. 2011, Copenhagen: European Environment Agency

    Google Scholar 

  4. 4.

    Hewitt G: The genetic legacy of the Quaternary ice ages. Nature. 2000, 405 (6789): 907-913.

    PubMed  CAS  Google Scholar 

  5. 5.

    Taberlet P, Fumagalli L, Wust-Saucy A-G, Cosson J-F: Comparative phylogeography and postglacial colonization routes in Europe. Mol Ecol. 1998, 7 (4): 453-464.

    PubMed  CAS  Google Scholar 

  6. 6.

    Schmitt T: Biogeographical and evolutionary importance of the European high mountain systems. Front Zool. 2009, 6: 9-

    PubMed  PubMed Central  Google Scholar 

  7. 7.

    Hewitt G: Genetic consequences of climatic oscillations in the Quaternary. Philos Trans R Soc Lond Ser B Biol Sci. 2004, 359 (1442): 183-195.

    CAS  Google Scholar 

  8. 8.

    Médail F, Diadema K: Glacial refugia influence plant diversity patterns in the Mediterranean Basin. J Biogeogr. 2009, 36 (7): 1333-1345.

    Google Scholar 

  9. 9.

    Hewitt GM: Post-glacial re-colonization of European biota. Biol J Linn Soc. 1999, 68 (1–2): 87-112.

    Google Scholar 

  10. 10.

    Dupias G: Végétation Des Pyrénées. Notice Détaillée De La Partie Pyrénéenne De La Carte De La Végétation De La France au 1/200.000e. 1985, Toulouse, FR: CNRS

    Google Scholar 

  11. 11.

    Gómez D, Sesé JA, Villar L: The vegetation of the Alpine Zone in the Pyrenees. Alpine Biodiversity In Europe. Edited by: Nagy L, Grabherr G, Koörner C, Thompson DBA. 2003, Berlin, DE: Springer, 85-92. [Ecological Studies no.167]

    Google Scholar 

  12. 12.

    Villar L, Dendaletche C: Pyrenees. France, Spain and Andorra. Centres Of Plant Diversity A Guide And Strategy For Their Conservation. Edited by: Heywood VR, Davis SD, Hamilton AC. 1994, Cambridge: WWF and IUCN, 61-64.

    Google Scholar 

  13. 13.

    Puşcaş M, Choler P: A biogeographic delineation of the European Alpine System based on a cluster analysis of Carex curvula dominated grasslands. Flora Morphol Distrib Funct Ecol Plants. 2012, 207: 168-178.

    Google Scholar 

  14. 14.

    Charlesworth JK: The Quaternary Era: With Special Reference To Its Glaciation. 1957

    Google Scholar 

  15. 15.

    Frenzel B, Pecsi M, Velichko AA: Atlas of the Paleoclimates and paleoenvironments of the Northern Hemisphere. Late Pleistocene Holocene. Gustav Fisher, Budapest. 1992

    Google Scholar 

  16. 16.

    Vuia F: Studiul Reliefului Glaciar oei Periglaciar din România. Thesis raport, Babes-Bolyai University, Cluj-Napoca. 2005

    Google Scholar 

  17. 17.

    Messerli B: Die eiszeitliche und die gegenwärtige Vergletscherung im Mittelmeerraum. Geogr Helv. 1967, 22: 105-228.

    Google Scholar 

  18. 18.

    Schönswetter P, Stehlik I, Holderegger R, Tribsch A: Molecular evidence for glacial refugia of mountain plants in the European Alps. Mol Ecol. 2005, 14 (11): 3547-3555.

    PubMed  Google Scholar 

  19. 19.

    Dzialuk A, Muchewicz E, Boratyński A, Montserrat JM, Boratyńska K, Burczyk J: Genetic variation of Pinus uncinata (Pinaceae) in the Pyrenees determined with cpSSR markers. Plant Syst Evol. 2009, 277 (3–4): 197-205.

    CAS  Google Scholar 

  20. 20.

    Vargas P, Rosselló JA, Oyama R, Güemes J: Molecular evidence for naturalness of genera in the tribe Antirrhineae (Scrophulariaceae) and three independent evolutionary lineages from the New World and the Old. Plant Syst Evol. 2004, 249 (3–4): 151-172.

    CAS  Google Scholar 

  21. 21.

    Jiménez JF, Sánchez-Gómez P, Güemes J, Rosselló JA: Phylogeny of snapdragon species (Antirrhinum; Scrophulariaceae) using non-coding cpDNA sequences. Isr J Plant Sci. 2005, 53 (1): 47-54.

    Google Scholar 

  22. 22.

    Vargas P, Carrió E, Guzmán B, Amat E, Güemes J: A geographical pattern of Antirrhinum (Scrophulariaceae) speciation since the Pliocene based on plastid and nuclear DNA polymorphisms. J Biogeogr. 2009, 36 (7): 1297-1312.

    Google Scholar 

  23. 23.

    Carrió E, Forrest AD, Güemes J, Vargas P: Evaluating species nonmonophyly as a trait affecting genetic diversity: A case study of three endangered species of Antirrhinum L. (Scrophulariaceae). Plant Syst Evol. 2010, 288 (1): 43-58.

    Google Scholar 

  24. 24.

    Wilson Y, Hudson A: The evolutionary history of Antirrhinum suggests that ancestral phenotype combinations survived repeated hybridizations. Plant J. 2011, 66 (6): 1032-1043.

    PubMed  CAS  Google Scholar 

  25. 25.

    Andalo C, Cruzan MB, Cazettes C, Pujol B, Burrus M, Thébaud C: Post-pollination barriers do not explain the persistence of two distinct Antirrhinum subspecies with parapatric distribution. Plant Syst Evol. 2010, 286 (3): 223-234.

    Google Scholar 

  26. 26.

    Khimoun A, Cornuault J, Burrus M, Pujol B, Thébaud C, Andalo C: Ecology predicts parapatric distributions in two closely related Antirrhinum majus subspecies. Evol Ecol. 2013, 27 (1): 51-64.

    Google Scholar 

  27. 27.

    Tastard E, Ferdy JB, Burrus M, Thébaud C, Andalo C: Patterns of floral colour neighbourhood and their effects on female reproductive success in an Antirrhinum hybrid zone. J Evol Biol. 2012, 25 (2): 388-399.

    PubMed  CAS  Google Scholar 

  28. 28.

    Rambaut A, Drummond AJ: Tracer 1.5.0. 2007, Available at:

    Google Scholar 

  29. 29.

    Rogers AR, Harpending H: Population growth makes waves in the distribution of pairwise genetic differences. Mol Biol Evol. 1992, 9 (3): 552-569.

    PubMed  CAS  Google Scholar 

  30. 30.

    Alarcón M, Vargas P, Sáez L, Molero J, Aldasoro JJ: Genetic diversity of mountain plants: Two migration episodes of Mediterranean Erodium (Geraniaceae). Mol Phylogenet Evol. 2012, 63 (3): 866-876.

    PubMed  Google Scholar 

  31. 31.

    Zhang LB, Comes HP, Kadereit JW: The temporal course of Quaternary diversification in the European high mountain endemic Primula sect. Auricula (Primulaceae). Int J Plant Sci. 2004, 165 (1): 191-207.

    Google Scholar 

  32. 32.

    Güemes J: Antirrhinum L. Flora Iberica. Edited by: Castroviejo S, Herrero A, Benedí C, Rico E, Güemes J. 2009, Madrid: CSIC, vol. 13

    Google Scholar 

  33. 33.

    Sutton DA: A revision of the tribe Antirrhineae. 1988, Oxford: Oxford University Press

    Google Scholar 

  34. 34.

    Walsh H, Kidd M, Moum T, Friesen V: Polytomies and the power of phylogenetic inference. Evolution. 1999, 53: 932-937.

    Google Scholar 

  35. 35.

    Magri D, Vendramin GG, Comps B, Dupanloup I, Geburek T, Gömöry D, Latałowa M, Litt T, Paule L, Roure JM, Tantau I, Van Der Knaap WO, Petit RJ, De Beaulieu JL: A new scenario for the quaternary history of European beech populations: palaeobotanical evidence and genetic consequences. New Phytol. 2006, 171 (1): 199-221.

    PubMed  CAS  Google Scholar 

  36. 36.

    Albre J, Gers C, Legal L: Molecular phylogeny of the Erebia tyndarus (Lepidoptera, Rhopalocera, Nymphalidae, Satyrinae) species group combining CoxII and ND5 mitochondrial genes: a case study of a recent radiation. Mol Phylogenet Evol. 2008, 47 (1): 196-210.

    PubMed  CAS  Google Scholar 

  37. 37.

    Kropf M, Kadereit JW, Comes HP: Differential cycles of range contraction and expansion in European high mountain plants during the late quaternary: insights from Pritzelago alpina (L.) O. Kuntze (Brassicaceae). Mol Ecol. 2003, 12 (4): 931-949.

    PubMed  CAS  Google Scholar 

  38. 38.

    Sánchez de Dios R, Benito-Garzón M, Sainz-Ollero H: Hybrid zones between two European oaks: a plant community approach. Plant Ecol. 2006, 187 (1): 109-125.

    Google Scholar 

  39. 39.

    Rivas Martinez S: Memoria Del Mapa De Series De Vegetaci6n De Espafia 1: 400000. 1987, Madrid: ICONA

    Google Scholar 

  40. 40.

    Vigo J, Ninot J: Los Pirineos. La Vegetación de España. Edited by: Peinado Lorea M, Rivas Martinez S. 1987, Alcalá de Henares: Col. Aula Abierta. Publ. Univ, 351-384.

    Google Scholar 

  41. 41.

    Dubreuil M, Riba M, Mayol M: Genetic structure and diversity in Ramonda myconi (Gesneriaceae): effects of historical climate change on a preglacial relict species. Am J Bot. 2008, 95 (5): 577-587.

    PubMed  CAS  Google Scholar 

  42. 42.

    Jalut G, Amat AE, Riera I, Mora S, Fontugne M, Mook R, Bonnet L, Gauquelin T: Holocene climatic changes in the western mediterranean: installation of the mediterranean climate. C R Acad Sci II Sci Terre Planetes. 1997, 325 (5): 327-334.

    Google Scholar 

  43. 43.

    Valero-Garcés BL, González-Sampériz P, Delgado-Huertas A, Navas A, Machín J, Kelts K: Lateglacial and Late Holocene environmental and vegetational change in Salada Mediana, central Ebro Basin, Spain. Quat Int. 2000, 73–74: 29-46.

    Google Scholar 

  44. 44.

    Gómez A, Lunt DH: Refugia within refugia: Patterns of phylogeographic concordance in the Iberian Peninsula. Phylogeography of Southern European Refugia. Edited by: Weiss S, Ferrand N. 2007, NL: Dordrecht, Springer, 155-188.

    Google Scholar 

  45. 45.

    Villar L, Sesé JA, Ferrández JV: Atlas de la flora del Pirineo Aragonés. II. (Pyrolaceae-Orchidaceae. Síntesis). 1997, Huesca, Spain: Consejo de Protección de la Naturaleza de Aragón & Instituto de Estudios Altoaragonenses

    Google Scholar 

  46. 46.

    Segarra-Moragues JG, Palop-Esteban M, González-Candelas F, Catalán P: Nunatak survival vs. tabula rasa in the Central Pyrenees: a study on the endemic plant species Borderea pyrenaica (Dioscoreaceae). J Biogeogr. 2007, 34 (11): 1893-1906.

    Google Scholar 

  47. 47.

    Pèrez-Obiol R, Julià R: Climatic change on the iberian peninsula recorded in a 30,000-yr pollen record from lake banyoles. Quat Res. 1994, 41 (1): 91-98.

    Google Scholar 

  48. 48.

    Salvador L, Alía R, Agúndez D, Gil L: Genetic variation and migration pathways of maritime pine (Pinus pinaster Ait) in the Iberian peninsula. Theor Appl Genet. 2000, 100 (1): 89-95.

    Google Scholar 

  49. 49.

    Lumaret R, Mir C, Michaud H, Raynal V: Phylogeographical variation of chloroplast DNA in holm oak (Quercus ilex L.). Mol Ecol. 2002, 11 (11): 2327-2336.

    PubMed  CAS  Google Scholar 

  50. 50.

    Sánchez Goñi M, Cacho I, Turon J, Guiot J, Sierro F, Peypouquet J, Grimalt J, Shackleton N: Synchroneity between marine and terrestrial responses to millennial scale climatic variability during the last glacial period in the Mediterranean region. Clim Dyn. 2002, 19 (1): 95-105.

    Google Scholar 

  51. 51.

    Moreno A, Cacho I, Canals M, Grimalt JO, Sánchez-Goñi MF, Shackleton N, Sierro FJ: Links between marine and atmospheric processes oscillating on a millennial time-scale. A multi-proxy study of the last 50,000 yr from the Alboran Sea (Western Mediterranean Sea). Quat Sci Rev. 2005, 24 (14–15): 1623-1636.

    Google Scholar 

  52. 52.

    García-Ramos G, Kirkpatrick M: Genetic models of adaptation and gene flow in peripheral populations. Evolution. 1997, 51 (1): 21-28.

    Google Scholar 

  53. 53.

    Cassel-Lundhagen A, Tammaru T, Windig JJ, Ryrholm N, Nylin S: Are peripheral populations special? Congruent patterns in two butterfly species. Ecography. 2009, 32 (4): 591-600.

    Google Scholar 

  54. 54.

    Lesica P, Allendorf FW: When are peripheral populations valuable for conservation?. Conserv Biol. 1995, 9 (4): 753-760.

    Google Scholar 

  55. 55.

    Losos JB, Glor RE: Phylogenetic comparative methods and the geography of speciation. Trends Ecol Evol. 2003, 18 (5): 220-227.

    Google Scholar 

  56. 56.

    Petit RJ, Brewer S, Bordács S, Burg K, Cheddadi R, Coart E, Cottrell J, Csaikl UM, van Dam B, Deans JD: Identification of refugia and post-glacial colonisation routes of European white oaks based on chloroplast DNA and fossil pollen evidence. For Ecol Manag. 2002, 156 (1): 49-74.

    Google Scholar 

  57. 57.

    Schmitt T: Molecular biogeography of Europe: pleistocene cycles and postglacial trends. Front Zool. 2007, 4 (11): 1-13.

    Google Scholar 

  58. 58.

    Templeton AR: Nested clade analyses of phylogeographic data: testing hypotheses about gene flow and population history. Mol Ecol. 1998, 7 (4): 381-397.

    PubMed  CAS  Google Scholar 

  59. 59.

    Ronikier M, Schneeweiss GM, Schönswetter P: The extreme disjunction between Beringia and Europe in Ranunculus glacialis s. l. (Ranunculaceae) does not coincide with the deepest genetic split - a story of the importance of temperate mountain ranges in arctic-alpine phylogeography. Mol Ecol. 2012, 21 (22): 5561-5578.

    PubMed  CAS  Google Scholar 

  60. 60.

    Albach DC, Schönswetter P, Tribsch A: Comparative phylogeography of the Veronica alpina complex in Europe and North America. Mol Ecol. 2006, 15 (11): 3269-3286.

    PubMed  CAS  Google Scholar 

  61. 61.

    Schönswetter P, Tribsch A, Barfuss M, Niklfeld H: Several Pleistocene refugia detected in the high alpine plant Phyteuma globulariifolium Sternb. & Hoppe (Campanulaceae) in the European Alps. Mol Ecol. 2002, 11 (12): 2637-2647.

    PubMed  Google Scholar 

  62. 62.

    Kerner A: Die natürlichen Floren im Gelände der deutschen Alpen. Die Dtsch Alpen für Einheimische und Fremde Geschildert. 1871, 1: 126-189.

    Google Scholar 

  63. 63.

    Aeschimann D, Rasolofo N, Theurillat JP: Analysis of the flora of the Alps. 1: historical account and biodiversity. Analyse de la flore des Alpes 1: Historique et biodiversité. 2011, 66 (1): 27-55.

    Google Scholar 

  64. 64.

    Pampanini R: Essai Sur La Géographie Botanique Des Alpes Et En Particulier Des Alpes Sud-Orientales. Imprimerie Fragnière frères. 1903

    Google Scholar 

  65. 65.

    Theissinger K, Bálint M, Feldheim KA, Haase P, Johannesen J, Laube I, Pauls SU: Glacial survival and post-glacial recolonization of an arctic-alpine freshwater insect (Arcynopteryx dichroa, Plecoptera, Perlodidae) in Europe. J Biogeogr. 2013, 40 (2): 236-248.

    Google Scholar 

  66. 66.

    Schmitt T, Hewitt GM: Molecular biogeography of the arctic-alpine disjunct burnet moth species Zygaena exulans (Zygaenidae, Lepidoptera) in the Pyrenees and Alps. J Biogeogr. 2004, 31 (6): 885-893.

    Google Scholar 

  67. 67.

    Petit RJ, Kremer A, Wagner DB: Finite island model for organelle and nuclear genes in plants. Heredity. 1993, 71 (6): 630-641.

    Google Scholar 

  68. 68.

    Ennos R, Sinclair W, Hu X, Langdon A: Using organelle markers to elucidate the history, ecology and evolution of plant populations. Systematics association special volume. 1999, 57: 1-19.

    Google Scholar 

  69. 69.

    Soltis DE, Gitzendanner MA, Strenge DD, Soltis PS: Chloroplast DNA intraspecific phylogeography of plants from the Pacific Northwest of North America. Plant Syst Evol. 1997, 206 (1): 353-373.

    Google Scholar 

  70. 70.

    Hansen AK, Escobar LK, Gilbert LE, Jansen RK: Paternal, maternal, and biparental inheritance of the chloroplast genome in Passiflora (Passifloraceae): implications for phylogenetic studies. Am J Bot. 2007, 94 (1): 42-46.

    PubMed  CAS  Google Scholar 

  71. 71.

    Azhagiri AK, Maliga P: Exceptional paternal inheritance of plastids in Arabidopsis suggests that low-frequency leakage of plastids via pollen may be universal in plants. Plant J. 2007, 52 (5): 817-823.

    PubMed  CAS  Google Scholar 

  72. 72.

    Drummond A, Ashton B, Buxton S, Cheung M, Cooper A, Heled J, Kearse M, Moir R, Stones-Havas S, Sturrock S: Geneious v5. 1. 2010,,

    Google Scholar 

  73. 73.

    Katoh K, Misawa K, Kuma KI, Miyata T: MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 2002, 30 (14): 3059-3066.

    PubMed  CAS  PubMed Central  Google Scholar 

  74. 74.

    Ronquist F, Huelsenbeck JP: MrBayes 3: bayesian phylogenetic inference under mixed models. Bioinformatics. 2003, 19 (12): 1572-1574.

    PubMed  CAS  Google Scholar 

  75. 75.

    Guindon S, Gascuel O: A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. Syst Biol. 2003, 52 (5): 696-704.

    PubMed  Google Scholar 

  76. 76.

    Goloboff PA, Farris JS, Nixon KC: TNT, a free program for phylogenetic analysis. Cladistics. 2008, 24 (5): 774-786.

    Google Scholar 

  77. 77.

    Drummond AJ, Rambaut A: BEAST: bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007, 7 (1): 214-

    PubMed  PubMed Central  Google Scholar 

  78. 78.

    Fernández-Mazuecos M, Vargas P: Congruence between distribution modelling and phylogeographical analyses reveals Quaternary survival of a toadflax species (Linaria elegans) in oceanic climate areas of a mountain ring range. New Phytol. 2013, 198: 1274-1289.

    PubMed  Google Scholar 

  79. 79.

    Bell CD, Soltis DE, Soltis PS: The age and diversification of the angiosperms re-revisited. Am J Bot. 2010, 97 (8): 1296-1303.

    PubMed  Google Scholar 

  80. 80.

    Fernández-Mazuecos M, Vargas P: Historical isolation versus recent long-distance connections between Europe and Africa in bifid toadflaxes (Linaria sect. Versicolores). PLoS One. 2011, 6 (7): e22234-

    PubMed  PubMed Central  Google Scholar 

  81. 81.

    Lemey P, Rambaut A, Drummond AJ, Suchard MA: Bayesian phylogeography finds its roots. PLoS Comput Biol. 2009, 5 (9): e1000520-

    PubMed  PubMed Central  Google Scholar 

  82. 82.

    Yu Y, Harris A, He X: Rasp (Reconstruct Ancestral State In Phylogenies) 2.0 Beta. 2011,,

    Google Scholar 

  83. 83.

    Ree RH, Smith SA: Maximum likelihood inference of geographic range evolution by dispersal, local extinction, and cladogenesis. Syst Biol. 2008, 57 (1): 4-14.

    PubMed  Google Scholar 

  84. 84.

    Clement M, Posada D, Crandall KA: TCS: a computer program to estimate gene genealogies. Mol Ecol. 2000, 9 (10): 1657-1659.

    PubMed  CAS  Google Scholar 

  85. 85.

    Kelchner SA: The evolution of non-coding chloroplast DNA and its application in plant systematics. Ann Mo Bot Gard. 2000, 87 (4): 482-498.

    Google Scholar 

  86. 86.

    Souquet P, Bilotte M, Canerot J, Debroas E, Peybernés B, Rey J: Nouvelle interprétation de la structure des Pyrénées. C R Acad Sci Paris. 1975, 281: 609-612.

    Google Scholar 

  87. 87.

    Souquet P, Mediavilla F: Nouvelle hypothèse sur la formation des Pyrénées. CR Acad Sci. 1976, 282: 2139-2142.

    Google Scholar 

  88. 88.

    Rivas Martínez S, Báscones Carretero JC, Díaz González TE, Fernández González F, Loidi Arregui J: Vegetación del Pirineo occidental y Navarra: VI Excursión Internacional de Fitosociología (AEFA). Itinera Geobotanica. 1991, 5: 5-456.

    Google Scholar 

  89. 89.

    Montserrat P: L’exploration floristique des Pyrénées occidentales. (The floristic exploration of West Pyrenees). Bol Soc Brot. 1974, 47: 227-241.

    Google Scholar 

  90. 90.

    Villar L: La vegetación del Pirineo occidental: estudio de Geobotánica Ecológica. Príncipe de Viana Suplemento de Ciencias. 1982, 263-434.

    Google Scholar 

  91. 91.

    Izard M: Le Climat. La Végétation des Pyrénees. Edited by: Dupias G. 1985, Paris: Editions du CNRS, 17-36.

    Google Scholar 

  92. 92.

    Librado P, Rozas J: DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009, 25 (11): 1451-1452.

    PubMed  CAS  Google Scholar 

  93. 93.

    Excoffier L, Smouse PE, Quattro JM: Analysis of molecular variance inferred from metric distances among DNA haplotypes: Application to human mitochondrial DNA restriction data. Genetics. 1992, 131 (2): 479-491.

    PubMed  CAS  PubMed Central  Google Scholar 

  94. 94.

    Excoffier L, Laval G, Schneider S: Arlequin (version 3.0): an integrated software package for population genetics data analysis. Evol Bioinformatics Online. 2005, 1: 47-

    CAS  Google Scholar 

  95. 95.

    Corander J, Marttinen P, Sirén J, Tang J: Enhanced Bayesian modelling in BAPS software for learning genetic structures of populations. BMC bioinformatics. 2008, 9: 539-

    PubMed  PubMed Central  Google Scholar 

  96. 96.

    Corander J, Tang J: Bayesian analysis of population structure based on linked molecular information. Math Biosci. 2007, 205 (1): 19-31.

    PubMed  Google Scholar 

  97. 97.

    Dupanloup I, Schneider S, Excoffier L: A simulated annealing approach to define the genetic structure of populations. Mol Ecol. 2002, 11 (12): 2571-2581.

    PubMed  CAS  Google Scholar 

  98. 98.

    Tajima F: Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989, 123 (3): 585-595.

    PubMed  CAS  PubMed Central  Google Scholar 

  99. 99.

    Fu YX: Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997, 147 (2): 915-925.

    PubMed  CAS  PubMed Central  Google Scholar 

  100. 100.

    Slatkin M, Hudson RR: Pairwise comparisons of mitochondrial DNA sequences in stable and exponentially growing populations. Genetics. 1991, 129 (2): 555-562.

    PubMed  CAS  PubMed Central  Google Scholar 

  101. 101.

    Harpending HC: Signature of ancient population growth in a low-resolution mitochondrial DNA mismatch distribution. Hum Biol. 1994, 66 (4): 591-600.

    PubMed  CAS  Google Scholar 

  102. 102.

    Excoffier L, Lischer HEL: Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010, 10 (3): 564-567.

    PubMed  Google Scholar 

Download references


We would like to thank to Jose L. Blanco Pastor and Mario Fernández-Mazuecos for helping in the analysis and the manuscript writing; Emilio Cano, for laboratory assistance and Jose Luis Benito and Virginia Valcárcel for help with the sample collection. We would like to thank Betina Perales for assistance with map creation by using ArcGIS; and Dr. Jukka Corander and Dr. Andrew Rambaut for help with the BAPS and DPA analyses performance respectively. This research was supported by a project “Evolution of the personate flower” (CGL2009-10031) of the Spanish Ministry of Science and Technology and by the Synthesys program (CS and MB).

Data Accessibility

Genbank Accession numbers for DNA sequences: trnS-trnG (KJ911400-KJ911861) and trnK-matK (KJ911294-KJ911399).

Author information



Corresponding author

Correspondence to Isabel M Liberal.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

IML generated the molecular data, performed the analyses and drafted the manuscript. PV conceived of the study, and participated in its design and coordination. Both authors participated in the interpretation of the data and co-wrote the manuscript. MB and CS contributes to the generation of molecular data. CT & MB contributed ideas and suggestions for improving the final manuscript. All authors participated in the obtaining of plant material. All authors read and approved the final manuscript.

Electronic supplementary material

Additional file 1: Comparison of the different taxonomic treatments of two studied taxa according to different authors.(XLS 26 KB)

Additional file 2: Phylogeographical patterns of alpine and montane angiosperms distributed in the Pyrenees. Review based on geographical origin and population dynamics of species during the Quaternary.(DOCX 21 KB)

Samples used for trnS-trnG and trnK-matK sequencing.

Additional file 3: Population numbers after species names in brackets. IML, Isabel Liberal’s collection numbers; EDB, Laboratoire Evolution et Diversité Biologique. For others Voucher see Vargas et al. [22]. Asterisks (*) after GenBank accession numbers refer to those from Vargas et al. [22]. Nomenclature from Vargas et al. [22] is followed except for A. controversum (=A. barrelieri). The samples here named A. majus subsp. striatum (populations 1-3) were treated like A. latifolium in Vargas et al. [22]. Abbreviations used: DEC = dispersal-extinction-cladogenesis analysis; S-DIVAS = Statistical Dispersal-Vicariance Analysis; DPA = Discrete Phylogeographic Analysis; SE = South-east; SW = South-west; NE = North-east; NW = North-west; Pyr. = Pyrenees; Pre-Pyr. = Pre-Pyrenees; Pyr. & adj. = Pyrenees and adjacent areas; S. French basin = Southern French basin; E. Pyrenees = Eastern Pyrenees; C. Pyrenees = Central Pyrenees; W. Pyrenees = Western Pyrenees. (XLS 78 KB)

Additional file 4: Bayesian phylogenetic tree constructed with plastid trnS-trnG sequences (Antirrhinum matrix).(PDF 37 KB)

Additional file 5: Bayesian phylogenetic tree constructed with plastid trnK-matK sequences (Antirrhinum matrix).(PDF 40 KB)

Additional file 6: Results of the Bayes Factor (BF) test for Antirrhinum matrix. Are indicated the Bayes Factors with well-supported rates of dispersal (values of BF > 3). (XLS 27 KB)

Additional file 7: Ancestral area reconstruction (Antirrhinum matrix) using DIVA method as implemented in the software RASP (ex S-DIVA; Yu et al. (2010)). At each node, the most likely inferred ancestral areas are drawn. The four areas are as follow: Iberian Peninsula (A), Pyrenees (B) Mediterranean Basin (C) and South-western Alps (D). Vicariant event among Alps (D) and the Iberian Peninsula (A) is indicated by * and node support. (PDF 44 KB)

DEC (Lagrange) Inference (

Additional file 8: Antirrhinum matrix). The optimal area reconstruction on the branch is represented by a two-letter code (i.e. A|A), being the area on the left the one inherited by the upper daughter branch, and the area on the right the one inherited by the lower daughter branch. (PDF 13 KB)

Haplotype network of the

Additional file 9: Antirrhinum matrix. For geographic abbreviations see [22]. (PDF 32 KB)

(correspondence with Figure

Additional file 10:  3). Distribution of the lineages and haplotypes on the grid (Pyrenees matrix). For each taxon, the first column indicate the number of the haplotypes and number of samples showing that haplotype (in brackets), while the second column indicates the lineage corresponding to each haplotype. The cells shaded in grey show areas harboring at least two lineages. (XLS 34 KB)

Additional file 11: Results of SAMOVA analyses (Pyrenees matrix).(PDF 26 KB)

Bayesian phylogenetic tree constructed with nrITS sequences published in

Additional file 12: 22 plus the five new nrITS sequences corresponding to the samples collected across South-western Alps for this study.(PDF 26 KB)

Haplotypes indentified in the interspecific crossing study.

Additional file 13: Abbreviations used: maj = A. majus; con = A. controversum; cha = A. charidemi; moll = A. mollissimum. (XLS 34 KB)

Additional file 14: Methods.(DOC 80 KB)

Authors’ original submitted files for images

Rights and permissions

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Liberal, I.M., Burrus, M., Suchet, C. et al. The evolutionary history of Antirrhinumin the Pyrenees inferred from phylogeographic analyses. BMC Evol Biol 14, 146 (2014).

Download citation


  • Antirrhinum
  • Phylogeny
  • Phylogeography
  • Pyrenees
  • Quaternary