Effect of geographic range discontinuity on taxonomic differentiation of Abies cilicica

Three populations of Abies cilicica subsp. isaurica and four of A. cilicica subsp. cilicica were analyzed using 35 morphological and anatomical needle characters with the implementation of multivariate statistical methods to verify the differences between subspecies. Moreover, the possible geographic differentiation of A. cilicica subsp. cilicica populations from the East Taurus and Lebanon Mountains was examined. Abies cilicica subsp. isaurica has been distinguished from A. cilicica subsp. cilicica by its glabrous young shoots and resinous buds. We detected that needles of A. cilicica subsp. isaurica are longer, broader and thicker, with a higher number of stomata rows, and larger cells of the epidermis, hypodermis and endodermis than A. cilicica subsp. cilicica. Additionally, A. cilicica subsp. isaurica needles have frequently rounded to obtuse-acute apex and resinous canals positioned more centrally inside the mesophyll than needles of A. cilicica subsp. cilicica. This indicates that a set of most of the tested needle characters can be used to distinguish the subspecies; however, any of characters enable that when used separately. Morphological and anatomical distinctiveness between these two taxa justify their recognition at the subspecies rank. Additionally, the populations of A. cilicica subsp. cilicica from the East Taurus and Lebanon are morphologically different. This geographic differentiation of populations is congruent with results provided by genetic analyses of nuclear microsatellites markers (nSSR).


Introduction
The observed geographic ranges of species are historically determined and have been formed together with the species evolution [1]. The Mediterranean region history has been altered with geological events and climate changes. The land movements connected with regression of Thetys [2][3][4][5], the Messinian "salt crisis" [6] and climate cooling during late Tertiary and Quaternary, with the Pleistocene climate oscillations [7][8][9] had a significant imprint on the plant evolution and migrations. These processes also concerned the oro-Mediterranean plant species [10], which evolved together with the formation of mountain ridges [8].
The East Mediterranean mountain systems have been formed mostly during Miocene [3,4]. Expansion of ancestors of the genus Abies is connected with this process [11]. The ancestor of contemporary A. cilicica (Antoine & Kotschy) Carrière probably appeared during Oligocene and Miocene [12] and settled first at Taurid and then also at Lebanese mountains [11]. It also had a somewhat broader geographic range during Miocene-Pliocene than at present ( [11] and Fig. 2 and Fig. 3 therein). The Pliocene climate cooling and Pleistocene climate oscillations were the reasons for the fragmentation of the geographic range of A. cilicica, including its divergence and the formation of the subspecies A. cilicica subsp. isaurica Cullen & Coode in the West Taurus [13]. The development of the Taurids in Anatolia and at the Lebanese mountains allowed A. cilicica to persist in these regions during Pleistocene. The species could migrate up during hot and down during cold periods [7,8]. However, the isolation of the mountain massifs and, more importantly to A. cilicica, the climate aridity during cold periods [14], has led to the reduction and strong fragmentation of the geographic range of the species. The early Holocene distribution of the genus Abies in the Mediterranean region was more abundant than at present [15]. The increased aridity during late Holocene together with extensive deforestation from the millennia concerning this region [16,17] formed the bases of the further fragmentation of oro-Mediterranean tree species, including the Cilician fir [18][19][20][21]. Finally, the historically broader geographic range of A. cilicica has been reduced [22] and the species is currently at risk of extinction due to aridity in its lower localities [21]. It was recognized as a near threatened species in Turkey, Syria and Lebanon [18,20,23]. It grows in the areas assumed to be glacial refugia of the Tertiary flora [24] in several dozen mountains isolated from each other [20,25].
The spatial isolation between the West (Isaurian) Taurus and East Taurus is assumed to be one of the reasons for the differentiation of A. cilicica into eastern A. cilicica subsp. cilicica and western A. cilicica subsp. isaurica, with pubescent versus glabrous young shoots, respectively [26][27][28]. These two subspecies were clearly distinguished using nSSR markers [13]. The disjunctive character of occurrence of the Cilician fir and genetic differentiation between populations form the West and East Taurus and Lebanon Mountains also suggests morphological and anatomical differences between them. Similar geographic pattern of phenotypic structure has been described for Juniperus excelsa M. Bieb. using morphological characters of cones and sprouts [29] and for Cedrus libani A. Rich. using morphological and anatomical characters of needles [30].
Thus, we hypothesized that (i) the long lasting spatial isolation between West Taurus and East Taurus and Lebanon mountains caused not only genetic but also morphological and anatomical differences between A. cilicica subsp. cilicica and A. cilicica subsp. isaurica, and (ii) the isolation between mountain massif within the geographical range of A. cicilica subsp. cilicica also involved phenotypic differences between populations from these distant regions. In the present study we verified these hypotheses applying biometric analyses of morphological and anatomical needle characters. Most of the characters used in our study are applied for the first time and were not considered before in the subspecies descriptions (except of L, MW, SW and ST; see [31]), nor for evidence of phenotypic differentiation in the geographical space.

Studied species
Abies cilicica is a large tree, attaining a height of 30-35(-42) m and diameter of 1(-2) m at 1.3 m above ground [22,27,32]. It grows in the mountains of the East Mediterranean region, in Turkey in the West and East Taurus and in the Amanos, in Syria on the Jebel Ansariye and in Lebanon on the J. Ammoua and the J. Ehden [18,19,25] (Fig. 1). In Turkey, A. cilicica occurs between 1150 m and the timberline at about 2000 m on the north facing slopes, and between 1450 and 2000 m on the south facing slopes of the Taurus, with optimal conditions between 1200 and 1800 m, mostly in the valleys [20,22]. The species forms pure, shady forests or mixed forests with Cedrus libani, and also with Pinus nigra J.F. Arnold subsp. pallasiana (Lamb.) Holmboe in the West Taurus [20,22,33]. Juniperus excelsa, J. foetidissima Willd. and J. drupacea Labill. frequently enter the Cilician fir forests and even replace A. cilicica when overexploited and/or overgrazed [33].

Plant material
The needles were collected from seven natural populations of A. cilicica, four representing A. cilicica subsp. cilicica from the East Taurus Mountains in Turkey and the Lebanon Mountains, and three representing A. cilicica subsp. isaurica from West Taurus in Turkey (Tab. 1). Thirty cone-bearing individuals, separated by a distance of about 50 m, were sampled from each population, with the exception of the LB_2, where only 12 individuals could be sampled. Studied individuals were ascribed to the subspecies based on morphology of the young shoots [26,27] and molecular identification [13]. Ten needles from the central part of a two-year-old shoot increment were collected from each individual, from the sunny, predominantly south-facing parts of the tree crown, at a height of about 2.0 to 5.0 m above ground level. Plant material was conserved in 70% alcohol and kept there until further preparation and measurements. In total, 192 individuals were examined, represented by 1920 needles.
Five needles from each individual were used to analyze morphology, and another five to measure anatomical characters from needle cross-sections (Tab. 2). The set of biometric characters, methods of preparation and measurements were based on previous investigations of West Mediterranean firs [34] and Turkish firs [31], and supplemented by characters of stomata occurring on the upper side of needles (Tab. 2). The CH was estimated in the following scale: discontinuous layer of single cells -0.5; continuous layer of single cells -1; continuous layer of single cells with additional discontinuous cell layer -1.5.

Statistical treatment
The normality of the frequency distribution of each character was verified using the Shapiro-Wilk W-test, and the homoscedasticity of the variance of the measured data using the Brown-Forsythe test. The evaluated characters (LSU, TA and LC) data were converted to percentages and arcsine transformed. Values of all characters were standardized before multivariate statistical analyses [35]. The Pearson correlation between characters was verified to avoid the most redundant ones, with |r| > 0.9.
A t-test (measured and ratio characters) and the Mann-Whitney U-test (evaluated characters) for independent samples were used to evaluate the significance of differences between the subspecies of A. cilicica and between Turkish and Lebanese populations of A. cilicica subsp. cilicica. Tukey's honest significant differences (HSD) post-hoc test and Kruskal-Wallis test for the characters with biased distribution were performed on average values of characters for individuals to test the significance of differences between populations, and, consequently, between subspecies and regions.
A forward stepwise discrimination analysis (FSDA) was performed to identify the discrimination power of each character, to eliminate the closely redundant ones and to detect the relationships between populations, and consequently between subspecies and regions. A set of cluster analyses on the shortest Euclidean distances and Mahalanobis' distances (after Ward's, UPGMA, WPGMA) were applied to verify the relationships between populations between taxa and regions. Afterwards, it was verified again using discrimination analysis, to detect fit differentiation of particular individuals from the populations representing each of the groups [35]. The statistical analyses were carried out using STATISTICA v. 9 (StatSoft PL).
The Mantel test [36] was implemented to verify the relationships between Euclidean distances among populations and the geographic distances. Geographic distances were retrieved from the geographic coordinates, using MapInfo 9.5 (Pitney Bowes). The significance of the correlation was tested with 9999 random permutations. PopTools v.3.2 software [37] was used in the calculations. N -number of individuals sampled; AMT -annual mean temperature; APR -annual average precipitation.

Variation and correlation of characters
The distribution of most of the characters was unimodal and normal or very close to normal. The evaluated characters LSU, TA and LC were the only exceptions. The latter data were arcsine-transformed and assumed to have a closeto-normal distribution, which allowed the application of multivariate tests. The data after transformation and standardization were homoscedastic or close to, which allowed the assessment of parametric tests.
The needle dimensional characteristics (A, P, and L) correlated positively with each other at very high significant level (r = 0.95, P < 0.01). The anatomical characters of the needle ST, SW, VCT and VCW, as well as WC and HC correlated significantly with each other at a similar level. The level of correlation was slightly different for each population, but generally the same pattern of relationships between measured characters was found. From the groups of the most closely correlated and thus redundant characters, only single ones were used for the multivariate analyses. The forward stepwise analysis of discrimination (FSDA) reduced the set of characters and only 22 from 48 previously measured/evaluated ones were the basis of the discrimination and clustering, which described the differentiation between populations, subspecies and regions. The fourteen needle characteristics discriminated between populations of A. cilicica s. l. at a significant level (P < 0.01; Tab. 2), but P, MT, NC, MC, LSU_1, LSU_2, TA_3, LC_3 and LC_45 were excluded from the dataset in the FSDA. The highest discriminant power had NSL, LC_1, RW_1 and A with values of partial Wilks' λ of 0.620, 0.811, 0.836 and 0.844, respectively.
The particular characters differed in the value of variation coefficients. NC was the most stable trait, completely without variation in several populations and V = 0.4% on average. Among the other characters, NSL, EW and HH had average values of V ≈ 7%. Apex forms (TA), position of resin canals (LC) and location of stomata on the upper side of the needle (LSU) were the most variable. Among the measured characters, A, P, L, BD, DV and CH had V between 20 and 40% (Tab. 2).

Phenotypic distinctiveness of subspecies
The average values of characters appeared to some degree to be specific for particular populations, but with generally overlapping frequency distribution between populations. Most of the analyzed needle characters differentiated between A. cilicica subsp. cilicica and A. cilicica subsp. isaurica at a statistically significant level (Tab. 2). The average values of needle characters were higher in populations classified as A. cilicica subsp. isaurica than in populations of A. cilicica subsp. cilicica, with the only exceptions being DV and NC. The ratio characters differentiated subspecies to a lesser extent, with only RW_1, RW_2, LMW and HS being significant at P < 0.01 (Tab. 2).
According to the Mann-Whitney U-test, significant differences (P < 0.01) between subspecies were observed in the TA_1 and LC_1 (Tab. 2). The HSD Tukey's test also revealed that the majority of morphological and anatomical characters differed at a significant level (P ≤ 0.01) between sampled populations of A. cilicica subsp. cilicica (TR_1, TR_2, LB_1 and LB_2) and A. cilicica subsp. isaurica (TR_3-5; Tab. 3). The characters DV, NC, LSU_3, TA_3, TA_5, LC_2, LC_3 and LC_4 were the only biometric characters that did not differ significantly between samples according to the results of the Tukey's test.
Based on the first three discriminant variables of FSDA, U 1 , U 2 , and U 3 , which explained 86% of the total variation, the analyzed populations formed three groups. The first group was composed of A. cilicica subsp. isaurica populations (TR_3, TR_4 and TR_5), while the other two were of A. cilicica subsp. cilicica (Fig. 2a-c). The first discrimination variable (U 1 ), which covered 60% of the total variation, was determined mostly by LC_1, NSL and NRL, the second (U 2 ), which covered 15% of the variation, was determined first of all by NML, NS and RW_2, while the third (U 3 ), which covered 11% of the variation, was determined by RW_2, CH and TA_1. The analyzed populations were discriminated by U1 at the subspecies level (Fig. 2a), while further grouping of the populations of A. cilicica subsp. cilicica was mostly determined by U 3 (Fig. 2b and Fig. 2c).
Afterwards, we verified how particular individuals from the populations representing each of subspecies fitted this differentiation. Again, we used FSDA with the characters: NC, MC, CS, LSU_1, LSU_2, TA_3, TA_4, TA_5, LC_2, LC_3 and LC_4 excluded from the model. From the remaining characteristics, 11 discriminated between individuals at a significant level. NSL, LC_1, CH and RW_1 had the highest discrimination power, with partial Wilks' λ values: 0.800, 0.846, 0.883 and 0.907, respectively. The total variation was divided between the first two discriminant variables, where U1 covered more than 81%. It was determined first of all by NRL, LC_1, NSL, TA_1, LC_5 and A. The second discrimination variable U2 was determined mostly by CH, RW_2 and HS. The individuals formed three groups on the dispersion diagram (Fig. 2d). The populations representing A. cilicica subsp. isaurica (TR_3, TR_4 and TR_5) formed a coherent group with only one individual outside of the 95% confidence interval, but included six individuals from A. cilicica subsp. cilicica (Fig. 2d). In summary, 95% of individuals of A. cilicica subsp. isaurica were correctly classified to the subspecies.
The cluster analysis on the shortest Euclidean distances according to Ward's method, divided all of the samples into two main groups. The populations assigned to A. cilicica subsp. cilicica formed the first cluster, while the populations classified as A. cilicica subsp. isaurica (TR_3, TR_4, and TR_5) comprised the second one (Fig. 3). Similar patterns of differences between A. cilicica subsp. cilicica and A. cilicica subsp. isaurica populations were detected using UPGMA and WPGMA cluster analyses on the Euclidean distances and analyses on Mahalanobis distances (data not shown).

Variation within subspecies
All populations of A. cilicica subsp. cilicica had a marginallower type of resin canal position (LC_1), while two types of resin canal positions were observed at a similar frequency in A. cilicica subsp. isaurica, namely marginal-lower (TR_3 and TR_5) and mesophyll-central (TR_5; Fig. 1). This subspecies was quite homogenous in terms of the location of stomata at the upper needle surface (LSU), while A. cilicica subsp. cilicica was more variable in this aspect (Fig. 1). Individuals of A. cilicica subsp. cilicica had indented (TA_1), or rounded (TA_2) types of needle apex, 90.5% and 9.5%, respectively, while in A. cilicica subsp. isaurica all of the types were observed, with prominent percentages of obtuse (TA_3) and acute (TA_4) types (Fig. 1).
The Mantel test detected positive and significant correlations between the Euclidean distance and geographic distances for populations (r 2 = 0.54, P = 0.012). The multivariate differences were found not only between subspecies, but also between populations of A. cilicica subsp. cilicica from the East Taurus (TR_1 and TR_2) and the Lebanon mountains (LB_1 and LB_2; Fig. 2a-c). The latter two groups of populations were determined mostly by the U 3 variable ( Fig. 2b and Fig. 2c). We used FSDA to verify how particular individuals of A. cilicica subsp. cilicica from the East Taurus and the Lebanon mountains fit the two geographic groups described. The FSDA detected that the compared individuals formed two partly intermingled groups on the dispersion diagram (Fig. 2d). Three individuals from the East Taurus and another three individuals from the Lebanon Mountains fall into the 95% confidential interval of A. cilicica subsp. isaurica. The individuals of A. cilicica subsp. cilicica from the Eastern Taurus formed a separate group from that representing the Lebanon Mountains; however, about 30% of the East Taurus individuals entered the Lebanese group at a 95% confidential interval (Fig. 2d). The correct classification of the Lebanese versus East Taurus individuals were at the level of 93% and 86%, respectively. The populations of A. cilicica subsp. cilicica from East Taurus and the Lebanon Mountains differed with respect to the type of needle apex. In the Lebanese populations only the indented type (TA_1) was observed, while in those from the East Taurus 20% of individuals had rounded (TA_2) type of needle apex (Fig. 1). A significant level of statistical differences (P < 0.01) was also observed between the Turkish and Lebanese A. cilicica subsp. cilicica populations for number of stomata (NSL), width of hypodermal cells (HW), shape of needle cross-section (NS), shape of hypodermis cells (HS) and shape of resin canal cross-section (CS; Tab. 2). The geographic differentiation of A. cilicica subsp. cilicica, however, has not been confirmed using the agglomeration method (Fig. 3).

Needle characteristics variation
Data on the morphological and anatomical variation of the needle characteristics of A. cilicica were scarce, with only the length and width of needle (L and MW) and sometimes the needle apex type (TA) reported. This results in a low level of differences between the Mediterranean taxa of the genus Abies on the needle characters known to date [38,39]. It is commonly known and generally accepted that cones are essential to correctly determine the Abies taxa (e.g., [27,28,[40][41][42]). This rule was also confirmed in the only known biometric study of the Turkish firs, but some differences between Turkish fir species in the needle characters were also described [31]. Comparing these data with our findings, it should be stressed that we found higher average values of needle width and height on the cross-section preparation (SW and ST, respectively) and diameter of resin canals (WC and HC) than reported by Bagci and Babaç [31]. The differences between our data and that of Bagci and Babaç [31] might be a result of different preparation and measurement procedures used in both studies and the higher number of individuals tested in our study. The comparison of Bagci and Babaç [31] and other accessible data with our results also stresses similarities in the data concerning the length and width of needles (Tab. 4).
Our data are based on the examination of a large number of individuals and thus shall be considered as bearing not only the real values of the examined characters, but also ranges of variation. Our results fill the gap in data and provide a broad set of needle characteristics of A. cilicica. We expect that some of them could also be used in palaeobotanical studies. Abies needles were detected several times in the Tertiary and Quaternary deposits and many of them have not been determined to the species level (e.g., [12,40,43]).

Infraspecific differentiation
Our study is the first where the differences between A. cilicica subsp. cilicica and A. cilicica subsp. isaurica are studied using wide spectrum of the needle characters [only L and MW (e.g., [28,41,42,44,45]), ST, SW and WC/HC were investigated before [31]].The biometric analyses reveal that the most of examined needle characteristics are suitable to distinguish A. cilicica subsp. cilicica from A. cilicica subsp. isaurica at a significant level (Tab. 2).  The differences between average values of the most of verified characters found in our study justify the taxonomic position of A. cilicica subsp. isaurica when compared with typical A. cilicica subsp. cilicica. Generally, the needles of A. cilicica subsp. cilicica are smaller than those of A. cilicica subsp. isaurica (compare characters A, P, L, MW; Tab. 2), have a smaller endodermis tube (SW and ST), slighter epidermis and hypodermis cells (EW, EH, HW, HH), lower values of resin canal width and height (WC, HC) and lower numbers of stomata rows and stomata (NRL and NSL). Abies cilicica subsp. isaurica could be distinguished using a set of these characters and the evaluated ones, which are types of location of stomata on the adaxial needle side (LSU), the needle apex type (TA) and position of resin canals (LC; Fig. 1). The average values of measured characters of A. cilicica subsp. isaurica are about 20-30% higher than detected for A. cilicica subsp. cilicica, which has not been described until now [26,27,31,32]. However, none of the mentioned characters allows distinction between subspecies solely, as the distribution ranges of the characters that may be used for distinguishing between subspecies overlap to some degree.

Geographic pattern of differentiation
The Mantel test result suggests an important role of spatial isolation in shaping the inter-population differentiation of the phenotypic characters. The pattern of intraspecific morphological and anatomical differentiation of A. cilicica s. l. documented in the present study based on the needle characteristics appeared similar to those described using nuclear microsatellite markers (compare Fig. 1-Fig. 3 and Fig. 1, Fig. 2 in [13]). On the other hand, the geographic differentiation among populations of A. cilicica subsp. cilicica was less evident in phenotypic characters than in molecular markers. This result is somehow in contrary with molecular evidence, because the genetic differences between Lebanese and Turkish populations of A. cilicica subsp. cilicica were even at a higher level than between subspecies (see [13] and Fig. 2 therein).
The pattern and significant level of genetic differentiation found between populations of A. cilicica subsp. isaurica and A. cilicica subsp. cilicica as well as between populations of the latter from the Lebanon Mountains versus East Taurus were interpreted as a result of a long-lasting isolation [13].
The spatial isolation and climate changes during glacial and interglacial periods of the Pleistocene [46,47] caused only local, vertical migrations in the mountains of the Mediterranean region [47,48], which reduced the possibility of gene exchange by seeds and pollen between populations and, consequently, were the reason for differentiation or even speciation processes (e.g., [49][50][51]). This also concerns the A. cilicica (or its ancestor) populations in the Taurids and Lebanese mountain systems [11]. A similar pattern of morphological differentiation to that mentioned above was detected in Juniperus drupacea [52], J. excelsa subsp. excelsa [29] and Cedrus libani [30]. All three species co-occur with A. cilicica [22,25,32,53]. Interestingly, the geographic differentiation on the morphological and/or anatomical characteristics of each of these three taxa resembled geographic structure on the genetic markers. Congruent genetic and phenotypic patterns of differences between the West Taurus, East Taurus and Lebanon Mountains populations were detected in Cedrus libani [30,54] and Juniperus excelsa [55], and differences between Lebanese and Turkish populations were found in J. excelsa [29,55]. This could indicate a more universal character of differentiation that resulted from the species history and ancient demographic processes for the oro-East-Mediterranean tree species.

Conclusion
The populations sampled as A. cilicica subsp. cilicica and A. cilicica subsp. isaurica could be clearly distinguished, but only using the set of morphological and anatomical characters of the needles. No one single needles character allowed distinguishing between them without any doubt. The differences were also detected between Lebanese and East Taurus populations of the typical subspecies A. cilicica subsp. cilicica. The geographic pattern of differentiation among populations based on the morphological and anatomical needle characters resembles those received with the nSSR markers [13]. The geographic differentiation between both subspecies and among populations in the East Taurus and Lebanon Mountains, detected using both nSSR markers and phenotypic characters, suggests local management of the A. cilicica woodlands, without seed exchange between the regions.