Skip to main content

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


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.


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.


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.


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.


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.


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.


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. turbinataJ. 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.

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).

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.

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 105 and 106 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.


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
figure 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

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
figure 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
figure 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
figure 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


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.


  • Adams RP (2014) Junipers of the world: the genus Juniperus, 4th edn. Trafford Publishing Co., Bloomington

    Google Scholar 

  • Adams RP, Schwarzbach AE (2013) Phylogeny of Juniperus using nrDNA and four cpDNA regions. Phytologia 95:179–187

    Google Scholar 

  • Adams RP, Padney N, Rezzi S, Casanova J (2002) Geographic variation in the random amplified polymorphic DNAs (RAPDs) of Juniperus phoenicea, J. p. var. canariensis, J. p. subsp. eumediterranea and J. p. var turbinata. Biochem Syst Ecol 30:223–229

    CAS  Article  Google Scholar 

  • Adams RP, Nguyen S, Achak N (2006) Geographic variation in Juniperus phoenicea (Cupressaceae) from the Canary Islands, Morocco and Spain based on RAPDs analysis. Phytologia 88:270–278

    Google Scholar 

  • Adams RP, Rumeu BR, Nogales M, Fontinha SS (2010) Geographic variation and systematics of Juniperus phoenicea L. from Madeira and the Canary Islands: analyses of SNPs from nrDNA and petN-psbM DNA. Phytologia 92:59–67

    Google Scholar 

  • Adams RP, Boratyński A, Arista M, Schwarzbach AE, Leschner H, Liber Z, Minissale P, Mataraci T, Manolis A (2013) Analysis of Juniperus phoenicea from throughout its range in the Mediterranean using DNA sequence data from nrDNA and petN-psbM: the case for the recognition of J. turbinata Guss. Phytologia 95:202–209

    Google Scholar 

  • Arista M, Ortiz PL, Talavera S (1997) Reproductive isolation of two sympatric subspecies of Juniperus phoenicea (Cupressaceae) in southern Spain. Plant Syst Evol 208:225–237

    Article  Google Scholar 

  • Blondel J, Aronson J (1999) Biology and wildlife of the Mediterranean region. Oxford University Press, New York

    Google Scholar 

  • Boratyński A, Lewandowski A, Boratyńska K, Montserrat JM, Romo A (2009) High level of genetic differentiation of Juniperus phoenicea (Cupressaceae) in the Mediterranean region: geographic implications. Plant Syst Evol 277:163–172

    Article  Google Scholar 

  • Browicz K, Zieliński J (1982) Chorology of trees and shrubs in South-West Asia and adjacent regions. Polish Scientific Publishers, Warszawa-Poznan

    Google Scholar 

  • Burban C, Petit RJ (2003) Phylogeography of maritime pine inferred with organelle markers having contrasted inheritance. Mol Ecol 12:1487–1495

    CAS  Article  PubMed  Google Scholar 

  • Charco J (2001) Guía de los árboles y arbustos del norte de África. Agencia Española de Cooperación Internacional, Madrid

    Google Scholar 

  • Christensen KI (1997) Juniperus L. In: Strid A, Tan K (eds) Flora hellenica, vol 1. Koeltz Scientific Books, Königstein, pp 10–14

    Google Scholar 

  • Conord C, Gurevitch J, Fady B (2012) Large-scale longitudinal gradients of genetic diversity: a meta-analysis across six phyla in the Mediterranean basin. Ecol Evol 2:2595–2609.

    Article  Google Scholar 

  • Corander J, Marttinen P (2006) Bayesian identification of admixture events using multilocus molecular markers. Mol Ecol 15:2833–2843

    Article  PubMed  Google Scholar 

  • do Amaral Franco J (1986) Juniperus L. In: Castroviejo S, Laínz M, López González G, Montserrat P, Muñoz Garmendia F, Paiva J, Villar L (eds) Flora iberica, vol 1. Real Jardín Botánico, CSIC, Madrid, pp 181–188

    Google Scholar 

  • Doyle JJ, Doyle JL (1987) A rapid DNA isolation procedure for small quantities of fresh leaf material. Phytochem Bull 19:11–15

    Google Scholar 

  • Dzialuk A, Mazur M, Boratyńska K, Monstserrat JM, Romo A, Boratyński A (2011) Population genetic structure of Juniperus phoenicea (Cupressaceae) in the western Mediterranean Basin: gradient of diversity on a broad geographical scale. Ann For Sci 68:1341–1350

    Article  Google Scholar 

  • Earl DA, von Holdt BM (2012) Structure Harvester: a website and program for visualizing Structure output and implementing the Evanno method. Conserv Genet Resour 4:359–361. (available at

    Article  Google Scholar 

  • Evanno G, Regnaut S, Goudet J (2005) Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol 14:2611–2620.

    CAS  Article  PubMed  Google Scholar 

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

    Article  PubMed  Google Scholar 

  • Excoffier L, Smouse PE, Quattro JM (1992) Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondiral DNA restriction data. Genetics 131:479–491

    CAS  PubMed  PubMed Central  Google Scholar 

  • Fady B, Conord C (2010) Macroecological patterns of species and genetic diversity in vascular plants of the Mediterranean Basin. Divers Distrib 16:53–64

    Article  Google Scholar 

  • Fady-Welterlen B (2005) Is there really more biodiversity in Mediterranean forest ecosystems? Taxon 54:905–910

    Article  Google Scholar 

  • Falush D, Stephens M, Pritchard JK (2003) Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies. Genetics 164:1567–1587

    CAS  PubMed  PubMed Central  Google Scholar 

  • Falush D, Stephens M, Pritchard JK (2007) Inference of population structure using multilocus genotype data: dominant markers and null alleles. Mol Ecol Notes 7:574–578.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  • Farjon A (2005) A monograph of Cupressaceae and Sciadopitys. Royal Botanic Gardens, Kew

    Google Scholar 

  • Felsenstein J (1989) PHYLIP—phylogeny inference package (version 3.2). Cladistics 5:164–166

    Google Scholar 

  • Greuter W, Burdet HM, Long G (1984) Med.-Checklist 1. Conservatoire et Jardin Botaniques de la Ville de Genève et Med.-Checklist Trust of OPTIMA, Genève

  • Hammer Ø, Harper DAT, Ryan PD (2001) PAST: paleontological statistics software package for education and data analysis. Palaeontol Electron 4:1–9

    Google Scholar 

  • Hamrick JL, Godt MJW, Sherman-Broyles SL (1992) Factors influencing levels of genetic diversity in woody plant species. New For 6:95–124

    Article  Google Scholar 

  • Jakobsson M, Rosenberg NA (2007) CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics 23:1801–1806.

    CAS  Article  PubMed  Google Scholar 

  • Jalas J, Suominen J (eds) (1973) Atlas florae europaeae 2. The Committee for Mapping the Flora of Europe and Societas Biologica Fennica Vanamo, Helsinki

  • Jiménez JF, Werner O, Sánchez-Gómez P, Fernández S, Guerra J (2003) Genetic variation and migration pathways of Juniperus thurifera L. (Cupressaceae) in the western Mediterranean region. Israel J Plant Sci 51:11–22

    Article  Google Scholar 

  • Jiménez JF, Sánchez-Gómez P, Cánovas JL, Hensen I, Aouissat M (2017) Influence of natural habitat fragmentation on the genetic structure of Canarian populations of Juniperus turbinata. Silva Fenn.

  • Jordano P (1993) Geographic ecology and variation of plant-seed disperser interaction: southern Spanish junipers and frugivorous thrushes. Vegetatio 107/108:85–104

    Google Scholar 

  • Koopman W (2005) Phylogenetic signal in AFLP data sets. Syst Biol 54:197–217.

    Article  PubMed  Google Scholar 

  • Korol L, Kara N, Isik K, Schiller G (1997) Differentiation among and within natural and planted Cupressus sempervirens L. eastern Mediterranean populations. Silvae Genet 46:151–155

    Google Scholar 

  • Krauss SL (2000) Accurate gene diversity estimated from amplified fragment length polymorphism (AFLP) markers. Mol Ecol 9:1241–1245

    CAS  Article  PubMed  Google Scholar 

  • Krutovskii KV, Erofeeva SY, Aagaard JE, Strauss SH (1999) Simulation of effects of dominance on estimates of population genetic diversity and differentiation. J Hered 90:499–502.

    Article  Google Scholar 

  • Lagercrantz U, Ryman N (1990) Genetic structure of Norway spruce (Picea abies): concordance of morphological and allozymatic variation. Evolution 44:38–53

    PubMed  Google Scholar 

  • Lebreton P, Pérez de Paz PL (2001) Définition du Genévrier de Phénicie (Juniperus aggr. phoenicea), reconsidéré à ses limites biogéographiques: Méditerranée orientale (Crète et Chypre) et Atlantique (Iles Canaries). Bull Mens Soc Linn Lyon 70:73–92.

    Article  Google Scholar 

  • Linares JC (2011) Biogeography and evolution of Abies (Pinaceae) in the Mediterranean Basin: the roles of long-term climatic change and glacial refugia. J Biogeogr 38:619–630.

    Article  Google Scholar 

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

    CAS  Article  PubMed  Google Scholar 

  • Lynch M, Milligan BG (1994) Analysis of population genetic structure with RAPD markers. Mol Ecol 3:91–99

    CAS  Article  PubMed  Google Scholar 

  • Mantel N (1967) The detection of disease clustering and a generalized regression approach. Cancer Res 27:209–220

    CAS  PubMed  Google Scholar 

  • Mao K, Hao G, Liu J, Adams RP, Milne RI (2010) Diversification and biogeography of Juniperus (Cupressaceae): variable diversification rates and multiple intercontinental dispersals. New Phytol 188:254–272.

    CAS  Article  PubMed  Google Scholar 

  • Mazur M, Klajbor K, Kielich M, Sowinska M, Romo A, Montserrat JM, Boratyński A (2010) Intra-specific differentiation of Juniperus phoenicea in the western Mediterranean region revealed in morphological multivariate analysis. Dendrobiology 63:21–31

    Google Scholar 

  • Mazur M, Minissale P, Sciandrello S, Boratyński A (2016) Morphological and ecological comparison of populations of Juniperus turbinata Guss. and J. phoenicea L. from the Mediterranean region. Plant Biosyst 150:313–322.

    Article  Google Scholar 

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

    Article  Google Scholar 

  • Médail F, Quézel P (1997) Hot-spots analysis for conservation of plant biodiversity in the Mediterranean basin. Ann Mo Bot Gard 84:112–127

    Article  Google Scholar 

  • Meloni M, Perini D, Filigheddu R, Binelli G (2006) Genetic variation in five Mediterranean populations of Juniperus phoenicea as revealed by inter-simple sequence repeat (ISSR) markers. Ann Bot 97:299–304

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  • Meudt HM, Clarke AC (2007) Almost forgotten or latest practice? AFLP applications, analyses and advances. Trends Plant Sci 12:106–117.

    CAS  Article  PubMed  Google Scholar 

  • Mueller UG, Wolfenbarger L (1999) AFLP genotyping and fingerprinting. Trends Ecol Evol 14:389–394

    CAS  Article  PubMed  Google Scholar 

  • Müller-Starck G, Baradat P, Bergmann F (1992) Genetic variation within European tree species. New For 6:23–47

    Article  Google Scholar 

  • Nei M (1978) Estimation of average heterozygosity and genetic distance from a small number of individuals. Genetics 89:583–590

    CAS  PubMed  PubMed Central  Google Scholar 

  • Nieto Feliner G (2014) Patterns and processes in plant phylogeography in the Mediterranean Basin. A review. Perspect Plant Ecol 16:265–278.

    Article  Google Scholar 

  • Nybom H (2004) Comparison of different nuclear DNA markers for estimating intraspecific genetic diversity in plants. Mol Ecol 13:1143–1155.

    CAS  Article  PubMed  Google Scholar 

  • Nybom H, Bartish I (2000) Effects of life history traits and sampling strategies on genetic diversity estimates obtained with RAPD markers in plants. Perspect Plant Ecol 3:93–114.

    Article  Google Scholar 

  • Padilla DP, González-Castro A, Nogales M (2012) Significance and extent of secondary seed dispersal by predatory birds on oceanic islands: the case of the Canary archipelago. J Ecol 100:416–427.

    Article  Google Scholar 

  • Panagiotopoulos K, Böhm A, Leng MJ, Wagner B, Schäbitz F (2014) Climate variability over the last 92 ka in SW Balkans from analysis of sediments from Lake Prespa. Clim Past 10:643–660.

    Article  Google Scholar 

  • Peakall R, Smouse PE (2006) Genalex 6: genetic analysis in Excel. Population genetic software for teaching and research. Mol Ecol Notes 6:288–295.

    Article  Google Scholar 

  • Petit RJ, Csaikl UM, Bordács S, Burg K, Coart E, Cottrell J, van Dam B, Deans JD, Dumolin-Lapegue S, Fineschi S, Finkeldey R, Gillies A, Glaz I, Goicoechea PG, Jensen JS, Konig AO, Lowe AJ, Madsen SF, Matyas G, Munro RC, Olalde M, Pemonge MH, Popescu F, Slade D, Tabbener H, Taurchini D, de Vries SGM, Ziegenhagen B, Kremer A (2002) Chloroplast DNA variation in European white oaks: phylogeography and patterns of diversity based on data from over 2600 populations. For Ecol Manag 156:5–26

    Article  Google Scholar 

  • Pritchard JK, Wen W (2004) Documentation for STRUCTURE software: version 2. University of Chicago, Chicago

    Google Scholar 

  • Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genotype data. Genetics 155:945–959

    CAS  PubMed  PubMed Central  Google Scholar 

  • Quézel P, Pesson P (1980) Biogeographie et ecologie des coniferes sur le poutrour mediterraneen. In: Pesson P (ed) Actualites d’ecologie forestiere. Gauthier-Villars, Paris, pp 205–255

    Google Scholar 

  • Quézel P, Barbero M, Benabid A, Rivas-Martines S (1992) Contributions à l’étude des groupements forestiers et pre-forestiers du Maroc oriental. Studia Bot 10:57–90

    Google Scholar 

  • Reyes-Betancort JA, Santos-Guerra A, Rosana-Guma I, Humphries CJ, Carine MA (2008) Diversity, rarity and the evolution and conservation of the Canary Islands endemic flora. Anales Jard Bot Madrid 65:25–45

    Google Scholar 

  • Rivas-Martínez S, Wildpret M, del Arco M, Rodríguez O, Pérez-de-Paz PL, García-Gallo A, Acebes JR, Díaz-González TE, Fernández-González F (1993) Las comunidades vegetales de la Isla de Tenerife (Islas Canarias). Itinera Geobotanica 7:169–374

    Google Scholar 

  • Rosenberg NA (2004) DISTRUCT: a program for the graphical display of population structure. Mol Ecol Notes 4:137–138.

    Article  Google Scholar 

  • Sánchez-Gómez P, Soriano MC, Sotomayor JA, García-Vallejo MC, Sáez F, Correal E (1993) Estudio del Aceite esencial de diversas Cupresáceas autóctonas del SE. de España. Actas del I Congreso Forestal Español 4:359–364

    Google Scholar 

  • Sánchez-Gómez P, Jiménez JF, Vera JB, Sánchez-Saorín FJ, Martínez JF, Buhagiar J (2013) Genetic structure of Tetraclinis articulata, an endangered conifer of the western Mediterranean basin. Silva Fenn.

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

    Article  Google Scholar 

  • Santos-Guerra A (2001) Flora vascular nativa. In: Fernández-Palacios JM, Martín Esquivel JL (eds) Naturaleza de las Islas Canarias. Ecología y Conservación, Santa Cruz de Tenerife, pp 185–192

    Google Scholar 

  • Sękiewicz K, Dering M, Sękiewicz M, Boratyńska K, Iszkulo G, Litkowiec M, Ok T, Dagher-Kharrat MB, Boratyński A (2015) Effect of geographic range discontinuity on species differentiation—East-Mediterranean Abies cilicica: a case study. Tree Genet Genomes 11:810.

    Article  Google Scholar 

  • Sobierajska K, Boratyńska K, Jasińska A, Dering M, Ok T, Douaihy B, Dagher-Kharrat MB, Romo A, Boratyński A (2016) Effect of the Aegean Sea barrier between Europe and Asia on differentiation in Juniperus drupacea (Cupressaceae). Bot J Linn Soc 180:365–385.

    Article  Google Scholar 

  • Stuessy TF, Jakubowsky G, Salguero-Gómez R, Pfosser M, Schluter PM, Fer T, Sun B-Y, Kato H (2006) Anagenetic evolution in island plants. J Biogeogr 33:1259–1265

    Article  Google Scholar 

  • Sytsma KJ, Givnish TJ, Smith JF, Hain WJ (1993) Collection and storage of land plant samples for macromolecular comparisons. Methods Enzymol 224:23–37

    CAS  Article  PubMed  Google Scholar 

  • Terrab A, Schönswetter P, Talavera S, Vela E, Stuessey TF (2008) Range-wide phylogeography of Juniperus thurifera L., a presumptive keystone species of western Mediterranean vegetation during cold stages of the Pleistocene. Mol Phylogenet Evol 48:94–102

    CAS  Article  PubMed  Google Scholar 

  • Thompson JD (2005) Plant evolution in the Mediterranean. Oxford University Press, Oxford

    Book  Google Scholar 

  • Van Andel TH (2002) The climate and landscape of middle part of Weichselian glaciation in Europe: the stage 3 project. Quat Res 57:2–8

    Article  Google Scholar 

  • Vekemans X (2002) AFLP-SURV version 1.0. Distributed by the author. Laboratoire de Génétique et Ecologie Végétale, Université Libre de Bruxelles, Bruxelles

    Google Scholar 

  • Vos P, Hogers R, Bleeker M, Reijans M, van de Lee T, Hornes M, Frijters A, Pot J, Peleman J, Kuiper M, Zabeau M (1995) AFLP: a new technique for DNA fingerprinting. Nucleic Acid Res 11:4407–4414

    Article  Google Scholar 

  • Zamora C, Grandal A (1997) Reconstrucción de la vegetación potencial del Campo de Cartagena a la luz de la documentación de su archivo municipal. An Biol 22:67–76

    Google Scholar 

Download references


Authors wish to thank to Jaime Güemes, Gianluigi Bachetta, Jorge Balibrea, Laura Aznar, Jesús Muñoz, Mohammad Al-Gharaibeh, and David López for field assistance in collecting samples.


This work has been supported by a Ministry of Economy and Competitiveness CGL2011-30099 project grant to Pedro Sánchez Gómez.

Author information

Authors and Affiliations


Corresponding author

Correspondence to Juan F. Jiménez.

Ethics declarations

Conflict of interest

The authors declare that they have no conflict of interest.

Additional information

Handling Editor: Ricardo Alia

Contribution of the co-authors

J. F. Jiménez has contributed as co-supervisor of the project, running the data analysis and editing the manuscript and figures. P. Sánchez-Gómez is responsible of designing the work and editing the manuscript. J. L. Cánovas has contributed running the data and editing manuscript, tables and figures. J. B. Vera, I. Hensen and M. Aouissat have contributed designing the work and running several analyses. All authors are also responsible of collecting samples in the field and writing the manuscript.

Electronic supplementary material

Fig. S1
figure 5

Map of sampled populations. Codes of populations are listed in Table 1 (PNG 1.28 mb)

Fig. S2
figure 6

Neighbour-joining phenogram based on Nei's unbiased genetic distance for the 46 populations. Bootstrap values (>70%) over loci (based on 1000 replicates) are indicated for each node (PNG 1.39 mb)

Table S1

Pairwise genetic differentiation among populations. Upper section, Fst values obtained with AMOVA. Lower section, Fst values obtained with AFLPsurv. (XLS 57 kb)

High resolution image (TIF 690 kb)

High resolution image (TIF 536 kb)

Rights and permissions

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Sánchez-Gómez, P., Jiménez, J.F., Cánovas, J.L. et al. Genetic structure and phylogeography of Juniperus phoenicea complex throughout Mediterranean and Macaronesian regions: different stories in one. Annals of Forest Science 75, 75 (2018).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • AFLP markers
  • Biogeography
  • Cupressaceae
  • Genetic diversity
  • Inter-populational differentiation
  • Taxonomy