A glacial survivor of the alpine Mediterranean region: phylogenetic and phylogeographic insights into Silene ciliata Pourr. (Caryophyllaceae)

Studied species

Silene ciliata Pourr. (subsect. Fruticulosae, Caryophyllaceae) is endemic to Europe and inhabits the main Mediterranean mountain ranges in the northern half of Mediterranean Basin countries spreading along the Iberian Peninsula, the Central Massif, the Apennines and the Balkan Peninsula (Tutin et al., 1993). It is an alpine, chamaephytic, perennial, cushion plant, which typically forms pulviniform rosettes of up to 2 cm in height and 15 cm in diameter with high variability in size. Each plant has an average of 13 ± 11 cm (mean ± SD) flowering stems that reach 15 cm in height and bear 1–5 flowers (Giménez-Benavides, Escudero & Iriondo, 2007).

Taxon selection

Twenty-five specimens of S. ciliata populations covering the species distribution range were sampled for this study (Fig. 1). Plant material was obtained from herbarium specimens or directly from the field and stored as silica gel-dried material (Table 1). All field studies carried out by the authors were conducted with the permission of “Junta de Castilla y León” and “Comunidad de Madrid” (approval code numbers: 20144360000894 and 10/117476.9/14, respectively). To assess possible intrapopulation cpDNA variation, DNA from four additional individuals of the Cen3 population was also extracted and amplified.

For the estimation of the phylogeny of the polymorphic cpDNA region, eight additional species of genus Silene, phylogenetically close to Silene ciliata, were included in the study. The selection of these species was based on the most recent phylogenetic studies of Sileneae (Sloan et al., 2009; Greenberg & Donoghue, 2011) and the availability of the required polymorphic cpDNA regions. The search was performed in the GenBank sequence database, and the species selected as outgroups were S. latifolia Poiret, S. uniflora Roth, S. vulgaris (Moench) Garcke and -phylogenetically closer to S. ciliata—S. acaulis (L.) Jacq, S. otites (L.) Wibel, S. nutans L., S. paradoxa L. and S. schafta S. G. Gmel. ex Hohen. Although S. borderei Jordan and S. legionensis Lag. are classified in the same section as S. ciliata in recent Floras (e.g., Flora Europaea, Tutin et al., 1993; Flora Iberica; Castroviejo et al., 1986–2001), we did not consider them because they are not included in the phylogenetic study of Sileneae by Greenberg & Donoghue (2011). The accession numbers of all outgroup regions are listed in Table S3.

DNA markers

In this phylogenetic and phylogeographic approach, we chose maternally inherited plastid DNA (hereafter cpDNA) as the focus of our study, since it provides a conservative and enduring record of plant migrational spread (McCauley, 1997; Irwin, 2002) compared to biparentally inherited nuclear markers that show recombination (Petit, Kremer & Wagner, 1993; Heuertz et al., 2004).

Each of the 25 extracted DNA samples was amplified for the rbcL, rps16 and trnL polymorphic cpDNA regions. These regions were selected from the 12 regions that had previously shown major variation and the best amplification profile (Shaw et al., 2005; Shaw et al., 2007).

DNA extraction and amplification

For DNA extraction, approximately 20 mg of dried leaf tissue of each plant sample were weighed. Extractions were performed following the protocol of Qiagen Plant DNA extraction kit (QIAGEN Inc., CA, USA) with some modifications. The primers used and the PCR conditions applied for each marker, as well as the primer sequences and references, are listed in Table S2. The PCR mix was prepared using PureTaq Ready-To-Go PCR beads (GE Healthcare, Uppsala, Sweden).

Data analyses

Sequencing results were evaluated and corrected manually before being subjected to multiple alignment. The manual corrections were made to check whether the differences found among some bases of the sequences were unique/repeated in some of the sequences and to ensure the presence of gaps. Contigs were assembled and edited with Sequencher 4.1.4 (Gene Codes Corp., MI, USA), Bioedit (Hall, 1999) and ClustalW (Thompson, Higgins & Gibson, 1994). In the latter, default settings were used.

The number of variation and informative sites of our aligned sequences was determined using DnaSP v.5.10.01 (Librado & Rozas, 2009). The phylogenetic analyses were performed using two different statistical approaches (“Bayesian inference” and “Maximum likelihood”) for verification reasons. In the Bayesian analysis, sequence data were first introduced to jModeltest (Posada, 2008) to determine the best fitting evolutionary model according to the AIC criterion. This process was followed to generate a dendrogram for each polymorphic cpDNA region, plus one dendrogram that included all polymorphic cpDNA regions together. The suggested model for rbcL was [HKY], for rps16 [GTR + G], for trnL [HKY + I] and for the tree including all markers [GTR + G]. These models were then inserted into MrBayes 3.1.2 (Huelsenbeck et al., 2001) and posterior probabilities (hereafter PP) were estimated using the Markov chain Monte Carlo (MCMC) method. Four Markov chains were run in parallel for 10,000,000 generations and sampled every 100 generations. The first 100 generations were set as the “burn-in” period, while the rest were used to calculate the 50% majority rule consensus phylogeny and posterior probability. The resulting dendrogram archives were revised with FigTree v. 1.3.1 (Rambaut, 2009). A maximum likelihood dendrogram including all the polymorphic cpDNA regions together was also generated with PhyML 3.0 (Guindon et al., 2010) under the same evolutionary model used for the Bayesian analysis. The reliability of the branches was calculated through bootstrapping, after producing 1,000 bootstrapped data sets. All outputs were compared and analysed to infer the evolutionary history of our study species.

Next, each group of polymorphic cpDNA region sequences was analysed with TCS 1.2.1 (Clement, Posada & Crandall, 2000) and classified according to statistically parsimonious haplotype groups. The haplotype groups were linked by the program, constructing a network of mutation steps, which visualized the genetic distance between them. For the construction of the haplotype networks, deletions were not treated as polymorphic sites, while the analysis was performed under the default of 95% connection limit. To facilitate interpretation, a total cpDNA haplotype network was created with this method. In addition, haplotype networks of rbcL, rps16 and trnL regions were obtained separately to check the congruence between markers. Likewise, total neighbour-net analysis network including all three cpDNA regions together was also designed using Splits Tree v. 4.13.1 (Huson & Bryant, 2006) and following the uncorrected p-distance between individuals. The support for each branch was tested using the bootstrapping method with 1,000 replicates. One final test was performed with the Bayesian Analysis of Population Structure 6 (BAPS, Corander et al., 2008). This Bayesian approach is conditioned on the geographical sampling information available. The actual analysis is performed using a systematic hierarchical Bayesian approach, where a Markov chain Monte Carlo (MCMC) estimation is used whenever the number of possible partitions is too large to be handled with exact calculations. We chose BAPS software to infer the best genetic structure, considering the coordinates of each sample, and ran a test of spatial clustering of individuals, with five replicates for each possible number of groups (K).

Chloroplast haplotype and intrapopulation variation

After multiple alignment evaluation of the three polymorphic cpDNA regions, the final length of the study region resulted in 564 nucleotides for rbcL, 756 nucleotides for rps16 and 509 nucleotides for trnL. Thus, the length of the combined matrix of an “all-marker” region was 1,829 nucleotides. The number of variable sites among chloroplast markers ranged from 4 to 25, while that of parsimony informative sites ranged from 3 to 16 (Table 2). Sequences were submitted to GenBank (accession numbers are available in Table S4).

The intrapopulation study showed no divergence for rbcL and inconsistent polymorphisms (i.e., only present in one individual and probably associated to sequencing errors) in one and two bases for rps16 and trnL, respectively. Therefore, we considered that the evidence for intrapopulation variation was not strong enough to require further testing.

Phylogeny, genetic distance analyses and population structure

No incongruence in results was found among the single markers. Therefore, we essentially used a combined study because all the markers are in the chloroplast genome. The resulting “all-marker” dendrogram from the Bayesian analysis (Fig. 2) revealed two distinct groups, one including all individuals in the western region (i.e., the Iberian Peninsula and France) and another one including all individuals in the eastern region (i.e., the Italian and Balkan Peninsulas). However, the calculated 65% PP for the “eastern group” did not provide a significant difference between the two groups. On the other hand, significant differentiation (100% PP) was found between S. ciliata individuals and the outgroups. Strikingly, two S. ciliata individuals, Pyr1 and Pyr4, were located between the outgroups and the rest of S. ciliata, and were significantly different from them as well as from each other. Both Pyr1 and Pyr4 branches were long, implying high substitution rates. One overarching clade was observed (99% PP) in the “eastern group”, and the Din population was the only one branching off this clade. The “western group” consisted of one clade (78% PP), but also had many separate individual branches. The maximum likelihood dendrogram obtained with the bootstrapping method did not differ, either in formation or in significance of branches support, from the Bayesian dendrogram.

In the haplotype network approach, the overall analysis (Fig. 3) corroborated the existence of two (eastern–western) groups and found higher haplotype diversification in the western group. Of all the eastern populations, the Din haplotype had the nearest position to the western group, which was in agreement with the results obtained from the dendrogram. Sequences assembled into 24 haplotypes, with 12 haplotypes including only western region sequences, 10 haplotypes including only eastern region sequences and two haplotypes outside the network (Pyr1 and Pyr4). Only eastern region sequences Ape1 and Ape2 shared the same haplotype pattern, and no shared haplotype patterns were found between the “eastern” and “western” groups. Moreover, the haplotype network revealed a close relationship between the haplotype pattern of Ape3 and some Balkan populations and among the haplotype pattern of Pyr2 and some Central System populations. The rbcL haplotype network (Fig. S1) was selected to visualise the geographic distribution of haplotypes by regions, as it showed the most representative and parsimonious patterns of the three cpDNA regions, when analysed separately (Fig. 5). In that network, Cen2 and Bal1 haplotypes were prevalent in the western and eastern regions, respectively.

The neighbour-net method suggested a grouping pattern that was in accordance with the one obtained using the haplotype network approach. Besides that, it provided a chance to delve deeper into the differences among S. ciliata sampled populations. The all-inclusive neighbour-net network (Fig. 4) confirmed the classification of all studied populations into a western and an eastern region, which was 91.6% statistically supported. Furthermore, some distances inside the network were noteworthy because they verified previous results. This is the case of the observed 98.4% difference in the distance between Cen1 and the rest of Central System populations (implied by the haplotype network). The Italian Ape3 showed a minor differentiation (77.1%) that was also noticed in the haplotype network and in the dendrograms. Last but not least, the eastern population Din was placed in the “western” group.

The Bayesian spatial clustering of populations resulted in an optimal grouping of K = 2. This supported the western–eastern region division of populations noted in previous analyses. Only the Balkan population Din deviated once more from this division, clustering with the western-region populations.

Genetic diversity in the cpDNA of S. ciliata: a comparative approach

This study reveals high haplotype variability and therefore supports the hypothesis of high cpDNA diversification among S. ciliata populations. Similar results have been reported in previous studies on other Silene species, such as S. latifolia (Ingvarsson & Taylor, 2002), S. vulgaris (Štorchová & Olson, 2004) and S. dioica (Hathaway, Malm & Prentice, 2009), among others. Yet, S. ciliata is ranked among the most varied. Low levels of cpDNA diversification and no diversification at all have been found in S. hifacensis (Prentice et al., 2003) and S. sennenii (López-Vinyallonga et al., 2012), respectively, although this may be due to their narrow distributions and low population sizes (Gitzendanner & Soltis, 2000; López-Pujol et al., 2009; López-Vinyallonga et al., 2012). Considering this, we suggest that the variation detected in S. ciliata is the outcome of a relatively ancient, wider distribution range, followed by a gradual splintering caused by a series of ice ages, as proposed for many other high-elevation species (reviewed by Nieto Feliner, 2014). An alternative explanation to vicariance is that diversification patterns resulted from dispersal. However, long-distance dispersal events most likely played a minor role given that S. ciliata seeds lack any specialized dispersal structure and, thus, most seeds are barochorously dispersed at distances of less than 1 m (Lara-Romero et al., 2014). Finally, although the step node of the species is dated around 10 million years ago (Sloan et al., 2009), and the diversification of the southern European alpine flora has been correlated with the climatic oscillations of the Pleistocene (Kadereit, Griebeler & Comes, 2004), we have no information about the date of the crown node. Therefore, we have no evidence of the period when the diversification into the two subspecies took place.

Interpreting the distinction of S. ciliata between western and eastern regions and their origin

No evidence was found against the classification of S. ciliata into a western and an eastern race (Blackburn, 1933; Tutin et al., 1993). Hence, we propose maintaining the names Silene ciliata subsp. ciliata and S. ciliata subsp. graefferi to describe the noted clustering of S. ciliata individuals into a western and eastern group, respectively. On the other hand, both dendrograms indicated a significant difference between S. ciliata individuals and the outgroups which, together with the nonessential divergence between populations, corroborates the monophyly of our species.

Tracing back to the species’ differentiation, we hypothesize that populations of an ancestor of S. ciliata dominated the Mediterranean Basin. At the onset of glacial period climatic oscillations in the late Tertiary and in the Quaternary period, these ancestral populations might have been forced to migrate to favourable areas, while those unable to encounter a glacial refugium because of distance, time or natural barriers perished. Given that we are dealing with an alpine species, S. ciliata populations could have migrated following the paths that constitute links between neighbouring mountains. The Alps mountain range system seems to have posed a persistent and significant hurdle for this species’ migration. A rigorous example supporting this theory is that during Quaternary glaciations, the Alps, in contrast to the Mediterranean mountains, were extensively covered with ice sheets (Hughes, Woodward & Gibbard, 2006; but see Stehlik et al., 2002). This is in accordance with previous phylogeographic studies (e.g., Taberlet et al., 1998; Hewitt, 2000) and may explain why S. ciliata populations have not been found there. Moreover, it would account for the observed disconnected distribution and division of the species into the western and eastern groups, since the geographical borders formed by the two groups coincide with the location of the Alps. A similar grouping pattern has been found in the Mediterranean for Androsace vitaliana (Vargas, 2003) and Heliosperma (Frajman & Oxelman, 2007), genera with the barrier shifting west and east of the Alps region, respectively. Disjunction in distribution, possibly resulting from the Alps and distinction into two subspecies has recently been proposed in Artemisia eriantha, another alpine plant distributed along the Alps and many Mediterranean mountains (Sanz et al., 2014), which also supports our hypothesis.

Evolutionary processes and geo-climatic effects on western and eastern populations

Apart from the significant difference found between eastern and western cpDNA sequences, further important diversification was found inside each group. Polyploidization is an evolutionary mechanism that has generated evolutionary lineages during the Pleistocene (Stebbins, 1984). The number of chromosomes of S. ciliata is mostly n = 24 or n = 48, although populations with much higher levels of ploidy have been described in S. ciliata subsp. ciliata (Blackburn, 1933; Küpfer, 1974; Tutin et al., 1993). Hence, we propose that this mechanism could explain some of the differentiation within S. ciliata. Nevertheless, considering available chromosome data, we found no relationship between chromosome number and the clustering obtained in our results. Variation may have also resulted from slow mutation events inside disjunct refugia during periods of adverse climatic conditions (Sanz et al., 2014) either as an outcome of elevational range shifts (lowland glacial refugia) or in situ endurance (nunataks). Additionally, other sources of diversification like genetic drift associated to low population sizes and prolonged isolation should be considered (Young, Boyle & Brown, 1996).

Regarding the western group, genetic diversity is apparent in the Pyrenees mountain range and has led to the genetic disaffiliation of the range into a western and an eastern section. This is in line with the genetic break found in Artemisia eriantha (Sanz et al., 2014). Another component of the western group diversification was introduced by the highly divergent Cen1 sequence of Serra da Estrela, suggesting high isolation of this population. On the other hand, the merging of Pyr2 sequence with Central System S. ciliata individuals may imply braided migrational paths between these populations during glacial-interglacial events.

Interestingly, the degree of divergence recorded in the eastern group of S. ciliata is higher than that in the western group. This observation has also been made for temperate trees and shrub taxa (Petit et al., 2003). This higher genetic diversity and the existence of more unique haplotypes, especially in the Balkan Peninsula, might be due to the additional effect of the complex orography and restricted territorial extent of existing refugia, which did not facilitate exchange among populations. More specifically, the various orientations of mountain chains in the Balkans may have acted as a barrier to migration (Tzedakis, 2004). The individuals from the western part of the eastern groups (e.g., Ari and Din) showed some important differences in certain analyses (see Figs. 3 and 4). This might be related to the nature of the east Balkan slopes, which have a more gentle relief compared to the steep west mountains (Reed, Kryštufek & Eastwood, 2004), thereby fostering higher levels of isolation. Further differentiation in Din could be because the Dinaric Alps were much less affected by glaciations than the rest of the Mediterranean mountain systems (Frajman & Oxelman, 2007), resulting in the maintenance of relict populations. Lastly, the close relationship of Italian Ape3 with some Balkan populations (especially Bal4) might result from the proposed land connection of the north Italian and Balkan Peninsulas during the early Holocene (approx. 20–16 ka BP) (Lambeck et al., 2004) which would have facilitated dispersal events between the two regions.

The Pyrenees case

The Bayesian and maximum likelihood analyses showed that Pyr1 and Pyr4 differed from the outgroups as well as from the rest of S. ciliata individuals and were situated in an intermediate position between them in the dendrogram (see Fig. 2). Similar results were found in the rest of the analyses. We surmised that this pattern could be another example of the Pyrenees range acting as a stable hybrid zone, as argued in Chorthipopus parallelus (Hewitt, 1993) and Saxifraga subsect. Triplinervium (Mas de Xaxars et al., 2015). At any rate, the rise of hybrid zones due to glaciations, and hence, the preservation of different species genomic information via gene flow (Harrison, 1990) are linked with high altitudes (Hewitt, 2001 and references therein). In the case of Pyr1 and Pyr4, their haplotype patterns may have resulted from chloroplast capture between S. ciliata and other congeneric, sympatric species (Rieseberg & Soltis, 1991). After all, the geographical contact of congeneric species causing chloroplast sharing has been reported in other studies including species in the same genus, like S. latifolia and S. dioica (Prentice, Malm & Hathaway, 2008), as well as in other plant groups (e.g., Gardner et al., 2004; Okuyama et al., 2005). Given the number of haplotype patterns detected inside the species, the alternative explanation of Pyr1 and Pyr4 resulting from lineage sorting is another option that cannot be readily rejected (Galtier & Daubin, 2008).

Conclusions and future prospects

Our results confirm the monophyly of S. ciliata due to the differences found between the studied populations and the outgroups and reveal a clear west-to-east division of S. ciliata populations with the borderline set in the region of the Alps. This division validates the past classification of the species into two subspecies; S. ciliata subsp. ciliata found west of the Alps (“Spanish race”, Blackburn, 1933) and S. ciliata subsp. graefferi located east of the Alps (“Italian race”, Blackburn, 1933). Major intraspecific variation is supported by all analyses, but none of them supports the occurrence of additional varieties or subspecies (according to Küpfer, 1974; Castroviejo et al., 1986–2001). In addition, we suggest that geographic and climatic factors may have played a central role in the evolutionary history of the species and the formation of the two subspecies. Further analyses including more individuals and cpDNA markers, as well as mitochondrial DNA (mtDNA) markers and nuclear ribosomal internal transcribed spacer (nrITS) regions, are encouraged to secure conclusions of this role and clarify the status of unsolved-incongruent populations. Molecular clocks, the inclusion of additional congeneric species and increased sampling effort are necessary to resolve the remaining questions.