Morphological versus molecular markers to describe variability in Juniperus excelsa subsp. excelsa (Cupressaceae).

BACKGROUND AND AIMS: Juniperus excelsa M.-Bieb. is a major forest element in the mountains of the eastern part of Mediterranean and sub-Mediterranean regions. This study comprises the first morphological investigation covering a large part of the geographical range of J. excelsa and aims to verify the congruency between the morphological results and molecular results of a previous study.
METHODOLOGY: We studied 14 populations sampled from Greece, Cyprus, Ukraine, Turkey and Lebanon, 11 of which have previously been investigated using molecular markers. Three hundred and ninety-four individuals of J. excelsa were examined using nine biometric features characterizing cones, seeds and shoots, and eight derived ratios. Statistical analyses were conducted in order to evaluate the intra- and inter-population morphological variability. PRINCIPAL
RESULTS: The level of intra-population variability observed did not show any geographical trends. The total variation mostly depended on the ratios of cone diameter/seed width and seed width/seed length. The discrimination analysis, the Ward agglomeration method and barrier analysis results showed a separation of the sampled populations into three main clusters. These results confirmed, in part, the geographical differentiation revealed by molecular markers with a lower level of differentiation and a less clear geographical pattern. The most differentiated populations using both markers corresponded to old, isolated populations in the high altitudes of Lebanon (>2000 m). Moreover, a separation of the northern Turkish population from the southern Turkish populations was observed using both markers.
CONCLUSIONS: Morphological variation together with genetic and biogeographic studies make an effective tool for detecting relict plant populations and also populations subjected to more intensive selection.

Introduction

Juniperus excelsa M.-Bieb. (Grecian juniper) is an arborescent juniper that can reach 20–25 m in height (Farjon 2005, 2010; Schulz ; Adams 2008). It is slow growing, monoecious or dioecious, and wind pollinated (Farjon 2005; Adams 2008), with seeds dispersed by gravity or at longer distances by birds and small mammals (Jordano 1992; Santos ). It is a pioneer species, light demanding, with a high resistance to severe drought, cold conditions and shallow, degraded soils (Zohary 1973; Browicz 1982; Mayer and Aksoy 1986; Quézel and Médail 2003; Magyari ; Ozkan ).

Juniperus excelsa is divided into two subspecies based on morphological data (Farjon 2005, 2010): J. excelsa subsp. excelsa, covering mountain and sub-mountain areas from the Balkan Peninsula in the west, through Anatolia, Syria and Lebanon to Crimea in the north and Iran in the east (Jalas and Suominen 1973; Browicz 1982; Greuter ; Boratyński ; Christensen 1997; Farjon 2005, 2010) (Fig. 1); and J. excelsa subsp. polycarpos (K. Koch) Takht., found further to the east with a Transcaucasian-Central-Asian distribution. Adams (2008), based on random amplified polymorphic DNA (RAPD) molecular markers, considers these two taxa as separate species, J. excelsa and J. polycarpos, respectively.

Fig. 1

Inter-population diversity of The geographical positions of the sampled populations are indicated on a global distribution map of the taxa1,2,3 (acronyms as in Table 1). The K-mean analysis assigned the populations to three groups (1–3). The three main barriers (a–c) as obtained by Barrier 2.2 are also shown. 1Jalas and Suominen (1973); 2Browicz and Zieliński (1982); 3Boratyński ).

Table 1

Sampled populations of J. excelsa subsp. excelsa.

Code	Country	Locality	Longitude	Latitude	Altitude (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 Njass	N34°19′49″	E36°03′16″	1870–2300
LB4		Jbab El Homr	N34°20′16″	E36°12′18″	1860–2061
LB5		Barqa	N34°11′48″	E36°08′15″	1600–2200
LB6		Afqa	N34°04′25″	E35°54′20″	1100–1600
TU1	Turkey	Ilgaz-Tosya	N40°53′04″	E33°42′24″	850
TU2		Eğirdir	N38°08′12″	E30°46′42″	950
TU3		Göltarla	N36°34′56″	E29°58′29″	1100
TU4		Akçali Dağlari	N36°19′36″	E33°00′44″	1200
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″	35
CR2		Crimea-Kolkhoznoe	N44°29′00″	E33°49′54″	500

Juniperus excelsa subsp. excelsa is a major mountain forest element in the East Mediterranean Basin and sub-Mediterranean region. It colonizes sites that vary from sub-humid to the adjacent semi-arid steppe zone of the Mediterranean region. The altitudinal range of J. excelsa subsp. excelsa is very wide. In the Anatolian peninsula, it is mainly found at elevations between 1000 and 1300 m, and in Lebanon between 1600 and 1800 m in the western and eastern slope of Mount Lebanon (Quézel 1973; Abi-Saleh ; Akman ; Quézel and Médail 2003). It forms the tree line in the East Mediterranean Basin with old, sparse populations reaching elevations of 2100 m in Greece, and some individuals can be found at elevations of 2700–2800 m in the Taurus (Quézel 1973; Abi-Saleh ; Akman ; Browicz 1982; Barbero ).

The regions of contemporary occurrence of J. excelsa subsp. excelsa are situated around the Pleistocene refugial areas of the tertiary floras in the East Mediterranean Basin (Comes 2004; Tzedakis 2004; Weiss and Ferrand 2007, ; Médail and Diadema 2009). It seems possible that the Grecian juniper survived the glacial periods of the Pleistocene in places close to its present localities. Unfortunately, the pollen of junipers was not determined to the species level in palynological studies (e.g. Elenga ; Eastwood 2004; Tzedakis 2004). This makes a direct analysis of species migration during the Pleistocene/Holocene temperature oscillations impossible. In spite of that, the occurrence of the species during the last glacial maximum (LGM) was confirmed by macro-fossils from Eastern parts of the Balkan Peninsula (Magyari ). This could reflect a certain level of stability in the Eastern Mediterranean Basin during Pleistocene climatic oscillations that favoured the conservation of a high level of genetic and probably also morphological diversity of tree species (Fady-Welterlen 2005; BouDagher-Kharrat ; Fady ; Fady and Conord 2010; Douaihy ).

Morphological data are important in the comprehension of life cycles, geographical and ecological distributions, evolution, conservation status, as well as species delimitation (Kaplan 2001). However, with the rapid rise and advancement of molecular techniques, the role of the morphological data in phylogenetic studies was put into question and has raised an ongoing scientific debate (Stuessy ; Jenner 2004; Lee 2004; Wiens 2004; Wortley and Scotland 2006). The phenotypic variation of plants does not always follow the genetic pattern of variation and diversity of plant populations. The lack of congruence between morphological and genetic diversity was reported (e.g. Smissen and Heenan 2010; Ayele ). On the other hand, partial congruence in the geographical patterns of genetic and phenotypic diversity has been described several times (e.g. Jang ; Ruisi ). This is also the case for juniper taxa like Juniperus thurifera (Gauquelin ; Barbero ; Adams ; Romo and Boratyński 2007; Terrab ) and Juniperus pheonicia (Adams 2008; Boratyński ).

The morphological variability of J. excelsa across its large geographical range has not been extensively studied. The previous studies were always carried out on a limited number of populations (Barbero ; Christensen 1997; Mazur ; Farjon 2005; Marcysiak ). Mazur ) showed that the multivariate differences between two populations of J. excelsa from Crimea and one from the Balkan Peninsula correlated with geographical distance.

A recent genetic study on a large geographical range of J. excelsa subsp. excelsa based on nuclear microsatellites (Douaihy ) showed a high level of genetic diversity within this taxon, with a clear clustering of populations into three centres. The most differentiated populations corresponded to old vestigial stands found at the tree line (>2000 m) in Lebanon. The lower-altitude Lebanese populations clustered together, separated from the populations from Turkey, Cyprus, Greece and Ukraine (Crimea).

The main aims of the present study are (i) to present a first extensive morphological investigation of J. excelsa subsp. excelsa in the East Mediterranean Basin and (ii) to verify the hypothesis that J. excelsa subsp. excelsa has an intra-specific differentiation at the morphological level that corresponds to the differentiation described using molecular markers.

Materials and methods

Plant material and characters studied

Plant material was sampled from 14 populations from Greece, Cyprus, Ukraine (Crimea), Turkey (Anatolia) and Lebanon (Fig. 1). The samples of cones and twigs from the last ramification were collected only from arborescent adults bearing seed cones and were gathered separately from the southerly exposed parts of the individuals, at a height of ∼1.0–5.0 m above ground level, as described by Mazur . Three hundred and ninety-four individuals of J. excelsa were examined. Around 30 individuals per population were examined except for the Lebanese population LB3, where only 18 individuals with cones were found (Table 1). Each individual from every population tested was represented by 10 cones and 10 pieces of twig. Some individuals in the analysed populations had <10 twigs.

Nine morphological characters and eight ratios were examined (Table 2) in order to assess the variation within and between populations. The measurements of the shoots were performed under a stereoscopic microscope of ×8 magnification with a scaled ocular (LN and ST), as described by Marcysiak . Cone dimensions (CL and CD; Table 2) were measured using a numerical calliper (0.01 mm) and the numbers of cone scale rows and cone scales (CSR and CSN) were counted manually. Seed dimensions (SL and SW) were measured automatically using the WinSEEDLE software and the number of seeds in a cone (SN) was counted manually.

Table 2

Average value with standard deviation (SD), minima, maxima and coefficients of variation of analysed characters of cones, seeds and leaves of J. excelsa subsp. excelsa.

Code	Character	Mean ± SD	Minimum	Maximum	Coefficient of variation	N
CSR	Number of cone scale rows*	4.00 ± 0.06	4	6	1.50	3957
CL	Length of cone (mm)	9.19 ± 0.83	6.3	13.5	9.03	3956
CD	Diameter of cone (mm)	9.46 ± 0.98	6.1	14.4	10.36	3956
CSN	Cone scale number	6.00 ± 0.57	4	10	9.50	3956
SN	Number of seeds	5.64 ± 1.07	1	13	18.97	3923
SL	Length of seed (mm)	4.75 ± 0.35	3.08	6.5	7.37	3914
SW	Width of seed (mm)	2.89 ± 0.28	1.8	4.87	9.69	3914
LN	Number of leaves per 5-mm apical section of ultimate lateral branchlet	22.73 ± 2.91	12	36	12.80	3287
ST	Thickness of the last ramification shoot with leaves	0.72 ± 0.07	0.3	1.05	9.72	3286
CL/CD	Ratio of length of cone/diameter of cone	0.97 ± 0.06	0.72	1.28	6.19	3956
CD/CSR	Ratio of diameter of cone/number of cone scale rows	2.36 ± 0.25	1.44	3.61	10.59	3956
CSN/CL	Ratio of cone scale number/length of cone	0.66 ± 0.06	0.37	1.17	9.09	3956
SL/SW	Ratio of mean length of seed/mean width of seed	1.66 ± 0.14	1.08	2.46	8.43	3914
CD/SN	Ratio of diameter of cone/number of seeds	1.84 ± 0.40	0.61	10.09	21.74	3922
SW/SN	Ratio of mean width of seed/number of seeds	0.57 ± 0.15	0.19	4.46	26.32	3913
CD/SW	Ratio of diameter of cone/mean width of seed	3.32 ± 0.39	1.9	5.58	11.75	3913
ST/LN	Ratio of thickness of the last ramification shoot with leaves/number of leaves on the 5 mm of the last ramification shoot	0.03 ± 0.01	0.01	0.07	33.33	3286

N, number of measurements for a particular character.

*The scales of juniper cones are decussate or exceptionally ternate, alternately arranged and forming four or six rows, respectively.

Statistical analysis

The symmetry and unimodality of the distribution frequency of the measured character were verified using Shapiro–Wilks' W-test to assess the possibility of conducting a statistical analysis (Tabachnik and Fidell 1996; Zar 1999; Sokal and Rohlf 2003). The main statistics (arithmetic means, standard deviation, coefficient of variation) were calculated for the particular characters for individuals and populations. The correlations between the measured characters were verified using Pearson's correlation coefficient to avoid the most redundant ones. Tukey's T-test was performed to verify the influence of particular characters on the differentiation between individuals within populations and between populations. The level of intra-population variation was assessed using Ward's agglomeration on the shortest Euclidean distances among individuals (Zar 1999; Sokal and Rohlf 2003).

The discrimination analysis was performed to identify the discriminate power of each character and to determine the inter-population variation (Tabachnik and Fidell 1996; Sokal and Rohlf 2003). A dendrogram of the closest Euclidean distances on the basis of the unweighted pair-group method using arithmetic averages was constructed to check the affinities revealed in the discrimination analysis (Zar 1999; Sokal and Rohlf 2003). The discrimination analysis was calculated on the characters obtained from the ratios, except for the ratio of cone diameter to the number of cone scale rows (CD/CSR), because of the stable characteristic cone scale rows. All the above-mentioned statistical analyses were performed using STATISTICA 8 (StatSoft Poland, Tulsa, USA).

The geographical distances between the populations were calculated using MapInfo (Pitney Bowes Software Inc., New York, USA). The geographical structure of the populations was further analysed using Monmonier's maximum difference algorithm, implemented in BARRIER 2.2 software (Manni ). This analysis reveals discontinuities in morphological differentiation in relation to the populations' geographical arrangement. The barriers can be interpreted as breaks between adjacent populations in their morphological construction.

A K-mean analysis was performed on the calculated characters. This analysis reveals the number of K-clusters, which optimally illustrates the differentiation between populations. A classification matrix was constructed to show the percentage of individuals from each population that could be properly included in particular groups and to assess the number of individuals that fell into another K-group (Sneath and Sokal 1973).

Results

Evaluation of characters

The distribution frequency of the examined characters was normal or only slightly left- or right-biased (data not shown), which enabled further statistical analyses.

The coefficients of variation (CV) of the characters at the population level varied between 1.5 and 33 % (Table 3). The average values of cone length (CL) and cone diameter (CD) were ≈9 mm. These two characters varied from 6.3 to 13.5 mm (CV = 9 %) and from 6.1 to 14.4 mm (CV = 10 %). The average value of the ratio CL/CD was almost equal to one and varied between 0.72 and 1.28 (CV = 6 %). The largest cones were found in the populations from southern Lebanon and Crimea (samples LB5, LB6 and CR2, respectively), while the smallest ones were from the high-altitude population in Lebanon (sample LB3) (Table 3). The majority of the cones measured had four scale rows (CSR) and only a few cones (20 % of the trees from TU2) had six scale rows. The cone scale number (CSN) varied between four and 10 (CV = 10 %), and averaged 6.00 ± 0.57. The mean number of seeds per cone (SN) was close to 6, with cones containing 1–13 seeds (CV = 19 %) (Table 2). As for the cones, the highest value of SN was found in the sample from the eastern side of Mount Lebanon (sample LB5) and the lowest value was found in the sample from the high mountains of this country (sample LB3). The seeds had a mean length (SL) and width (SW) of 4.75 ± 0.35 and 2.89 ± 0.28 mm, respectively. Seed length varied between 3.08 and 6.5 mm (CV = 7 %), and SW between 1.8 and 4.87 mm (CV = 10 %). The mean ratio SL/SW was 1.66 ± 0.14 with CV = 8 %. The seed dimensions did not show the same geographical trend as for the cone dimensions: the smallest seeds were found in CR1 and LB1 (Table 3). The ratios CD/SN and SL/SW had a high level of variation (CV = 22 and 26 %, respectively). The ratios of the cone dimension over the seed dimension averaged 3.32 ± 0.39 (CD/SW, CV = 12 %).

Table 3

Descriptive statistics of the measured morphological traits at the population level.

Population	N	CSR		CL (mm)		CD (mm)		CSN		SN		SL (mm)		SW (mm)		LN		ST (mm)
		M	CV	M	CV	M	CV	M	CV	M	CV	M	CV	M	CV	M	CV	M	CV
LB1	30	4.00	0.00	9.46	7.36	10.04	8.85	5.92	8.44	6.16	19.50	4.54	7.48	2.82	10.45	21.78	11.28	0.74	7.33
LB2	30	4.00	0.00	9.45	7.12	9.63	6.37	6.19	5.88	5.40	15.49	4.74	7.13	2.79	10.76	21.45	8.64	0.71	7.58
LB3	18	4.00	0.00	7.80	7.41	7.70	7.07	5.34	7.81	4.23	21.50	4.78	7.44	2.93	7.73	25.10	13.19	0.78	5.59
LB4	27	4.00	0.00	8.93	7.24	9.25	7.99	5.97	4.33	4.92	16.87	4.90	6.05	2.91	8.67	25.84	12.29	0.82	7.82
LB5	30	4.00	0.00	10.10	7.12	10.48	7.90	6.01	3.59	6.28	17.70	4.93	7.68	2.88	11.27	22.49	11.96	0.75	5.77
LB6	30	4.00	0.00	9.61	7.04	9.83	7.55	6.01	9.40	6.16	16.44	4.57	7.57	2.69	10.26	22.23	11.42	0.70	8.67
TU1	31	4.01	0.90	8.93	7.40	9.27	9.64	5.65	7.11	5.45	14.10	4.70	5.77	2.95	6.94	22.76	11.22	0.74	10.16
TU2	30	4.05	5.93	8.81	7.96	8.97	7.73	6.02	14.16	6.00	17.64	4.86	6.49	3.08	7.86	21.32	13.23	0.67	8.28
TU3	30	4.01	0.91	9.08	7.36	9.49	6.61	6.22	6.45	5.73	12.11	4.76	7.44	2.95	7.60	23.23	8.06	0.70	6.96
TU4	29	4.01	0.93	9.41	7.16	9.90	6.62	6.39	11.17	6.22	12.26	4.82	6.75	2.81	7.78	23.53	8.87	0.66	8.53
GR	32	4.00	0.00	9.09	8.69	9.59	9.48	6.23	11.79	5.99	12.29	4.69	5.65	3.02	8.27	20.25	10.89	0.72	7.00
CY	30	4.00	0.00	8.72	6.67	8.46	6.99	5.55	9.63	5.46	16.25	4.86	6.98	2.77	10.14	24.11	12.69	0.70	6.91
CR1	30	4.00	0.00	8.75	5.85	8.84	7.47	5.91	5.39	4.60	18.18	4.47	7.73	2.84	7.66	22.63	13.82	0.69	9.30
CR2	30	4.00	0.00	9.87	7.02	10.23	6.48	6.27	7.13	5.84	17.21	4.97	5.68	3.07	8.02	23.04	11.33	0.71	6.41

N, number of individuals representing population; M, average value; CV, coefficient of variation (%).

The average number of leaves per 5-mm apical section of the ultimate lateral branchlet (LN) was 22.73 ± 2.91, and varied between 12 and 36 (CV = 13 %) (Table 2). The average value of the thickness of the last ramification shoot with leaves (ST) was 0.72 ± 0.07 mm and varied between 0.3 and 1.05 mm (CV = 6 %). All the populations had almost the same mean ST except for LB3 and LB4, which had a higher mean of 0.8 mm (Table 3). The ratios ST/LN, SW/SN and CD/SN had the highest levels of variation.

Cone length and cone diameter (CL and CD) are highly and significantly positively correlated (R2 = 0.98, P< 0.01), as well as seed length (SL) and seed width (SW) (R2 = 0.74, P< 0.01) (Table 4). The two features of the cone dimensions (CL and CD) significantly affected the number of seeds per cone (SN) (R2 = 0.77, P< 0.01 and R = 0.8, P< 0.01, respectively). Likewise, CL and CD were positively correlated with the number of cone scales (CSN) (R = 0.83, P< 0.01 and R= 0.86, P< 0.01, respectively). We obtained a high, but less significant, correlation coefficient (P< 0.05) between the number of seeds per cone (SN) and the number of cone scales (CSN). This correlation was most likely derived from the higher correlations noted earlier between these two characters and the cone dimension features.

Table 4

Correlation coefficients between nine characters of J. excelsa subsp. excelsa from all populations sampled; character acronyms as in Table 2.

Character	CSR	CL	CD	CSN	SN	SL	SW	LN
CL	−0.64*
CD	−0.55*	0.98**
CSN	−0.6*	0.83**	0.86**
SN	−0.16	0.77**	0.8**	0.56*
SL	−0.63*	0.53	0.43	0.4	0.24
SW	−0.56*	0.38	0.36	0.38	−0.01	0.74**
LN	−0.23	−0.13	−0.17	0.02	−0.32	0.58*	0.52
ST	−0.39	0.01	0.04	0.04	−0.36	0.32	0.52	0.53

*Significance at P< 0.05, **significance at P< 0.01.

Tukey's post-hoc test [Additional Information—File 1] showed that the highest number of statistically significant differences between the populations was found for the ratio of cone diameter/seed width (CD/SW), the ratio of cone diameter/number of cone scale rows (CD/CSR) and cone diameter (CD) and cone length (CL). The ratio of cone length/cone diameter (CL/CD), cone diameter/number of seeds per cone (CD/SN) and the seed length (SL) were significantly different between only a few populations, but mainly between the Turkish population (sample TU2) and some of the others. On the other hand, seed width (SW) showed significant differences only between Crimean and Turkish populations (CR2 and TU2). Moreover, the thickness of the last ramification shoot with leaves (ST), the number of leaves (LN) and the ratio of cone diameter/number of cone scale rows (CD/CSR) only differed significantly (P< 0.01) between two or three of the sampled populations.

All the calculated characters used in this study significantly discriminated between the samples at the level of P< 0.01. The discriminatory powers of the characters were very close, with values of partial Wilks' λ varying between 0.8 and 0.9 (Table 5).

Table 5

Discriminant power testing for the calculated characters of J. excelsa subsp. excelsa.

Character	Partial Wilks' lambda	P
Ratio of length of cone/diameter of cone (CL/CD)	0.904	0.000
Ratio of cone scale number/length of cone (CSN/CL)	0.848	0.000
Ratio of mean length of seed/mean width of seed (SL/SW)	0.741	0.000
Ratio of diameter of cone/number of seeds (CD/SN)	0.840	0.000
Ratio of mean width of seed/number of seeds (SW/SN)	0.817	0.000
Ratio of diameter of cone/mean width of seed (CD/SW)	0.800	0.000
Ratio of thickness of the last ramification shoot with leaves/number of leaves on the 5 mm of the last ramification shoot (ST/LN)	0.892	0.000

Intra-population variation

Sampling within each population was relatively uniform, with comparative numbers of individuals (Table 1). The within-population differentiation using Ward’s agglomeration method was generally found at similar levels, with the first split into two groups of individuals at a distance between 11 and 17 [Additional Information—File 2]. We note that the smallest distances for this split are observed in the most fragmented and marginal populations, such as LB3, CY, CR1 and GR. The maximum separation distances of individuals within populations varied between 6 and 9, and were observed for individual 22 from TU3, 13 from CY and 3 from GR.

Inter-population variation

In the discrimination analyses using the calculated characters (without CD/CSR), all of the individuals were grouped together, except for the high-mountain Lebanese population (LB3) and the Cypriote population (CY) (data not shown). The discrimination analysis at the population level shown in Fig. 2 reveals a clustering of the populations. According to the variable U1, which explained ∼53 % of the total variation and mostly depended on the ratio of cone diameter/seed width (Table 6), the high-altitude Lebanese population of Wadi El Njass (LB3) was separated from the other Lebanese populations (LB1, LB2, LB5 and LB6). The second high-altitude Lebanese population from Jbab el Homr (LB4) was dislocated at a middle distance in between. The populations most closely grouped were from Qammouaa, Barqa and Afqa in Lebanon (LB1, LB5 and LB6, respectively). All the other compared populations formed one group, although more dispersed.

Fig. 2

Discrimination analysis results for The results obtained by the two main discriminant variables (U1, U2) based on seven ratios are shown (acronyms as in Table 1).

Table 6

Coefficients of determination between discrimination variables and analysed characters of J. excelsa subsp. excelsa (character acronyms as in Table 2). Bold values correspond to the highest coefficient values.

	U1	U2	U3
CL/CD	3.65	1.74	0.07
CSN/CL	3.08	3.33	0.51
SL/SW	0.13	7.40	1.82
CD/SN	0.42	0.72	0.57
SW/SN	7.74	0.16	1.35
CD/SW	32.89	3.78	0.35
ST/LN	0.42	0.45	3.45

According to the variable U2, which was responsible for ∼20 % of the total variation and mostly depended on the ratio between seed length and width (SL/SW), all the populations except for three of the Turkish populations and the Greek one (TU1, TU2, TU3 and GR, respectively) formed one group. However, the separation of the high-altitude Lebanese population (LB3), and the populations from Cyprus and one from Turkey, can be recognized in the space between the second and third discrimination variables U2 and U3. These two variables also differentiated between all Lebanese and Turkish, and Greek and Crimean populations [Additional Information—File 3] (Table 6).

According to the cluster analysis by the Ward method, the sampled populations could be divided into three main sub-clusters (Fig. 3). The first included the high-altitude population from Lebanon (LB3) with the coastal population from Crimea (CR1) and the population from Cyprus (CY). The second cluster grouped the southern Turkish populations (TU2, TU3 and TU4) with the population from Greece (GR). The third cluster included all the other populations from Lebanon (LB1, LB2, LB4, LB5 and LB6) with the northern Turkish population (TU1) and the mountain population from Crimea (CR2).

Fig. 3

Dendrogram constructed by the Ward method of cluster analysis on the Euclidean distances between samples of ).

The K-grouping method also revealed the most probable split of the populations into three groups (Fig. 1). The classification matrix (Table 7) showed that the level of congruence was the highest for the Lebanese high-altitude, the Cypriote and the coastal Crimean populations (LB3, CY and CR1, respectively), where 67, 63 and 57 % of the individuals, respectively, were included in the correct population group (Table 7).

Table 7

Classification matrix for individuals of J. excelsa subsp. excelsa as a result of a K-means cluster analysis for calculated characters; acronyms as in Table 1.

Population	Level of conformability (%)*	Number of individuals classified to the populations
		A/N	LB1	LB2	LB3	LB4	LB5	LB6	TU1	TU2	TU3	TU4	GR	CY	CR1	CR2
LB1	31	29/30	9	3	0	0	2	4	2	0	0	0	6	0	2	1
LB2	33	30/30	2	10	0	2	0	3	0	0	1	3	1	0	4	4
LB3	67	18/18	0	0	12	0	0	0	0	0	0	0	0	6	0	0
LB4	36	25/27	1	1	1	9	0	0	1	1	0	3	1	3	2	2
LB5	36	28/30	2	4	0	0	10	7	0	0	0	1	0	0	0	4
LB6	47	30/30	2	3	0	0	4	14	0	0	3	1	0	1	1	1
TU1	38	29/31	4	0	1	1	2	0	11	2	2	1	1	2	2	0
TU2	33	24/30	0	0	0	0	0	2	6	8	5	0	1	1	0	1
TU3	23	30/30	2	1	0	2	0	1	2	2	7	3	4	1	1	4
TU4	41	29/29	0	1	0	0	2	2	1	0	7	12	0	1	0	3
GR	38	29/32	0	1	1	1	0	0	4	1	3	1	11	0	2	4
CY	63	27/30	0	1	1	1	0	1	0	3	0	0	0	17	2	1
CR1	57	30/30	0	2	0	2	0	0	0	3	3	1	0	1	17	1
CR2	37	30/30	1	2	0	1	5	0	3	0	2	2	2	0	1	11

*Percentage of individuals which were correctly classified to the population; A, number of analysed individuals; N, total number of individuals.

The first three barriers revealed by Monmonier's maximum difference algorithm applied by using BARRIER 2.2 software confirmed the separation of the high-mountain Lebanese population from all of the others and in the second level the differentiation of the Cypriote population (Fig. 1).

Discussion

Variation of characters

The values of the J. excelsa subsp. excelsa characters placed them directly within the characteristics of the taxon and did not differ drastically from the data reported in published taxonomic work and basic floras on samples from Bulgaria, Greece, Turkey and Crimea (Coode and Cullen 1965; Barbero ; Christensen 1997; Mazur ; Farjon 2005; Marcysiak ) [Additional Information—File 4].

We found that the values of the dimensional characters of the cones and seeds were positively correlated at a statistically significant level (P< 0.01). A high correlation between the cone dimensions (CL and CD) as well as between the cone and seed dimensions was also observed in previous studies on Juniperus species (Klimko , 2007; Mazur , 2010; Marcysiak ). On the other hand, we obtained a much higher correlation of seed number with the cone dimensions (CL and CD) and the CSN compared with the previous studies.

The characters that significantly differentiated between the largest number of populations were CD/SW, CD, CL and SN (Fig. 2). The latter two differed between the high-altitude population from Lebanon (LB3) and the other populations. The smallest cones found in that population, which had CL and CD <7 mm and the lowest number of seeds, with an SN of 4.2 on average, can be interpreted as resulting from the following three factors: (i) The harsher environmental conditions at the altitudinal line of the species. But we note that the mean cone dimensions in the other high-altitude population from the eastern side of Mount Lebanon (LB4) were as big as those in the middle-altitude populations. (ii) The marginal position of this population in the species range, as was observed in other taxa (Yamada and Miyamura 2005). However, the biggest cones were also found in the marginal populations from southern Lebanon and Crimea (samples LB5, LB6 and CR2, respectively). (iii) A genetic differentiation of this population.

The ontogenesis of J. excelsa seed cones has not been studied in detail. Our results on the mature cones indicate that typical dimerous cones strongly prevailed in all compared populations. The trimerous cones were found sporadically in three populations from Turkey. This pattern of cone arrangement is also typical for J. phoenicea (Schulz ; Mazur ) and J. oxycedrus (Klimko ). The number of seeds per cone (SN) reported for J. excelsa subsp. excelsa ranged between four and six, and sometimes from two to eight [Additional Information—File 4]. We obtained an average value of SN close to six, but this ranged from one to thirteen in particular populations. An average seed number higher than six for some populations (Table 3) suggests a more frequent occurrence of two whorls of seed scales, one with two ovules and the second with one ovule on each scale, or alternatively three ovules on one dimerous whorl of seed scales, as was described for other junipers (Schulz ). The two unexpected seed cones containing 13 seeds, found in populations from Greece and Turkey (GR and TU2, respectively), were dimerous. This suggests at least three whorls of fertile cone scales with at least one of them bearing three ovules.

Intra-population variability

The level of differentiation between individuals within populations, evaluated using Ward’s agglomeration analysis [Additional Information—File 2], was generally similar, except for the populations from the margins of the species’ geographical range, such as GR, CY and CR1, and also LB2 and LB3. The lower level of individual differentiation within these populations can be explained as resulting from (i) the long-lasting geographical isolation and adaptation to the local environmental conditions in Lebanese (LB2 and LB3) and Cypriote (CY) populations, (ii) the lower number of individuals representing Lebanese population LB3 or (iii) the possible origin from the low number of founders, such as in the case of samples GR and CR1. The latter two populations were also characterized by a lower level of observed heterozygosity than the other populations (Douaihy ), which could be a trace of an ancient founder effect.

We did not find a connection between genetic and morphological levels of variation when comparing the genetic heterozygosity with coefficients of variation of studied morphological characteristics of particular populations. This confirms the independence of the genetic and morphological markers, as reported earlier for several vascular species (Smissen and Heenan 2010; Ayele ).

Inter-population variability

The morphological multivariate differentiation of the studied populations of J. excelsa subsp. excelsa did not show a clear geographical pattern. The clustering of populations based on molecular data (Douaihy ) had a clearer geographical pattern with a strong clustering of the middle–altitude—1000–1900 m—Lebanese populations (LB1, LB2, LB5 and LB6) separated from the Turkish, Greek, Cypriote and Crimean populations that were grouped together (compare Fig. 2 in our paper with Douaihy : fig. 3). Nevertheless, some main congruencies are found between the two studies.

The discrimination, clustering and barrier analyses results based on morphological data (Figs 2, 3 and 1, respectively) showed a high differentiation of the high-altitude Lebanese population (LB3). This same population, along with another high-altitude Lebanese population from Aarsal (not included in this study), were found to be the most differentiated populations by Douaihy . Moreover, the separation of the northern Turkish (TU1) from the southern Turkish populations obtained by the Ward dendrogram based on the morphological data was also observed based on molecular data (Douaihy : figs 2–4).

The high level of morphological differences between geographically close populations in the Lebanese mountains could be interpreted as resulting from different abiotic conditions. However, genetic data also showed a separation of the high-altitude populations (Douaihy ). Hence, these differences could have resulted from a lack of gene flow and/or a different population history during glacial and postglacial migrations.

The junipers, as wind-pollinated plants, are expected to produce large amounts of pollen grain and could be successfully pollinated from long distances. Long-distance pollination plays a major role in reducing the differences among tree species populations (Hamrick ; Fady-Welterlen 2005; Fady ). The pollen grains of the juniper species are relatively small, characterized by slow setting velocity (Huntley and Birks 1983; Moore ), and their production is not very high. The dispersion of J. excelsa pollen grains in northern Iran was found to be significantly smaller than those of wind-pollinated broadleaved trees, such as Alnus, Carpinus, Quercus or Zelkova, and restricted to the areas close to juniper woodlands (Djamali ).

The similarity between Turkish, Crimean and Greek populations is very interesting (Fig. 2). The inclusion of CR1 and CR2 in the Turkish group of populations could have resulted from a common origin, which seems to be very possible when the history of the biota around the Black Sea during glaciation is taken into account. The level of the Black Sea was much lower during the glacial period than it is at present (Yena , 2005), making plant migrations along the coast from Anatolia to Crimea, and vice versa, much more possible.

The situation of the Greek population, marginal in the group of Turkish–Crimean populations (Fig. 2), could also have resulted from another palaeohistory event. The presence of J. excelsa in the Eastern Mediterranean Basin during the LGM (Magyari ) indicates possible migrations between the Anatolian and Balkan peninsulas during the glacial periods of the Pleistocene, and can explain the affinity between Turkish and Greek populations of this species. The different routes of migration to the Crimean and Balkan peninsulas can be a reason for the slightly higher differences between the populations of J. excelsa from these two centres, which also confirms the previous finding (Mazur ), and is partly consistent with the results of genetic analyses (Douaihy ). The J. excelsa populations differed significantly among each other but without a clear geographical clustering. A similar result was obtained in previous studies on: Juniperus seravschanica Kom. (Sultangaziev ), Juniperus oxycedrus L. subsp. oxycedrus from the Balkan Peninsula (Brus ) and for J. excelsa subsp. excelsa in the Lake district of the central Anatolia (Yücedağ ). The lack of a clear geographical clustering in these studies was explained by a restricted sampling area as well by a recent fragmentation that did not allow the appearance of a geographical structuring of the morphological traits.

The Cypriote population of J. excelsa subsp. excelsa differed morphologically from the geographically closest Turkish and Lebanese populations, but was genetically similar to the populations from Turkey, Crimea and the Balkan Peninsula (Douaihy ). This observation could have resulted from a variation in the environmental conditions or could suggest a selection process with a more rapid phenotypic than genetic differentiation in plants (McKay and Latta 2002). Interestingly, this was observed in the genus Cedrus, where the Cyprus population and the Turkish populations, phenotypically distinct, were found to be genetically very close (Bou Dagher-Kharrat ). We can similarly explain the morphological differences between core and marginal populations, such as Greek, Crimean and Lebanese.

Conclusions and forward look

The results of the multivariate biometric analyses follow at least partly the genetic differentiation of J. excelsa subsp. excelsa, also indicating the selection processes and/or variability in response to environmental differences in some localities, frequently marginal ones. It is concluded that examination of the morphological variation, together with genetic and biogeographic studies, should be treated as an effective tool for detecting relict plant populations and also populations with more intensive selection, where it is important to conserve the morphological characteristics. Finally, it should be stressed that plant populations with high levels of morphological diversity are key to establishing adequate strategies of biodiversity protection—crucial for conservation in the Mediterranean region. This is especially important in the case of relict, southern-most populations of J. excelsa in the mountains of Lebanon and in Cyprus, but also the northern-most ones in Crimea and on the Balkan Peninsula.

Additional information

The following additional information is available in the online version of this article –

File 1. Table. Results of the Tukey T-test for 17 characters of 14 samples of J. excelsa; *significance at level P= 0.05; **significance at level P= 0.01 (population acronyms as in Table 1; characters as in Table 2).

File 2. Figure. Intra-population variability of 14 populations of J. excelsa subsp. excelsa analysed using the Ward agglomeration on the shortest Euclidean distances among individuals (population acronyms as in Table 1).

File 3. Figure. Discriminate analysis results based on the three main variables (U1, U2 and U3).

File 4. Table. The morphological variation between two populations of J. excelsa subsp. excelsa from the Crimea and one from the Balkan Peninsula.

Sources of funding

The research was conducted within the statutory research programme of the Institute of Dendrology, Polish Academy of Sciences, Kórnik, Poland, and supported by the Research Council of Saint-Joseph University of Beirut.

Contributions by the authors

All the authors contributed to a similar extent overall.

Conflict of interest statement

None declared.