Genetic structure and phylogeography of Juniperus phoenicea complex throughout Mediterranean and Macaronesian regions: different stories in one

Abstract

Key message

The genetic structure of Juniperus phoenicea in the Mediterranean Basin is inferred using amplified fragment length polymorphism markers (AFLP) markers. As other Mediterranean conifers, J. phoenicea populations show moderate levels of genetic diversity and interpopulational differentiation. The pattern of distribution of genetic diversity seems highly influenced by the climatic fluctuations which occurred in the Pleistocene.

Context

It has been stated that the genetic structure of Mediterranean conifers is mediated by the historical climatic changes and the geological rearrangements which occurred in the Mediterranean Basin. J. phoenicea provides an excellent example to test how its genetic structure is influenced by these events.

Aims

In this work, we study the amount and distribution of genetic diversity of J. phoenicea complex, in order to evaluate its taxonomic status and to reveal underlying phylogeographic patterns.

Methods

The molecular diversity was analyzed for 805 individuals from 46 populations throughout its distribution range using AFLP markers. Principal coordinate analysis, analysis of molecular variance (AMOVA), and Bayesian-based analysis were applied to examine the population structure.

Results

AFLP markers revealed moderate levels of intrapopulation genetic diversity, pairwise genetic differentiation, and a clear pattern of isolation by distance. Bayesian analysis of population structure showed five clusters related to the taxonomic status of J. phoenicea and J. turbinata, and a geographic pattern of genetic structure in J. turbinata.

Conclusion

All the analysis separate J. phoenicea from J. turbinata. For J. turbinata, up to four groups can be distinguished from a phylogeographic point of view. The genetic structure of J. turbinata seems highly influenced by climatic and geologic fluctuations occurring since the Oligocene.

1 Introduction

The Mediterranean Basin is considered a biodiversity “hot spot” (Médail and Quézel 1997), because it contains approximately 10% of all higher plant species on Earth, despite only covering an area equivalent to less than 2% of the world’s land masses (Fady-Welterlen 2005). The main causes of this species richness are the great diversity of habitats that it shows (Blondel and Aronson 1999; Thompson 2005) arising from the drastic climatic and geologic changes that have occurred since the Tertiary, which has led to population fragmentation, the creation of glacial refugia for Tertiary species and subsequent speciation, and the creation of a great number of endemic species (Médail and Diadema 2009). The diversity of Mediterranean woody species is also very high, representing three times the number of European temperate tree species (Fady-Welterlen 2005). With this background in mind, a study of the genetic diversity of Mediterranean gymnosperms is of great interest, because it provides and explanation of the evolutionary strategies followed until now. Unlike conifers in general, Mediterranean conifers show a much larger span of gene diversity and genetic differentiation (Fady-Welterlen 2005). Furthermore, they exhibit some characteristics that provide a great opportunity for understanding how geological and climate changes might have influenced the evolution of different lineages, because they are long-lived outcrossing species that arose mostly in ancient times and dominated the surface of the Earth before the last palaeoclimate changes that occurred in the Mediterranean Basin.

Juniperus grex phoenicea (family Cupressaceae) is a small dioecious or rarely monoecious Mediterranean tree whose distribution encompasses the Mediterranean Basin and Macaronesian regions, from Madeira and Canary Islands in the west to Jordan and Saudi Arabia in the east, although it is most prevalent throughout the Iberian Peninsula and North Africa (Browicz and Zieliński 1982; Christensen 1997; do Amaral Franco 1986; Farjon 2005; Greuter et al. 1984; Jalas and Suominen 1973; Quézel and Pesson 1980). The species comprises two subspecies, J. phoenicea L. subsp. phoenicea and J. phoenicea subsp. turbinata (Guss) Nyman, which were defined based on taxonomic (Farjon 2005), morphological (Mazur et al. 2010; Mazur et al. 2016), and molecular studies (Adams et al. 2002; Boratyński et al. 2009; Dzialuk et al. 2011). Alternative interpretations define two independent species, J. phoenicea s. str. and J. turbinata Guss., with J. turbinata being discussed by several authors (see Farjon 2005), and its recognition as a species supported by recent molecular studies (Adams 2014; Adams et al. 2013; Adams and Schwarzbach 2013). As such, J. phoenicea ≡ J. phoenicea subsp. phoenicea is generally considered as being endemic to the Iberian Peninsula, south of France and Norhtwestern Italy, while J. turbinata ≡ J. phoenicea subsp. turbinata is widespread throughout the Mediterranean Basin and Macaronesia. However, Macaronesian populations were formerly described as J. canariensis Guyot, but were later combined as a new subspecies, J. turbinata subsp. canariensis by Rivas-Martínez et al. (1993), with the subspecific range being widely accepted in recent years. However, its taxonomic status has recently been challenged based on several molecular markers (RAPDs, Adams et al. 2006; cpDNA, Adams et al. 2010; nrDNA, and cpDNA, Adams et al. 2013), leading to the conclusion that J. turbinata subsp. canariensis should not be supported at the subspecific level.

Although several studies exist about the genetic structure of J. phoenicea based on different molecular markers (e.g., Boratyński et al. 2009; Dzialuk et al. 2011; Meloni et al. 2006), there is no one in which the whole distribution area has been sampled. For this purpose, we have chosen amplified fragment length polymorphism markers (AFLP) (Vos et al. 1995). Since its development, this technique has been applied to address questions regarding genetic relatedness among individuals, population structure, phylogenetic relationships, and mapping of quantitative trait loci (QTL) (Mueller and Wolfenbarger 1999; Meudt and Clarke 2007) The power of AFLP analysis derives from its ability to quickly generate large numbers of marker fragments for any organism, without prior knowledge of genomic sequence. In addition, AFLP requires only small amounts of starting template and it exhibits high reproducibility. Nevertheless, it still suffers from several weaknesses, such as its dominant nature, making difficult distinguish between homozygotes and heterozygotes, and problems of homology assignments of fragments (Koopman 2005), which can cause the variance of estimators, leading to significant biases (Lynch and Milligan 1994; Krauss 2000; Krutovskii et al. 1999; Meudt and Clarke 2007). Anyway, AFLP markers have proven to be very useful addressing the genetic diversity and structure of populations in related species. Therefore, in this work, we have used these molecular markers with the aim to elucidate how the genetic variation of populations is distributed, evaluate the taxonomic status of the two species (J. phoenicea and J. turbinata) and compare our results with previous studies about J. phoenicea and other Mediterranean conifers.

2 Material and methods

This article is a continuation of a partial study of J. turbinata (Jiménez et al. 2017) and includes data previously used (Canarian populations and three Moroccan populations).

2.1 Plant sampling

The sampling strategy covered the whole J. grex phoenicea geographic range, so that material was collected from 46 populations belonging to Mediterranean Basin (38 populations) and the Canary Islands (8 populations). A total of 810 individuals were included in the study. The number of individuals sampled per population varied from 10 to 20, with the exception of one (MURC, Murcia) where only 5 mature individuals could be sampled. The locations of populations are shown in Table 1 and supplementary material Fig. S1. Transects were made collecting specimens sufficiently far apart (wherever possible) to avoid collecting very closely related individuals. Only small branches (under 10 cm) were sampled to avoid damaging specimens. To avoid the degradation of plant tissues, all plant material was labeled and kept in sealed bags with silica gel as recommended by Sytsma et al. (1993), until DNA extraction.

Table 1 Geographical location of 46 sampled populations of Juniperus grex. phoenicea (3 belonging to J. phoenicea and 43 belonging to J. turbinata) and genetic diversity estimates for every population and groups of populations showed by BAPS and Structure.

2.2 DNA analyses

DNA was extracted from leaves using the × 2 cetyltrimethylammonium bromide (CTAB) method (Doyle and Doyle 1987). Total DNA extracts were quantified using a Nanodrop 2000 spectrophotometer. Amplified fragment length polymorphism (AFLP) analysis was carried out following Vos et al. (1995) with minor modifications, using fluorescence-labeled primers (VIC, NED, PET, Applied Biosystems, Madrid, Spain) instead of radioactive labeled primers. Fragments were selectively amplified with the primer pairs EcoR1-ATG/MseI-CGT, EcoRI-ACA/MseI-CGT, and EcoRI-ACG/MseI-CAC. Multiplex products were run for 4 h on an ABI 377 sequencer to separate fragments together with an internal size standard (GeneScan 600 LIZ, ABI). Sizing and peak identification were performed using the Genemapper 4.0 software (Applied Biosystems). The presence of a band was scored as 1, whereas the absence of a band was scored as 0 in the finally assembled binary matrix.

As AFLP markers are dominant, we assumed that null bands were homologous and that populations were in Hardy-Weinberg equilibrium (Lynch and Milligan 1994) in order to compute diversity indices [percentage of polymorphic markers (PLP) and Nei’s (1978) unbiased expected heterozygosity (He)] and genetic distance among populations (Fst). These parameters were calculated for all the populations sampled and for groups of populations depicted by Bayesian methods. Additionally, to eliminate the influence of uneven sample sizes per population, He was calculated for a standardized population size of ten randomly selected individuals, and a Tukey test was performed in order to observe significant differences among He and standardized He. These parameters were inferred with AFLP-surv 1.0 (Vekemans 2002). Significance of Fst values was determined using 1000 bootstrapped data sets. In order to test influence of geography on genetic diversity pairwise correlations between Nei’s heterozygosity (He) and geographic parameters (latitude, longitude, and altitude), were made for Mediterranean populations using Kendall’s τ estimate. The significance of the correlation was tested using a permutation procedure (10,000 permutations) with Past (Hammer et al. 2001). In addition, an ANOVA and post hoc HSD Tukey test were performed with SPSS 19.0, in order to address significant differences in He among the groups of populations detected by Bayesian methods.

Analysis of molecular variance (AMOVA, Excoffier et al. 1992) was conducted to estimate variance components at several hierarchical levels, partitioning the variation among predefined geographic groups (those proposed by BAPS and Structure analyses), among populations, and among individuals within populations. Pairwise Fst values between all populations were also calculated and tested for significance by resampling with 1000 random permutations. Furthermore, these pairwise Fst are shown as a graphic using the R graphical tool implemented in Arlequin 3.5 (Excoffier and Lischer 2010), because it allows genetic affinities between populations to be observed easily.

A Bayesian model-based analysis was performed to infer population structure with Structure version 2.2 (Falush et al. 2007; Pritchard et al. 2000). The F model, based on an admixture ancestry model with correlated allele frequencies, was imposed to estimate the posterior probabilities [LnP(D)] of K groups (Pritchard and Wen 2004) and the individual percentages of membership assigned to them according to their molecular multilocus profiles (Falush et al. 2003; Falush et al. 2007). Probabilities for a range of K were examined starting from one to the number of sampled populations plus one (K = 1–47), using a burn-in period and run length of the Markov chain Monte Carlo (MCMC) of 10⁵ and 10⁶ iterations, respectively, replicated 20 times. Results were uploaded into STRUCTURE HARVESTER (Earl and von Holdt 2012), which estimates the most likely K value (ΔK), following Evanno et al. (2005). We used CLUMPP 1.1.2 (Jakobsson and Rosenberg 2007) to reach a consensus on the results of the independent runs for the optimal K. For the consensus, we used the Greedy option with random input order and 10,000 repeats. The consensus was visualized in DISTRUCT 1.1 (Rosenberg 2004). Bayesian analyses of the genetic structure were also conducted with BAPS 5.4 [Bayesian Analysis of Population Structure, Spatial Clustering of Groups (Corander and Marttinen 2006)], using stochastic optimization instead of MCMC to find the optimal partition. We performed a mixture analysis of individuals without using the geographic origin of the samples as a prior informative (clustering of individuals). BAPS simulations were run with the maximal number of groups (K) set from 1 to 47. Each run was replicated 10 times, and results were averaged according to the resultant likelihood scores. The output of the mixture analyses was used as input for population admixture analysis (Corander and Marttinen 2006), with the default settings in order to detect admixture between clusters.

In addition, a principal coordinate analyses (PCoA) was performed to illustrate overall similarity among individuals using multivariate statistical package software (MVSP version 3.12d). PCoA and a neighbor-joining (NJ) tree with the PHYLIP packages (NEIGHBOR and CONSENSE; Felsenstein 1989) was inferred from the pairwise Nei’s genetic distance (Nei 1978) between all pairs of AFLP phenotypes. Branch support was assessed by bootstrap analysis (1000 replicates).

Mantel tests (Mantel 1967) was carried out to assess correlation among genetic and geographic distances. Thereafter, the fine-scale genetic structure among populations was investigated using spatial autocorrelation analysis in GenAlEx 6.5 (Peakall and Smouse 2006). The spatial autocorrelation coefficient (r) was computed using the geographic distance and the genetic distance between individuals. Tests for statistical significance were conducted by random shuffling (999 times) of individual geographic locations to define the upper and lower bounds of the 95% confidence interval for each distance class, and estimating 95% confidence intervals around mean r values by bootstrapping pairwise comparisons within each distance class (1000 repeats). Analyses were performed using the variable distance classes’ option based on the geographical distances estimated in each population. Mantel test analyses were also undertaken using GenAlEx 6.5.

Data availability

The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

3 Results

3.1 Genetic diversity and differentiation

The three selective AFLP primer combinations amplified a total of 757 fragments. In total 99.7% of the fragments were polymorphic and all 810 individuals had unique AFLP profiles. The proportion of polymorphic loci (PLP) for a single population ranged from 39.6% (MEN1) to 66.3% (CERD). Expected heterozygosity values showed that CRET was the most diverse population (He = 0.228), whereas CAD1 showed the lowest within-population diversity (0.147). The average gene diversity within populations (Hs) was 0.174, whereas the total diversity (Ht) value was 0.226 (Table 1). The fixation index was highly significant (Fst = 0.232, P < 0.000), suggesting moderate to high genetic differentiation among populations. Standardized He values are listed in Table 1. The Tukey test suggests that intrapopulation genetic diversity is not influenced by the number of individuals sampled (T = − 0.429, P = 0.669). There was no correlation between the genetic diversity and the latitude nor between genetic diversity and altitude (τ = − 0.032, P = 0.772; τ = 0.087, P = 0.443), but there was positive correlation with longitude (τ = 0.328, P = 0.003). When populations were grouped according to the clusters obtained with BAPS and Structure, we observed that easternmost populations showed higher Ht and Hs values than westernmost (Table 1). Furthermore, ANOVA test found significant differences in He values among these defined groups of populations (F = 13.191, P = 0.0001). AMOVA analysis structuring populations according to regional differences showed that the greatest percentage of variation was attributable to diversity within populations (63.14%). Among populations within regions accounted for 19.09% of the variation, whereas 17.77% of the variation was attributable to geographical differences (the four groups depicted by the PCoA). The Φst obtained with AMOVA (0.369) was higher than that showed by AFLP-SURV, suggesting high differentiation among populations. The graphic of the pairwise Φst obtained with AMOVA (Fig. 1) showed that J. phoenicea and the Canarian populations were clearly differentiated from the remaining populations, and also the Murcia population was. Pairwise Fst among Canarian populations were the lowest.

Fig. 1

Matrix of pairwise Φst of 46 populations of J. grex phoenicea obtained with Arlequin 3.5. Codes of populations are listed in Table 1

3.2 Identification of genetic structuring

The Bayesian analysis of genetic structure of J. phoenicea populations conducted with Structure found the highest LnP(D) and ΔK values for K = 5 (Fig. 2). The genetic structure revealed by BAPS was highly congruent with that obtained with Structure, identifying five clusters (p = 0.999). However, only 10.87% of individuals were admixed (Fig. 2). Whereas, J. phoenicea populations were assigned to one cluster, the Canarian populations were grouped in the second cluster. The third cluster grouped Algerian populations with high membership and several Mediterranean populations (Ibiza, Menorca, Italy) that were intermingled with the fourth cluster, in which eastern Mediterranean populations were grouped together with Mallorca, Sardinia and, surprisingly, Murcia. Finally, Moroccan and Iberian populations (except for Murcia population) were mainly assigned to the fifth cluster. PCoA analysis gave similar results (Fig. 3). PCoA plot revealed five clusters, one comprising the samples belonging to J. phoenicea, a second cluster that grouped Canarian populations of J. turbinata, a third cluster in which Algerian, Ibiza, Menorca, and Italian populations were grouped, a fourth group that included the remaining easternmost Mediterranean populations together with MALL, CERD, and MURC populations, and a fifth cluster comprising Moroccan and Iberian populations. For this analysis, the first three axes accounted for 43.20% of the variation (18.27, 16.33, and 8.59%, respectively). The NJ tree grouped populations into the former five clusters, which were supported with moderate to high bootstrap values, except for the admixed populations of Ibiza, Menorca, and Italy, that were located as sister clade of Murcia, Mallorca, and Chipre populations (Fig. S2).

Fig. 2

Location map and Bayesian analysis of the genetic structure of 810 individuals from 46 sampled populations of J. grex phoenicea. a Mean estimated log-normal (Ln) probability of the data in relation to the simulated number of clusters K. Vertical bars indicate the standard deviation among 10 replicates. b Delta K in relation to number of clusters K. c Plot showing the membership of individuals for five predefined clusters with c1, Structure and c2, BAPS. Number of populations are listed in Table 1. d Mediterranean populations. e Canarian populations. Rectangular bars represent the mean proportion of membership of each population to the predefined K = 5 clusters with Structure. Dotted lines separate populations belonging to each of the predefined clusters. f Geographic location of 46 sampled populations. Red circles, J. phoenicea, green circles, Mediterranean populations of J. turbinata, blue circles, Canarian populations of J. turbinata

Fig. 3

PCoA plot of 46 populations of J. grex phoenicea based on pairwise Nei’s (1978) genetic distances. Dark gray diamonds, J. phoenicea, white circles, Moroccan and Spanish J. turbinata populations, black circles, Algerian and Central Mediterranean J. turbinata populations, gray squares, Eastern and Central Mediterranean populations, gray triangles, Canarian populations of J. turbinata

Mantel test showed a great isolation by distance pattern (R = 0.295, P = 0.007). Spatial autocorrelation analysis showed positive isolation by distance at distances of 400 km approximately (Fig. 4).

Fig. 4

Correlogram from spatial autocorrelation analysis using the correlation coefficient r, and variable distance classes. 95% confidence error bars for r were estimated by bootstrapping over pairs of samples; dashed lines (U, L) represent upper and lower bounds of a 95% CI for r generated under the null hypothesis of random geographic distribution

4 Discussion

The genetic structure of J. phoenicea has already been partially studied using several molecular markers (ISSR, isozyme, RAPDs). However, an AFLP-based analysis of the genetic structure encompassing its whole distribution area (Macaronesian and Mediterranean populations) has not been performed before. Genetic diversity values obtained with molecular markers are difficult to compare between studies, because the statistical methods used to process the data often vary, particularly with dominant and co-dominant markers. However, parameters such as intrapopulation genetic markers may serve to establish comparisons between species (Nybom 2004). Regarding the He values found in our study, studies on J. grex phoenicea by Boratyński et al. (2009) and Dzialuk et al. (2011) produced similar results to ours, while Meloni et al. (2006) obtained lower Hs values. The intrapopulation genetic diversity observed in our study falls within the range of the values reported by Nybom (2004), or other European conifers (Conord et al. 2012; Müller-Starck et al. 1992). Numerous studies have shown that conifers, as long-lived and outcrossing species, typically exhibit high levels of within-population genetic variation (Dzialuk et al. 2011; Hamrick et al. 1992; Nybom and Bartish 2000). However, it is often difficult to predict variation patterns within the context of the complex historical, geological, and climatic changes that have occurred in the Mediterranean Basin. In fact, our results show moderate to high levels of genetic diversity in spite of the severe fragmentation and isolation of populations (Browicz and Zieliński 1982; Charco 2001; Quézel and Pesson 1980; Quézel et al. 1992). The level of intrapopulation genetic variability in J. phoenicea is thought to be associated with repeated long periods of survival of separate populations within the territories around the refugial areas under diverse environmental conditions, probably in the maritime mountain regions during the Pleistocene. Such levels of genetic diversity have been maintained because population sizes were big enough to avoid the effects of severe bottlenecks (Boratyński et al. 2009). Moreover, the genetic diversity values observed follow a statistically significant east-west gradient of decrease in intrapopulation genetic diversity, as suggested by Fady and Conord (2010) and Conord et al. (2012), because the western Mediterranean climatic conditions during the Quaternary Glaciations were more severe, involving colder and more arid periods than the eastern area (Panagiotopoulos et al. 2014; Sánchez-Goñi et al. 2002; van Andel 2002).

It has frequently been stated that conifer show low levels of genetic differentiation (Hamrick et al. 1992; Lagercrantz and Ryman 1990; Nybom and Bartish 2000). However, AFLP markers used in this study revealed great inter-population genetic differentiation. Our Fst value is similar to those obtained in other studies involving the genetic structure of J. phoenicea, which range from 0.106 to 0.43 (Boratyński et al. 2009; Dzialuk et al. 2011; Jiménez et al. 2017; Meloni et al. 2006). This variation in Fst values might be due to the fact that different populations, number of individual sampled per population, and types of marker have tended to be used in these studies. Nevertheless, interpopulation genetic differentiation is quite common in Mediterranean conifers, e.g., in Abies cilicica (Antoine and Kotschy) Carrière (Sękiewicz et al. 2015) Cupressus sempervirens L. (Korol et al. 1997), J. dupracea M. Bieb. (Sobierajska et al. 2016; J. thurifera L. (Jiménez et al. 2003; Terrab et al. 2008), and Tetraclinis articulata (Vahl) Masters (Sánchez-Gómez et al. 2013). The reason for this wide differentiation must be the fragmentation of populations caused by the profound geologic and climatic changes that occurred throughout the Mediterranean Basin. In spite of the fact that J. phoenicea show long distance dispersal ability (Jordano 1993; Padilla et al. 2012), we have observed genetic differentiation and a significant isolation by distance pattern up to 400 km.

As regards, the distribution of genetic variation, the PCoA showed five groups. The first one included the populations from eastern Spain (COPE, RUBIE, and GONT), which are considered as J. phoenicea s. str. Previous studies have demonstrated that this species, distributed through eastern Spain, south of France, and north of Italy, is clearly distinct morphologically, biochemically, and genetically from J. turbinata (Boratyński et al. 2009; Lebreton and Pérez de Paz 2001; Mazur et al. 2010, 2016; Sánchez-Gómez et al. 1993). One reason for the differentiation between J phoenicea and J. turbinata has been attributed to the colonization of current territories by J. phoenicea’s populations from more glacial refugia, which were isolated during longer periods of time than those from temperate areas (Fady-Welterlen 2005), where populations of J. turbinata also took refuge. Our results agree with those of Boratyński et al. 2009, who suggested isolation of the species could have occurred at the end of the Messinian salinity Crisis (about 5.3 million years ago). Also the phenological differentiation and absence of cross-pollination between species (Arista et al. 1997), could have favored the persistence of genetic differentiation. In this group, the COPE population showed a slight degree of genetic differentiation with respect to RUBIE and GONT. Not only is this population isolated from the others, but its climatic and ecology characteristics are also quite different. It appears at sea level, in more arid and warmer conditions.

The second group showed by PCoA, BAPS, and Structure includes the Canarian populations of J. turbinata, which formed a clearly separated group. As commented before, these populations were previously described as J. canariensis and latter combined as J. phoenicea subsp. canariensis, and, more recently, several authors have led to the conclusion that the subspecific status should not be supported (Adams et al. 2006; Adams et al. 2010; Adams et al. 2013). Results obtained in this work reopen the debate. The level of endemic flora in Canary Islands is more than 50% of native flora (Santos-Guerra 2001), and the main cause of this richness in endemic species is the effective radiation of lineages within the archipelago, which it has been promoted by an important ecological heterogeneity and geographical isolation between islands (Reyes-Betancort et al. 2008 and references therein). However, J. turbinata do not follow this pattern. Jiménez et al. (2017) suggested effective gene flow among islands due to the absence of genetic differentiation in Canarian populations and that these populations would be under anagenetic evolution, i.e., the survival and proliferation of an immigrant population over many generations without splitting events (Stuessy et al. 2006). Pairwise Fst values point to a lower degree of genetic differentiation between the Canarian populations and the rest of J. turbinata populations than that observed between J. phoenicea and J. turbinata, or even among several Mediterranean populations (Table S1). The genetic differentiation between Canarian and the rest of populations could be due to the early colonization of the islands and a later isolation from mainland populations. Furthermore, the historical active vulcanism of the archipelago the well-documented reduction of populations since the arrival of settlers and the effective gene flow among islands (Jiménez et al. 2017), might have played an important role in the genetic differentiation detected.

The third and fourth clusters showed by PCoA grouped Algerian and Eastern Mediterranean populations respectively, whereas Central Mediterranean Populations (MURC, IBIZ, MEN, MALL, CERD, and ITAL) are admixed between these two clusters (Fig. 2). Results depicted by Bayesian structure suggest that European Central Mediterranean populations might have followed a different migratory route from those of North Africa, as in the case of other Tertiary woody species, e.g., Quercus ilex (Lumaret et al. 2002). Subsequent gene flow might have occurred between some Central Mediterranean and Algerian populations. These relationships between Algerian and Central Mediterranean populations have been observed previously by other authors (e.g., Burban and Petit 2003; Linares 2011; Petit et al. 2002; Terrab et al. 2008).

Otherwise, the San Pedro population from Southeastern Spain (MURC), although included in this cluster, is clearly differentiated from the rest of the populations. This population suffered a drastic bottleneck during the sixteenth century, because the mayor of the town ordered every tree and shrub where populations of J. turbinata were located to be felled in order to avoid Berber pirates using these forests as refuge every time they attacked coastal villages (Zamora and Grandal 1997). As a result, only a small population formed by a few individuals survived, and the effects of genetic drift are still remarkable.

The last cluster depicted by PCoA and Structure grouped Moroccan and south-southwestern Iberian populations, all of them considered as J. turbinata. Dzialuk et al. (2011) suggested that genetic divergence between European and Moroccan populations was higher than that between subspecies. They also suggested that one of the reasons might be that the Strait of Gibraltar could have acted as an important barrier to gene flow between Moroccan and European populations, although they observed low Gst values (0.065). However, our results disagree with the above. We have observed low genetic differentiation among the populations included in this cluster (Table 1, Fst = 0.113). We even found several pairs of populations formed by one from the Iberian Peninsula and one from Morocco that showed no genetic differentiation, whereas Fst values between any population belonging to this cluster and J. phoenicea populations (COPE, RMO, and GON) were higher (Table S1), a result previously observed by Boratyński et al. (2009). A possible explanation for these differences could be the molecular marker used in this study and the different sampling strategy. Furthermore, Dzialuk et al. (2011) suggested that other markers showed the differences between J. phoenicea and J. turbinata better than RAPDs. In short, our results agree with Arista et al. (1997), who suggested effective reproductive isolation between subspecies, and effective gene flow (at least until recent times) among populations from Africa and Iberian Peninsula.

Whatever the case, previous works (Boratyński et al. 2009; Fady-Welterlen 2005; Thompson 2005) suggest that the origin of genetic differentiation in Mediterranean population is due to the abovementioned fragmentation which took place during Pleistocene times. However, for Tertiary species ancient palaeogeographic events should be taken into account (Nieto Feliner 2014). The origin of J. phoenicea s. l. seems located in Europe approx. 30–35 million years ago (Mao et al. 2010), suggesting a more complex scenario. This would imply a westward colonization through Mediterranean Basin with posterior fragmentation of populations during Pliocene cooling and Pleistocene glaciations. Survival populations took shelter in several isolated glacial refugia, from which they re-colonized other territories, this being clear from the different phylogeographic groups detected. Elucidate where survivals populations during last glacial maximum were located could help to understand the present distribution of genetic diversity in J. phoenicea. Paleodistribution modeling in conjunction with population genetic analyses could predict the past distributions and aid in locating Pleistocene refugia of J.grex phoenicea. However, it requires other techniques (for example, ecological niche model studies) that should be taken into account in further studies.