High genetic diversity with moderate differentiation in Juniperus excelsa from Lebanon and the eastern Mediterranean region.

BACKGROUND AND AIMS: Juniperus excelsa is an important woody species in the high mountain ecosystems of the eastern Mediterranean Basin where it constitutes the only coniferous species found at the tree line. The genetic diversity within and among J. excelsa populations of the eastern Mediterranean Basin is studied in the light of their historical fragmentation.
METHODOLOGY: Nuclear microsatellites originally developed for Juniperus communis and J. przewalskii were tested on 320 individuals from 12 different populations originating from Lebanon, Turkey, Cyprus, Greece and the Ukraine. PRINCIPAL
RESULTS: Among the 31 nuclear microsatellite primers tested, only three produced specific amplification products, with orthology confirmed by sequence analysis. They were then used for genetic diversity studies. The mean number of alleles and the expected heterozygosity means were N(a)=8.78 and H(e)=0.76, respectively. The fixation index showed a significant deviation from Hardy-Weinberg equilibrium and an excess of homozygotes (F(IS)=0.27-0.56). A moderate level of genetic differentiation was observed among the populations (F(ST)=0.075, P<0.001). The most differentiated populations corresponded to old vestigial stands found at the tree line (>2000 m) in Lebanon. These populations were differentiated from the other populations that are grouped into three sub-clusters.
CONCLUSIONS: High levels of genetic diversity were observed at species and population levels. The high level of differentiation in the high-mountain Lebanese populations reflects a long period of isolation or possibly a different origin. The admixture observed in other populations from Lebanon suggests a more recent separation from the Turkish-southeastern European populations.

Introduction

The Mediterranean Basin is one of 34 world biodiversity hotspots (Myers ; Myers 2003). Its floristic richness is primarily due to unique climatic and habitat heterogeneity, historical factors and the different origins of the flora itself (Quézel 1995). It is considered a global refuge with many Pleistocene refugia, which are key areas for the conservation of plant biodiversity (Médail and Diadema 2009).

The evaluation of genetic diversity and its distribution within and among populations is essential for providing information to implement strategies in breeding and genetic resource conservation programmes of plant taxa (Petit ; Bruschi ; Meloni ). The genetic studies on woody plant populations in the eastern Mediterranean Basin show a high level of within-population diversity (Fady and Conkle 1993; Bou Dagher-Kharrat ; Fady ). It has been suggested that these woody formations may have experienced weak population bottlenecks (Fady and Conord 2010) due to mild climatic oscillations during the Quaternary (Sanchez-Goni ; Van Andel 2002), and that this may explain the higher genetic diversity observed in the eastern basin as compared with that observed in woody species in the western part of the basin. On the other hand, the woodlands on the eastern side of the Mediterranean have historically suffered from strong human impacts that have led to significant forest fragmentation, as in the Juniperus excelsa woodlands, a major element of the mountainous conifer forests in the eastern Mediterranean Basin (Barbero ).

The old juniper woodlands, especially at the higher elevations, are of great biogeographical interest, being the remnants of more widespread ancestral pre-glacial juniper woodlands (Jalut ; Quézel and Médail 2003; Eastwood 2004). Moreover, they also have a high ecological value, being frequently the only tree species able to grow in semi-arid environments and therefore playing an important role in soil protection.

Juniperus excelsa sensu lato has a wide distribution (Athanasiadis 1986; Boratyński ; Christensen 1997). Two subspecies are currently recognized (Farjon 2005, 2010). The first is a western taxon, J. excelsa subsp. excelsa, which extends from the eastern Mediterranean Basin to Crimea in Ukraine, and more to the east reaching the Caucasus mountains, the Elburtz mountains in Iran and the Kopet mountains in Turkmenistan (Fig. 1). The second subspecies is an eastern taxon, J. excelsa subsp. polycarpos (Farjon 2005) or J. polycarpos K. Koch (Adams 2008), which has a Trans-Caucasian–Central Asian distribution.

Fig. 1

Distribution range of J. excelsa subsp. excelsa (after Browicz 1982; Boratyński ; supplemented with data of Farjon 2005).

In this study, we focused on the western taxon J. excelsa subsp. excelsa.

Juniperus excelsa is found in rare and fragmented woodlands in the southern and central Balkans and the Cyprus mountains (Milios ). The species is regionally widespread and continuously distributed along the Taurus chain in southern Anatolia while severely fragmented on the Anatolian plateau and along the Syrian and Lebanese mountains (Quézel and Médail 2003). It is present in two of the Mediterranean region's biodiversity hotspots, the south Anatolia and Cyprus hotspot, and the Syria–Lebanon–Israel/Palestine hotspot. This species grows in the mountains and is considered as putative glacial refugia in the Taurus–Ammanus chain and Cyprus–Lebanon mountains (Médail and Diadema 2009).

Juniperus excelsa is an arborescent juniper that can reach 20–25 m in height. It is slow growing, monoecious or dioecious, and wind pollinated (Farjon 2005; Adams 2008), with seeds being dispersed by gravity or at longer distances by birds and small mammals (Jordano 1993; Santos ). The species shows a wide range of climatic plasticity, colonizing sites that vary from sub-humid to the adjacent semi-arid steppe zone of the Mediterranean region (Quézel 1973; Abi-Saleh ; Akman ; Quézel and Médail 2003). Juniperus excelsa is most frequently found in cold Mediterranean zones (Barbero ) in lower sub-humid climates, at elevations between 1000 and 1300 m in the Anatolian forests, and between 1600 and 1800 m on the eastern slope of Mount Lebanon (Fig. 2A). It is capable of tolerating severe drought and cold conditions and can grow on shallow, degraded soils. In the oro-Mediterranean zone, J. excelsa is the dominant tree species found at the tree line in the eastern Mediterranean region with very sparse pure vestigial stands. It can reach elevations of 2100 m in Greece and some individuals can be found at elevations of 2700–2800 m in the Taurus and Mount Lebanon (Quézel 1973; Abi-Saleh , 1996; Akman ; Browicz 1982; Barbero ) (Fig. 2B).

Fig. 2

(A) Dense formation at 1600 m altitude and (B) old, sparse formation on the tree line at 2300 m altitude.

No previous studies have been performed to characterize the genetic variability of J. excelsa populations across the eastern Mediterranean Basin. Despite being the world's most widespread conifer genus, genetic studies on Juniperus in general remain scarce. The most extensive genetic investigations on Juniperus were focused on the phylogeny and phylogeography of the genus, studying inter-specific (Adams 2008; Mao ) and intra-specific (Opgenoorth ) differentiation based on cpDNA and internal transcribed spacer (ITS) markers. Only five juniper species have been studied for their within-species genetic diversity: Juniperus thurifera using amplified fragment length polymorphism (AFLP) (Jiménez ; Terrab ), J. phoenicia using inter-simple sequence repeat (ISSR) and isozymes (Meloni ; Boratyński ), J. rigida and J. coreana using enzyme electrophoresis (Huh and Huh 2000), and J. communis using AFLP and nuclear simple sequence repeat (nSSR) (Van Der Merwe ; Michalczyk 2008). These studies showed high levels of genetic diversity in the Juniperus species.

Since microsatellites are considered the markers of choice for studying intra-specific genetic diversity and genetic structure (Varshneya ; Morgante and Olivieri 2008), we used these markers to (i) investigate the transferability of microsatellite markers among Juniperus species, (ii) evaluate patterns of genetic variation within and among different J. excelsa populations and (iii) study the phylogeography of J. excelsa in the eastern Mediterranean. This information is essential to provide information on the historical processes shaping genetic diversity and to plan successful conservation strategies and reforestation programmes.

Materials and methods

Plant material

Twelve populations of J. excelsa were sampled encompassing almost the entire natural range of the species (Fig. 3, Table 1). The most eastern populations (e.g. Eastern Turkey, Armenia, Iran, etc.) were not included. The plant material collected in Lebanon was stored at −80 °C while samples from the other countries were air dried. Genomic DNA was extracted from the leaf material of 15–30 trees per population using a modified cetyltrimethyl ammonium bromide protocol (Bou Dagher-Kharrat ).

Fig. 3

Proportional assignment of individuals to the three genetic clusters as detected in an admixture analysis using a Bayesian model in the program Structure™. Inside the dashed box: second-order rate of change of the log likelihood of the data (ΔK) as a function of K, the number of clusters.

Table 1

Geographic locations of the 12 analysed J. excelsa populations

Population label	Country	Locality	Longitude	Latitude	Elevation (m)
LB1	Lebanon	Qammouaa	N34°29′34′	E36°15′14′	1450–1800
LB2		Danniyeh	N34°23′17′	E36°05′60′	1600–1850
LB3		Wadi El Njassa	N34°19′49′	E36°03′16′	1870–2300
LB4		Barqa	N34°11′48.4′	E36°8′15′	1600–2200
LB5		Afqa	N34°4′25′	E35°54′20′	1100–1600
LB6		Aarsala	N34°4′57′	E36°28′34′	2180
TU1	Turkey	IlgazTosya	N40°53′04′	E33°42′24′	850
TU2	Turkey	Eğirdir	N38°08′12	E30°46′42′	950
GR	Greece	Askion Oros	N40°15′58′	E21°37′26′	1000
CY	Cyprus	Troodos Oros	N34°55′20′	E33°05′55′	1500
CR1	Ukraine	Crimea Mys Aja	N44°25′18′	E33°39′57′	30–40
CR2	Ukraine	Crimea Kolkhoznoe	N44°29′00′	E33°49′54′	500

aHigh-elevation populations.

Microsatellites

Thirty-one microsatellite primer pairs isolated and characterized in other Juniperus species were tested: 27 microsatellites from J. communis (Bérubé ; Michalczyk ; I. M. Michalczyk, Philipps-University of Marburg, Germany, pers. comm.) and four from J. przewalskii (Zhang [Additional information]. Polymerase chain reaction (PCR) amplification was performed in a total volume of 12.5 μL, containing 25 ng of DNA template, 1 × PCR buffer (5× Green Go Taq® Reaction Buffer), 500 μM each deoxynucleotide triphosphate, 1 U Taq polymerase (Go Taq DNA polymerase (5 U μL-1) Promega), 1.5 mM MgCl2 and 0.5 μM each primer. Polymerase chain reaction amplification was performed using a Gene Amp PCR System 9700 set with an initial denaturation step at 94 °C for 5 min followed by 40 cycles at 94 °C for 30 s, 55 °C for 50 s, 72 °C for 50 s and a final extension step at 72 °C for 7 min. The amplification products were separated by capillary electrophoresis using a Megabace1000 (GE Healthcare) automated sequencer. Alleles were sized using FRAGMENT PROFILER version 1.2 (GE Healthcare).

The genotyping was done at the Institute of Plant Genetics, Consiglio Nazionale delle Ricerche of Florence, Italy.

Genetic analyses

GenAlEx™ software (Peakall and Smouse 2005) was used to estimate the following genetic diversity parameters: mean number of alleles (Na), effective number of alleles (Ne=1/(1− He)) (Brown and Weir 1983), expected (He) and observed (Ho) heterozygosity according to Nei (1987), and total genetic diversity Ht=(1−∑p2) (p2=mean allelic frequency of the allele (i)). Null allele frequencies were estimated using FreeNA (Chapuis and Estoup 2007). We used the excluding null alleles (ENA) correction (Chapuis and Estoup 2007) method to estimate FST in the absence of null alleles. The allelic richness per locus and population (Ar) was also estimated based on the minimum sample size of the data set using the rarefaction method of El Mousadik and Petit (1996) and the software Fstat (http://www.unil.ch/izea/softwares/fstat.html). The inbreeding coefficient (FIS=1−(Ho/He) and the deviation from Hardy–Weinberg expectations were determined according to Hedrick (2000).

The mean values of the different genetic diversity parameters were compared among the populations using the Kruskal–Wallis test. An analysis of molecular variance was performed to evaluate the total amount of variance explained by differences among populations. The genetic similarities between populations were inferred following a model-based Bayesian assignment using the software Structure™ (Pritchard ). To verify the consistency of results among different runs for a given K value, we calculated the similarity coefficient (SC) between run pairs as described by Rosenberg . Furthermore, to identify the number of clusters (K) that best explain the data, we calculated the rate of change of L(K) (ΔK) between successive K values, as shown by Evanno . Computation of ΔK and SC values was carried out according to Ehrich (2006) using the R package structure-sum (www.nhm.uio.no/ncb). A principal coordinate analysis (PCA) was performed on the pairwise FST distances following the algorithm of Orloci (1978). A dendrogram with 100 bootstraps was constructed based on Nei's genetic distances. Furthermore, to estimate the correlation between geographic and genetic distance, a Mantel test (Smouse and Long 1992) based on 9999 permutations was performed on the matrix of FST and Nei genetic distance values and the geographic distance matrix.

Results

Cross-species amplifications using microsatellite primers had limited success. When transferred from J. przewalskii to J. excelsa, one out of four primer pairs tested revealed amplification and the amplification product was monomorphic. On the other hand, only six out of 27 microsatellite primers succeeded in generating amplification products when transferred from J. communis to J. excelsa [Additional information].

Only three primer pairs (Jc031, Jc037 and Jc166) were polymorphic in J. excelsa (Table 2). The sequencing of these three transferred microsatellites confirmed the presence of the repeat stretches in J. excelsa (Table 2).

Table 2

The three nuclear microsatellite primers used for J. excelsa genotyping

Locus name	Primer sequence (5′–3′)	Motifa	Allele size range (bp)
			J. communis	J. excelsa
Jc031	F: FamCCTAATGTTGTAATCACGTATATCT R:TGACCTTGGGCGTATAGATT	(CA)15	174–242	150–196
Jc037	F: FamGGCAATTAGTAAGGCACAAG R:TAAGGTGGATATCACCAAGG	(TG)9(AG)22	176–222	153–197
Jc166	F: HexATTTGTTTTCTTGTGGATGC R: GCACTGACACCTATATGCAC	(TG)14	170	153–191

aRepeat numbers as observed in the J. excelsa sequences are reported.

The number of alleles (Na) observed for J. communis at the locus Jc031 ranged between 10 and 25, and at the locus Jc037 between 25 and 35 (Michalczyk 2008). These values are significantly higher than those obtained in our study for J. excelsa with Na between 4 and 12 at the locus Jc031, and between 11 and 16 at the locus Jc037. A higher number of alleles are usually observed for the focal species, for which the primers were developed (Jarne and Lagoda 1996).

Moreover, transferred microsatellites are expected to show a substantial level of null alleles (Oddou-Muratorio ). The estimated null allele frequency was intermediate to high for the locus Jc031, ranging from 20 (LB6) to 36 % (GR), and low to intermediate for Jc166 and Jc037, ranging from 0 (LB5) to 21 % (GR) and from 5 (CR2) to 26 % (LB1), respectively.

Genetic diversity

Considering the 12 populations included in this study, the mean number of alleles over the three microsatellites (Na) was 8.78 and ranged between 6.33 (CR1) and 11 (LB6). The mean effective number of alleles (Ne) ranged between 4.17 (CR2) and 6.48 (LB1). The rarefied allelic richness (Ar) ranged between 5.20 (CR2) and 7.82 (LB6). The mean value of observed heterozygosity (Ho) was 0.46 and expected heterozygosity (He) and total genetic diversity (Ht) were both high with mean values of 0.76 and 0.84, respectively (Table 3). The Kruskal–Wallis test showed no significant differences between the populations for all the genetic diversity parameters. A significant deviation from the Hardy–Weinberg equilibrium was observed at the population level, with an excess of homozygotes. Considering the mean value over the three loci, the inbreeding coefficient FIS ranged from 0.27 (LB6) to 0.56 (GR) (Table 3). After removing the potential null heterozygotes from the data set, FIS showed lower values varying from −0.092 for LB3 to 0.188 for LB5. Although FIS still tends to be positive, the values obtained were, in most cases, not significant.

Table 3

Genetic diversity, inbreeding coefficient and null allele frequency estimates

Population label	Locus	N	Na	Ne	Ar	Ho	He	FIS	FIS without null alleles	NA
LB1	Jc166	30	6	2.86	4.47	0.53	0.65	0.18 ns	−0.013 ns	0.08
	Jc31	25	9	4.61	6.70	0.28	0.78	0.64***	−0.077 ns	0.29
	Jc37	26	16	11.97	10.92	0.42	0.92	0.54***	−0.115 ns	0.26
	Mean ± SD		10.3 ± 2.9	6.48 ± 2.79	7.36 ± 1.89	0.41 ± 0.07	0.78 ± 0.077	0.45	−0.068
LB2	Jc166	27	5	2.68	4.17	0.41	0.63	0.35 ns	0.199 ns	0.14
	Jc31	21	7	4.22	5.99	0.29	0.76	0.63***	−0.032 ns	0.27
	Jc37	30	12	8.61	8.65	0.77	0.88	0.13 ns	0.070 ns	0.05
	Mean ± SD		8 ± 2.0	5.17 ± 1.78	6.27 ± 1.30	0.49 ± 0.14	0.76 ± 0.074	0.37	0.079
LB3	Jc166	29	11	3.06	6.03	0.59	0.67	0.13***	0.059**	0.01
	Jc31	26	8	4.68	6.64	0.19	0.79	0.76***	−0.186 ns	0.34
	Jc37	23	11	7.96	8.51	0.52	0.87	0.4**	−0.147 ns	0.19
	Mean ± SD		10 ± 1.0	5.23 ± 1.44	7.06 ± 0.75	0.43 ± 0.12	0.78 ± 0.058	0.43	−0.092
LB4	Jc166	29	7	3.41	5.13	0.52	0.71	0.27**	0.086**	0.09
	Jc31	23	9	4.12	6.20	0.3	0.76	0.6***	0.202 ns	0.25
	Jc37	23	16	11.38	10.96	0.65	0.91	0.29***	0.021*	0.14
	Mean ± SD		10.7 ± 2.7	6.3 ± 1.55	7.43 ± 1.79	0.49 ± 0.10	0.79 ± 0.062	0.38	0.103
LB5	Jc166	27	4	2.59	3.78	0.59	0.62	0.04 ns	0.036 ns	0
	Jc31	29	7	3.51	4.88	0.17	0.72	0.76***	0.472*	0.32
	Jc37	29	13	7.51	8.44	0.66	0.87	0.24*	0.056 ns	0.11
	Mean ± SD		8 ± 2.6	4.54 ± 1.51	5.70 ± 1.41	0.47 ± 0.15	0.73 ± 0.073	0.35	0.188
LB6	Jc166	30	10	3.4	6.23	0.67	0.71	0.06*	−0.170	0.07
	Jc31	25	12	8.93	9.15	0.52	0.89	0.41***	0.015	0.2
	Jc37	22	11	5.8	8.08	0.55	0.83	0.34***	−0.102**	0.17
	Mean ± SD		11 ± 0.6	6.04 ± 1.6	7.82 ± 0.85	0.58 ± 0.05	0.81 ± 0.053	0.27	−0.086
TU1	Jc166	26	4	3.17	3.58	0.65	0.68	0.04 ns	0.044 ns	0.01
	Jc31	25	7	3.37	5.48	0.24	0.7	0.66***	0.219 ns	0.27
	Jc37	28	13	6.27	8.58	0.5	0.84	0.41**	0.036 ns	0.18
	Mean ± SD		8 ± 2.6	4.27 ± 1.00	5.88 ± 1.46	0.47 ± 0.15	0.74 ± 0.049	0.37	0.100
TU2	Jc166	13	4	3.41	3.91	0.39	0.71	0.46*	0.101 ns	0.19
	Jc31	9	6	4.38	6.00	0.11	0.77	0.86***	0.200 ns	0.37
	Jc37	15	12	8.33	9.56	0.6	0.88	0.32*	0.057 ns	0.15
	Mean ± SD		7.3 ± 2.4	5.38 ± 1.50	6.49 ± 1.65	0.37 ± 0.14	0.79 ± 0.050	0.54	0.119
GR	Jc166	22	5	3.03	4.06	0.32	0.67	0.53***	0.184 ns	0.21
	Jc31	17	7	3.78	5.62	0.12	0.74	0.84***	0.310 ns	0.36
	Jc37	25	12	8.68	8.88	0.6	0.89	0.32*	0.053 ns	0.15
	Mean ± SD		8 ± 2.1	5.16 ± 1.77	6.19 ± 1.42	0.35 ± 0.14	0.76 ± 0.063	0.56	0.182
CY	Jc166	22	4	2.99	3.65	0.46	0.67	0.32***	−0.056 ns	0.14
	Jc31	21	9	5.04	6.68	0.43	0.8	0.47***	0.130**	0.21
	Jc37	22	15	9.13	10.13	0.73	0.89	0.18 ns	−0.018 ns	0.09
	Mean ± SD		9.3 ± 3.2	5.72 ± 1.81	6.82 ± 1.87	0.54 ± 0.10	0.79 ± 0.065	0.32	0.019
CR1	Jc166	21	4	2.65	3.63	0.48	0.62	0.24***	−0.049 ns	0.11
	Jc31	20	7	2.79	5.22	0.1	0.64	0.84***	0.378 ns	0.34
	Jc37	22	14	8.13	9.56	0.64	0.88	0.27 ns	−0.085 ns	0.13
	Mean ± SD		8.3 ± 3.0	4.52 ± 1.81	6.14 ± 1.77	0.4 ± 0.16	0.71 ± 0.082	0.45	0.081
CR2	Jc166	18	3	2.76	3.00	0.56	0.64	0.13 ns	0.021 ns	0.06
	Jc31	15	4	3.46	3.98	0.13	0.71	0.81***	0.467 ns	0.34
	Jc37	19	12	6.28	8.61	0.74	0.84	0.12 ns	0.010 ns	0.05
	Mean ± SD		6.3 ± 2.8	4.17 ± 1.08	5.20 ± 1.73	0.48 ± 0.18	0.73 ± 0.059	0.36	0.166

N, sample size; Na, number of alleles; Ne, effective number of alleles; Ar, allelic richness after rarefaction to the smallest population size; Ho, observed heterozygosity; He, expected heterozygosity; FIS, inbreeding coefficient; ns=not significant; NA, null allele frequency; SD, standard deviation.

*P<0.05, **P<0.01, ***P<0.001.

Genetic differentiation

The overall genetic differentiation at the nuclear microsatellite loci was relatively moderate (FST=0.069), and highly significant (P<0.001). After ENA correction, the overall FST was only slightly lower (FST ENA=0.063). This result indicates a non-significant effect of null alleles on genetic differentiation estimates. The highest level of divergence was observed between LB3 and GR (FST=0.129), and the lowest FST value was scored between LB1 and LB2, LB1 and LB4, and between TU2 and CY (FST=0) (Table 4). Despite the low level of genetic differentiation, a geographic pattern was observed with congruent results obtained by Structure™, PCA and the dendrogram (Figs 3–5).

Table 4

Pairwise FST between the 12 J. excelsa populations analysed using nSSR. Population acronyms as in Table 1

Population label	LB1	LB2	LB3	LB4	LB5	LB6	TU1	TU2	GR	CY	CR1	CR2
LB1	0.000
LB2	0.012	0.000
LB3	0.076	0.093	0.000
LB4	0.012	0.020	0.058	0.000
LB5	0.020	0.028	0.073	0.013	0.000
LB6	0.076	0.086	0.023	0.060	0.082	0.000
TU1	0.055	0.070	0.085	0.064	0.073	0.091	0.000
TU2	0.030	0.047	0.070	0.036	0.049	0.071	0.033	0.000
GR	0.038	0.038	0.105	0.048	0.065	0.095	0.039	0.027	0.000
CY	0.035	0.040	0.083	0.035	0.044	0.088	0.035	0.016	0.016	0.000
CR1	0.042	0.053	0.110	0.057	0.070	0.093	0.045	0.043	0.036	0.051	0.000
CR2	0.036	0.047	0.092	0.044	0.050	0.090	0.040	0.027	0.033	0.032	0.026	0.000

Fig. 4

Principal coordinates analysis via a covariance matrix with data standardization. Projection of the barycentre on the first and second axis. Three clusters can be observed: the first cluster groups the populations of J. excelsa from Turkey, Greece, Ukraine and Cyprus; the second cluster includes four populations from Lebanon; and the third cluster groups the two remaining populations from Lebanon (sample acronyms as in Table 1).

Fig. 5

Dendrogram based on Nei standard genetic distance (1972) with bootstrap support values (population acronyms as in Table 1).

The best subdivision obtained with the Bayesian approach using Structure™ corresponded to three clusters (K=3) with some admixture observed in the different populations (Fig. 3). Two populations from Lebanon (LB3 and LB6) grouped together with a high probability (70–80 %) to belong to cluster 1. The two populations from Turkey (TU1 and TU2) had a probability >50 % of belonging to cluster 3. The other populations displayed a stronger admixture with similar assignment probabilities for two different clusters. For the four populations from Lebanon (LB1, LB2, LB4 and LB5), cluster 2 was the most prevalent with probabilities of 40–60 %. In the populations from Greece, Cyprus and Ukraine, 45–48 % of samples were assigned to either cluster 2 (CR1 and CR2) or cluster 3 (GR and CY) and cluster 1 was less represented than in the Lebanese populations.

The PCA (Fig. 4) showed a clustering of the populations into three distinct groups and the first two coordinates explained 38.81 and 27.51 % of the total variation, respectively. The first group is composed of populations from Turkey, Greece, Cyprus and Ukraine; the second includes the four Lebanese populations LB1, LB2, LB4 and LB5; and the third group includes the remaining two Lebanese populations LB3 and LB6.

The dendrogram (Fig. 5) showed a significant separation of the high-altitude Lebanese group LB3 and LB6 (bootstrap value 67 %) from all the other populations that are grouped into two sub-clusters. The first sub-cluster separated the northern Turkish population (TU1) with moderate support (bootstrap value 54 %). The second sub-cluster grouped the four Lebanese populations (LB1, LB2, LB4 and LB5) together with the southern Turkish population (TU2), and the populations from Cyprus, Greece and Crimea (GR, CY, CR1 and CR2). The northern Lebanese populations (LB1 and LB2) were significantly (bootstrap value 54 %) separated from the more southern ones (LB3 and LB4).

The Mantel test showed no significant correlation between geographic and genetic distances between the populations, using either FST (R2=0.07, P=0.1) or Nei genetic distances (R2=0.02, P=0.15).

Discussion

To assess genetic variability and differentiation among different populations of J. excelsa from the eastern Mediterranean Basin, we undertook nuclear microsatellite analysis. High levels of genetic diversity were observed at species and population levels. While the differentiation level among the analysed population was low, high-elevation Mount Lebanon populations were more strongly differentiated than populations from the rest of the range.

The transferability of nSSR among conifers is difficult (e.g. González-Martínez ). The limited transferability success of nSSR from J. communis and J. przewalskii to J. excelsa is not surprising when considering the phylogenetic relationships between the three species. In fact, J. communis and J. excelsa belong to two different sections, Juniperus and Sabina, respectively, while J. przewalskii and J. excelsa belong to two different subsections I and IV (Mao ). According to the classification of Mao based on nuclear ITS (nrITS) and combined nrITS/cpDNA data, the divergence time between the sections Sabina and Juniperus is estimated at 55 and 47 million years ago between subsections I and IV. On the other hand, the three selected microsatellite primers originally developed for J. communis were demonstrated to be useful in genetic diversity studies in J. excelsa.

During the Holocene period and until recent times, a rise in temperature coupled with human impacts resulted in the fragmentation of the ancestral juniper woodlands (Quézel and Médail 2003; Eastwood 2004). This species suffers nowadays from the negative impact of human agro-pastoral activities such as intensive grazing, cutting and pruning, and land transformation. The exploitation of these forests began >1400 years ago with the establishment of agro-pastoral communities and continues today, with wood extraction for daily energy needs, house building and tool manufacture (Salamé 1957; Talhouk ; Romo and Boratyński 2005). On the eastern slopes of Mount Lebanon, the intensive wood extraction caused a drastic fragmentation and a reduction of 75 % of the juniper forests between 1965 and 1998 (Jomaa et al. 2007). Forest fragmentation is considered a major factor, inducing the loss of genetic diversity (Ellstrand and Elam 1993; Young ). Despite the fact that all J. excelsa populations sampled in this study are fragmented (except for TU2, the south Turkish population), a high degree of genetic diversity is still present.

A comparatively high degree of genetic diversity has also been observed in the Cedrus genus across its eastern Mediterranean Basin distribution. The average gene diversity (H) based on AFLP markers ranged between 0.175 and 0.350 (Bou Dagher-Kharrat ). The preservation of this high level of genetic diversity may be related to the fact that the eastern Mediterranean species and populations did not suffer from severe demographic and genetic bottlenecks (Fady-Welterlen 2005). The high within-population diversity might also be explained by a recent origin of the fragmentation coupled with a high initial level of diversity, as was suggested for Cedrus spp. (Bou Dagher-Kharrat ). In addition, some life-history traits of juniper species such as outcrossing, wind pollination, possible long-distance seed dispersal, long life span and a wide distribution may have promoted and maintained high genetic diversity (Hamrick ; Austerlitz ).

Considering these life-history traits, the high inbreeding coefficient estimated for this species was unexpected. This result could be explained by the underestimation of the number of heterozygotes because of the presence of null alleles (Jarne and Lagoda 1996; Oddou-Muratorio ). Nevertheless, the occurrence of selfing/inbreeding in J. excelsa is possible. Field observations in Lebanon have shown that J. excelsa is almost always monoic with male and female cones present on the same tree with no spatial separation within the crown. Interestingly, a high frequency of empty seed, estimated by radiography (observation of B. Douaihy, unpubl. res.), was observed in the Lebanese populations, especially in the high-altitude populations (e.g. LB3). Inbreeding depression is common in conifer species and often results in the early abortion of the embryo followed by the formation of empty seeds (Kormutak and Lindgren 1996; Williams and Savolainen 1996). It should also be stressed that J. excelsa in Lebanon occurs at low-density populations, with 165 trees ha−1 as a maximum in LB2 and 45 trees ha−1 as a minimum in LB6 (authors' field observations). The individual selfing rates for Mediterranean conifer species were found to be negatively correlated to the stand density (Restoux ).

More field and genetic investigations are needed to better understand the mating system of J. excelsa.

The moderate level of genetic differentiation among populations estimated for J. excelsa (FST=7.5 %, P<0.001) is in accordance with the values reported for other conifer species, the mean value being estimated up to 10 % for conifer species in general (Petit ) and up to 9.5 % for Mediterranean conifer species with a wide distribution in particular (Fady-Welterlen 2005). The pattern of geographical differentiation shows a significant separation of the two Lebanese high-altitude (1800–2300 m) populations (LB3 and LB6). The pattern also shows the sub-clustering of the lower-altitude Lebanese populations that are distinct from the populations from Turkey, Cyprus, Greece and Ukraine. The high degree of differentiation between the two Lebanese populations (LB3 and LB6) and all the other populations can be explained by either a long period of isolation or their possible different origins, i.e. from different refugial areas. The admixture observed in the other populations from Lebanon (LB1, LB2, LB4 and LB5) could reflect a more recent separation from the populations of Turkey, Cyprus, Crimea and Greece. A similar result was obtained for Cedrus libani (Bou Dagher-Kharrat ; Fady ) with a significant separation between Lebanese and Turkish populations, which additionally supports the hypothesis of relatively recent, separate refugia in Lebanon. The Syrian Desert was probably a barrier that could not be crossed by J. excelsa. The separation of the northern Turkish population from the southern population and the southeastern European populations is difficult to explain: it could be due to a conserved ancestral character of the northern Turkish population. A precise palaeo-botanical explanation of the observed pattern is difficult to formulate since Juniperus pollen cannot be identified at the species level (Elenga ; Carrión 2002; Tzedakis ). It is suggested that the tree-like Juniperus species of the section Sabina originated from an ancestral taxon, widespread in Europe during the Tertiary (Kvaček 2002), and that the migration to the Mediterranean region in the late Tertiary and early Quaternary led to the separation of the two vicariant species, J. excelsa in the East and J. thurifera in the West Mediterranean (Barbero ; Jiménez ; Marcysiak ). Given the ecological characteristics of J. excelsa, its geographical distribution could have been fragmented into several isolated populations, each occupying different mountain chain massifs. The areas occupied by these populations were probably smaller during warm periods and broader during cold periods of the Pleistocene, similar to what has been described for forest trees in Europe (Hewitt 1996, 2004). During the dry periods of the Holocene, J. excelsa could have replaced other more humidity-demanding tree species (Elenga ; Mudie ; Cordova ; Jalut et al. 2009) and consequently its distribution could have increased. The positive reaction of Juniperus to arid periods during the late Pleistocene/Holocene has been reported in the Iberian Peninsula (Uzquiano and Arnanz 1997; Carrión , 2004). Warm, arid environments have been described across the Mediterranean coast, the Balkans and Anatolia during the Last Glacial Maximum (Carrión 2002; Van Andel 2002; Eastwood 2004) and were also likely during the entire Pleistocene (Fady-Welterlen 2005).

Conclusions and forward look

The current conservation status of J. excelsa is ‘lower risk/least concern’ (International Union for Conservation of Nature 2010). Despite the high level of genetic diversity in J. excelsa at both the species and the population levels, the protection of J. excelsa woodlands is essential for their long-term persistence. In addition, the negative impacts of land exploitation on the persistence of this species are increased by its slow growth and generally low germination rates of juniper seeds. In the future, continued losses of the oldest individuals and a lack of regeneration could lead to a severe, genetically deleterious effect of fragmentation on local population diversity.

The preservation of J. excelsa in the different geographic regions will also ensure the preservation of a rich and diversified gene pool.

A significant level of genetic differentiation between the Lebanese and the Turkish–southeastern European populations of J. excelsa was observed. This result highlights the biogeographical importance of J. excelsa stands in the different mountain refuges of the eastern Mediterranean. Future studies should include samples from Syrian mountains to confirm the presence of significant genetic differentiation between Lebanese and Turkish populations. Adding populations from the eastern taxon, J. excelsa subsp. polycarpos (Farjon 2005) would also allow a better genetic description of J. excelsa sensu lato.

Finally, the strong differentiation observed in remnant populations at high altitudes should be investigated to better understand the factors shaping this pattern. Individuals from similar ecological habitats in geographic regions other than the Lebanese mountains, e.g. the Taurus, should be genotyped to confirm the altitudinal pattern observed in this study. Given that these populations are growing at the altitudinal limit of their distribution and with a near absence of natural regeneration, their long-term persistence is at risk. Due to the extremely limited cone production (authors’ observations) with few filled healthy seeds, seeds for reforestation can only be harvested from populations at lower altitudes. But the lower-altitude populations experience different ecological conditions, which increases the risks associated with the use of non-adapted reproductive material. To avoid the loss of the genetic pools found in the oldest populations, ex situ conservation could be considered since a high level of genetic diversity is still preserved in these populations.

Additional information

Additional information is available in the online version of this article.

Sources of funding

This work was supported by the Research Council of Saint-Joseph University of Beirut. Material collection from populations outside of Lebanon was supported by Institute of Dendrology of Kórnik, Poland.

Contribution by the authors

All the authors contributed equally.

Conflict of interest statement

None declared.