Mediterranean Islands Hosting Marginal and Peripheral Forest Tree Populations : The Case of Pinus brutia Ten . in Cyprus

Mediterranean islands have served as important Tertiary and glacial refuges, hosting important peripheral and ecologically marginal forest tree populations. These populations, presumably harboring unique gene complexes, are particularly interesting in the context of climate change. Pinus brutia Ten. is widespread in the eastern Mediterranean Basin and in Cyprus in particular it is the most common tree species. This study evaluated genetic patterns and morphoanatomical local adaptation along the species geographical distribution and altitudinal range in Cyprus. Analysis showed that the Cyprus population of P. brutia is a peripheral population with high genetic diversity, comprised of different subpopulations. Evidence suggests the presence of ongoing dynamic evolutionary processes among the different subpopulations, while the most relic and isolated subpopulations exhibited a decreased genetic diversity compared to the most compact subpopulations in the central area of the island. These results could be the consequence of the small size and prolonged isolation of the former. Comparing populations along an altitude gradient, higher genetic diversity was detected at the middle level. The phenotypic plasticity observed is particularly important for the adaptive potential of P. brutia in an island environment, since it allows rapid change in local environmental conditions.


Introduction
The islands of the Mediterranean basin comprise one of the 36 terrestrial biodiversity hot spots of the world and are characterized by high diversity of landscape and vegetation types due to the complex historical biogeography and the profound environmental heterogeneity [1,2].Mediterranean islands contain a significant component of Mediterranean biodiversity, notably a number of range-restricted species and unusual vegetation types [3,4].The vegetation types usually considered as "typically Mediterranean" are the evergreen and sclerophyllous shrublands and forests under semi-arid or subhumid bioclimates, corresponding to thermo-mediterranean and meso-mediterranean vegetation belts [4].Although the majority of Mediterranean islands are "continental islands" (that is, they became progressively isolated from the mainland and from each other by a complex combination of tectonic and glacio-eustatic processes), Forests 2018, 9, 514 2 of 27 there are also cases of Mediterranean "oceanic islands" in the geological sense which were part of mainland (such as islands that emerged from the bottom of the sea) (see previous papers [3,5]).Thus, "oceanic islands" often present lower richness of biodiversity elements compared to "continental islands" [6], while the genetic background of their wild populations could be restricted, due to isolation, small population, founder effects, bottlenecks, low effective population sizes, and genetic drift [7].
The Mediterranean islands have served as important Tertiary and glacial refuges, and, hence, Mediterranean islands possess highly polymorphic species and vicariant endemic plants which emerged from more or less recent speciation events [6,7].The geographic isolation and the environmental heterogeneity of the Mediterranean basin have favored diverse evolutionary processes of gradual speciation of plants, such as genetic drift or adaptive radiation (see previous papers [4,8,9]).These features indicate the role of the wild populations of flora and fauna species in the Mediterranean islands as important peripheral (and marginal) populations.Currently, the value of peripheral forest tree populations is of particular interest in the context of climate change [10].These populations may concurrently be those where the most significant evolutionary changes will occur; those with increasing extinction risk; the source of migrants for the colonization of new habitats at leading edges; or the source of genetic variation for reinforcing existing genetic variation in various parts of the range [10].
Several authors argued that demographic and evolutionary processes shape peripheral populations differently, compared to populations at the core of the distribution, depending on their situation in the geographic space [11][12][13].Thus, whether leading edge populations are diverse enough to efficiently contribute to colonization will depend on their interpopulation gene flow and the amount of gene flow from core populations [10].The Mediterranean island forest tree populations are identified as geographically peripheral populations, since they are found at the rear edge of distribution areas; indeed an increasingly unfavorable climate may lead to their ecological marginalization, with drastic consequences for their survival [10].
Brutia pine (Pinus brutia Ten.), along with Aleppo pine (Pinus halepensis Mill.) form a distinct group (Group Halepensis) within the Eurasian hard pines; their combined geographic distribution reflects their prominence among low-elevation Mediterranean forest species [14].Pinus brutia Ten. is a coniferous species confined mainly to the Eastern Mediterranean region (incl.Greece, Turkey, Cyprus, Syria, and Lebanon), and can also be found in Iraq and Iran [15,16].P. brutia grows under several variations of the Mediterranean climate [15], on a wide range of soil types [17,18], while it is recognized for its adaptation to drought and alkaline soils [19,20].In addition, P. brutia is able to form stable vegetation associations with broad-leaved species [21,22], a characteristic which has led to an increased interest in the species for commercial plantations, illustrated by breeding and provenance trails carried out in various countries in the Mediterranean region [14,23,24] and even in wider geographical ranges [20,25,26].The P. brutia forests correspond to habitat type "9450 Mediterranean Pine Forests with Endemic Mesogean pine" according to the European Directive 93/42/EC.Despite the fact that the largest area covered by P. brutia occurs in Turkey (3.8-5.4 million ha) [16,27,28], the conservation interest nowadays focuses on the geographically peripheral and ecologically marginal populations of this species.
The present study examines the genetic and ecological processes acting on the peripheral population of P. brutia in Cyprus, an oceanic island as defined above, located at the southern edge of the species' distribution.Thus, the current study examines possible genetic and morphoanatomical responses of P. brutia in relation to (i) its geographical distribution within the island of Cyprus and (ii) its distribution in different altitudes.Therefore, this study is an attempt to pursue a genetic and morphoanatomic analysis at a landscape scale, following a refined longitudinal, latitudinal, and altitudinal sampling within an island environment.Its outcomes will be invaluable in delineating a conservation strategy for mitigating adverse effects of global warming on the potential growth and survival of P. brutia forests in Cyprus and their genetic resources.The outcomes from this study will also contribute in obtaining more knowledge on the evolutionary capacity of geographically peripheral and potentially ecologically marginal populations of coniferous species, under the island Forests 2018, 9, 514 3 of 27 environment (biogeography) and their adaptability to changing conditions, as well as the impact of altitude gradients on population genetic structure within the island ecosystems.

The P. brutia Forest in Cyprus
In Cyprus the thermophilous pine forests with P. brutia is the most extensive and widespread forest type, occurring in all mountainous areas from dry to subhumid climates (0-1.400m), covering 66% of island-wide forest land (~88,790 ha) [29].The Troodos range is well covered with dense pine forests, which attain their best development in Pafos forest (60,159 ha), where the largest unfragmented and best conserved P. brutia forests are found.According to the Rivas-Martínez bioclimatic classification, Cyprus has a Mediterranean Mesophytic to Xerophytic-Oceania bioclimate with zones ranging from Thermo-Mediterranean-semi-arid (lowlands) to supra-Mediterranean-humid (Troodos) [30].

Sampling Design
For the purposes of the current study, sampling was carried out in the central and main mountain range of the island, namely the "Troodos Mountain range".Sampling covered a distance of 90 km longitudinal and 45 km latitudinal (Figure 1) and was implemented in the six forests, as these are defined by the Department of Forests, namely: Akamas forest, Pafos forest, Troodos forest, Adelphi forest, Limassol forest, and Macheras forest (referred to as "subpopulations", see Table 1 & Figure 1).Plant tissues were collected from adult trees 50-70 years old (see Table 1) at a distance of 200 m in order to avoid genetic kinship.Within Pafos forest, sampling adopted ecological parameters, namely: (i) altitudinal subpopulations (altitude zones of 400 m, 800 m, and 1200 m-referred to as "altitudinal subpopulations") and (ii) different aspects, namely the northeast and the southwest orientation (referred to as "aspect subpopulations").P. brutia forest is the dominant vegetation in Pafos forest (size: 60,159 ha), shaping the best growing P. brutia forest on the island, covering a large area of forest from near-sea-level up to the peak of Tripylos at 1352 m.These characteristics allow a detailed assessment of genetic structure reflecting the dynamic effect of differential adaptation within Pafos forest (which is distinguished in three altitudinal subpopulations and two aspect subpopulations).The altitudinal subpopulations within Pafos forest were defined based on the Rivas-Martínez bioclimatic classification of Cyprus' bioclimatic zones [31], where: (i) 0-400 m a.s.l.corresponding to the Hot Arid to Mild Arid bioclimatic belt, (ii) 400-800 m a.s.l.corresponding to the semi wet mild bioclimatic belt and (iii) 800-1300 m a.s.l.corresponding to the semi wet cool to cool wet bioclimatic belt.} Distribution area corresponds to the area that is covered by sampled "subpopulation" and estimated either in ha in case of extended forest area or in km in case of specific altitude zone.* For the Aka and PaZ.800 the 150 megagametophytes resulted by the collection of six (6) seeds per tree from 25 different trees.For the rest "subpopulations" the sampling of seeds for genetic analysis was done as described in Section 2.3.Abbr.: Abbreviation.
PaZ.800 the 150 megagametophytes resulted by the collection of six ( 6) seeds per tree from 25 different trees.For the rest "subpopulations" the sampling of seeds for genetic analysis was done as described in Section 2.3.Abbr.: Abbreviation.

Figure 1.
The distribution of sampled subpopulations of Pinus brutia in Cyprus.

Assessment of Genetic Diversity and Structure
For assessing the patterns of genetic diversity at both intra-and intersampling levels, three windpollinated cones were collected from the middle of the trees' crown (Table 1), based on the assertion that at the middle-to-high range of a tree's crown the possibility of autogamy (self-fertilization) is much lower (practically absent) than at the low level of the crown [32].In continuation all sampled seeds from each "subpopulation" were bulked, and, hence, bulked seed material was obtained for further genetic analyses (Table 1).The genetic diversity in P. brutia from Cyprus was assessed based on haploid megagametophytes from germinated seeds assayed by isoenzyme horizontal starch gel electrophoresis.Isoenzyme analysis provides a nonrandom sampling of expressed genomic sequences and has proven invaluable in population genetic analysis over many decades of implementation.The protocols of Conkle et al. [33] and Cheliak and Pitel [34] were used to study the following enzyme systems: aspartate aminotransferase (AAT; EC 2.6.In order to assess Mendelian inheritance in P. brutia in Cyprus, six megagametophytes from each sampled tree of the Akamas (Aka) subpopulation and the altitudinal subpopulation PaZ.800 were used to derive the genotype of the (maternal) tree.Furthermore, the observed heterozygosity (Ho) and the expected heterozygosity (He) for these subpopulations were also estimated, based on the formulas by Nei [35,36].In addition, the Hardy-Weinberg equilibrium (HWE) was calculated for these two subpopulations, as well as their genetic heterozygosity.The significance level of HWE was estimated based on the differentiation between observed and expected frequency of genotypes (as this was observed in each of the two subpopulations) using the chi-squared (X 2 ) statistic test.These

Assessment of Genetic Diversity and Structure
For assessing the patterns of genetic diversity at both intra-and intersampling levels, three wind-pollinated cones were collected from the middle of the trees' crown (Table 1), based on the assertion that at the middle-to-high range of a tree's crown the possibility of autogamy (self-fertilization) is much lower (practically absent) than at the low level of the crown [32].In continuation all sampled seeds from each "subpopulation" were bulked, and, hence, bulked seed material was obtained for further genetic analyses (Table 1).The genetic diversity in P. brutia from Cyprus was assessed based on haploid megagametophytes from germinated seeds assayed by isoenzyme horizontal starch gel electrophoresis.Isoenzyme analysis provides a nonrandom sampling of expressed genomic sequences and has proven invaluable in population genetic analysis over many decades of implementation.The protocols of Conkle et al. [33] and Cheliak and Pitel [34] were used to study the following enzyme systems: aspartate aminotransferase (AAT; EC 2.6.In order to assess Mendelian inheritance in P. brutia in Cyprus, six megagametophytes from each sampled tree of the Akamas (Aka) subpopulation and the altitudinal subpopulation PaZ.800 were used to derive the genotype of the (maternal) tree.Furthermore, the observed heterozygosity (Ho) and the expected heterozygosity (He) for these subpopulations were also estimated, based on the formulas by Nei [35,36].In addition, the Hardy-Weinberg equilibrium (HWE) was calculated for these two subpopulations, as well as their genetic heterozygosity.The significance level of HWE was estimated based on the differentiation between observed and expected frequency of genotypes (as this was observed in each of the two subpopulations) using the chi-squared (X 2 ) statistic test.These subpopulations were chosen since both showed specific ecological and demographic features Forests 2018, 9, 514 5 of 27 according to the national forest inventories of the Department of Forests (1981;1991;2001;2011): The Aka subpopulation represents the driest and most degrading pine forest in Cyprus, while the subpopulation PaZ.800 corresponds to the altitude zone which is classified as the best ecological niche of the species growing in Cyprus.
For all sampling levels (subpopulations, altitude subpopulations, and orientation subpopulations) the multilocus intra-level genetic variation was assessed by the: percentage of polymorphic loci (PPL), observed number of allelic (Na), effective number of allelic (Ne), Shannon's index (I), and genetic diversity (H E ).Regarding the Pafos forest, the above measures were calculated as the mean values from the different subpopulations (i.e., altitude subpopulations and orientation subpopulations).The software GenAlEx 6.5 [37] was used for the calculation of the above interpopulation measures.
The intersampling level genetic diversity was assessed at three levels: (i) all sampled subpopulations, (ii) range-wide subpopulations (including all sampled subpopulations but without the altitude subpopulations), and (iii) altitudinal subpopulations.Furthermore, the subpopulations genetic structure was investigated using a Bayesian model-based clustering analysis [38], as implemented in the Structure v 2.3.4 software [38].Bayesian analysis was performed using the admixture and the frequency-independent allele models with 50,000 Markov chain Monte Carlo (MCMC) steps and 10,000 burn-in periods.The number of K was set (i) from 1 to 15 when all subpopulations were included in the analysis, (ii) from 1 to 9 when the range-wide subpopulations were included, and (iii) from 1 to 5 when the three altitudinal subpopulations were included; for all runs each value of K for each case was run by three replicates.Post-processing of Structure software's results, for selecting the optimum number of clusters (K) based on the Evanno method [39] and producing the graphical output, was implemented using the software CLUMPAK [40].
The genetic structure was also examined using a hierarchical analysis of molecular variance (AMOVA) as this is applied in software GenAlEx 6.5 [37], while its significance level was computed using 999 permutations.As in the previous investigation, the AMOVA was performed at three different levels (all sampled subpopulations, range-wide subpopulations and altitudinal subpopulations).At the range-wide subpopulations level AMOVA was implemented by subdividing the sampling locations into six "Groups" (i.e., each group corresponding to each sampling forest), while for the range-wide subpopulation the analysis was performed into two hierarchical levels.In addition, the hierarchical levels of genetic structure were also investigated at the marginal subpopulations, where two separate AMOVA runs were carried out: one for the three altitude gradient subpopulations, and one for the two orientation subpopulations.
The genetic relationships among sampled subpopulations (all sampled subpopulations, range-wide subpopulations, and altitude gradient subpopulations), were analyzed by means of a cluster analysis based on the unweighted pair group method with arithmetic mean (UPGMA) and the genetic distance of Nei [41], using the software TFPGA (version 1.3) [42].Bootstrap values for the dendrogram were generated using the same software, after 10,000 replications over individuals.Visualization of the genetic structure, at multivariate space was carried out based on a Principal Coordinate Analysis (PCoA) using the GenAlEx 6.5 [37].

Assessment of Morphoanatomical Diversity and Structure
The morphoanatomical variation of P. brutia in Cyprus was investigated using needles and cones.Three branches and three cones were collected from each sampled tree.For cones the following morphological traits were measured: cone length (CLen), cone width (CWid), and the ratio of cone length/width (CLen/CWid).From each of the measured cones three seeds (from the middle-part of the cone) were selected for measuring the morphological traits: seed length (SLen), seed width (SWid), length of seed's wing (SWing), and the total length of seed and wing (SLenWing).Needles were collected from the north side and middle part of tree crowns and three needles (two-years-old) from each sampled tree were used to measure 13 morphoanatomical traits: length of needle sheath (NShLen), needle length (NLen), needle width (NWid), needle thickness (NThic), number of resin ducts (internal side-NResIn), number of resin ducts (dorsal side-NResDo), total number of resin ducts per needle (NResTot), number of stomata rows (dorsal side-NStomDo), number of stomata rows (internal side-NStomIn), total number of stomata rows per needle (NStoRow), number of stomata per row (NSto/Row), total number of stomata per 1 cm 2 of needle (NStom), and number of needle teeth (NTeh).For 10 out of the 13 of the needle morphoanatomical traits, that is, apart from NLenSh, NLen, and NWin, measurements were carried out at the middle of the needle length, while the anatomic traits were measured using a stereoscope (magnification: 2 × 40).
The software SPSS 20.0 (SPSS Inc. ® ) was used to assess morphoanatomical trait variation from all sampled subpopulations at the intra-and intersampling location levels.For each trait and in each sampled subpopulation the following descriptive statistics were calculated: mean (µ), standard deviation (SD), and coefficient of variation (CV).The Spearman test was used to assess the correlation between the morphological and anatomical traits and to evaluate the correlation between traits and altitude within the Pafos forest.Furthermore, to assess morphoanatomical variation at the multivariate space level a principal component analysis (PCA) was used.The new independent components that were formed with eigenvalues above unity (>1) were used for the estimation of Euclidean distances among the sampled subpopulations.In order to visualize the classification patterns, in each of the three levels of sampling (all sampled subpopulations, range-wide sampled subpopulations, and altitudinal subpopulation) an unweighted pair group method with arithmetic average (UPGMA) dendrograms was constructed based on the morphological Euclidean distances (morphoanatomical distance), using the software NTSYS-pc 2.0 [43].In addition, the hypothesis that trees belonging to their original sampled subpopulation are morphologically and anatomically similar was tested, using back-grouping, a nonparametric classification method analogous to discriminant analysis [44,45].Finally, a Mantel test [46,47] was performed in order to investigate the possible relationship between morphoanatomical distance and genetic distance using also NTSYS-pc 2.0 [43].

Genetic Diversity and Subpopulations Structure of P. brutia in Cyprus
Eight enzyme systems encoded by 10 loci (enzyme systems AAT and MDH encoded two loci each: AAT-1, AAT-2, MDH-1, and MDH-2) were analyzed.Nine out of 10 loci were found to be polymorphic.Six loci (AAT-1, AAT-2, LAP, GDH, MDH-1, and IDH) presented two alleles, two loci (PGI and MNR) displayed three alleles, and one locus (6PGD) showed four alleles.Remarkably, one out of the three allelic detected in locus MNR was found in subpopulation PaN and in the altitudinal subpopulation PaZ.800, and one out of the two alleles of locus IDH was found in subpopulation PaN and in the altitudinal subpopulation PaZ.1200.In the locus-by-locus analysis the genetic diversity (H E ) ranged from 0 (for IDH in PaZ.400 & PaZ.800 & the monomorphic MDH-2) to 0.66 (for 6PGD in PaZ.800); the latter exhibited the highest mean genetic diversity (H E = 0.631) (Table S1).However, the highest average number of allelic (Na = 2.667) was detected in locus PGI (Table S1).
Mendelian inheritance was verified for seven out of the nine polymorphic loci (see Table S2 for more details), where X 2 and p > 0.05 at a 99% CI were tested.The outcomes from this analysis support the absence of segregation distortion of the tested loci (except for LAP where segregation distortion was found) (see Table S2).The Hardy-Weinberg equilibrium (HWE) was also calculated for the two subpopulations (Akamas and PaZ.800) where the Mendelian inheritance was investigated.Both subpopulations were under HWE for most of the study loci, since loci MNR-1 and LAP-1 showed nonsignificant HWE for the Aka subpopulation and the locus LAP-1 was not significant in PaZ.800.Based on individuals' genotype, the overall genetic diversity of Aka subpopulation was estimated: observed heterozygosity (H o ) = 0.164, expected heterozygosity (H e ) = 0.216, and the value of inbreeding Forests 2018, 9, 514 7 of 27 (F is ) = 0.241.The level of genetic diversity for PaZ.800 was slightly higher (H o = 0.216) and the H e = 0.242, while the F is was almost half (F is = 0.107).
Assessment of multilocus genetic diversity within each of the sampled subpopulations (Table 2) showed that PaN and PaZ.1200 have the highest PPL value (90%), whereas Mach, Lim, and Tro the lowest PPL value (60%).PaN was the subpopulation with the highest number of allelic per locus (Na = 2.300) and Lim the one with the lowest (Na = 1.700).The effective number of alleles (Ne) was the lowest in Lim (Ne = 1.300), but relatively similar for two out of the ten subpopulations (PaN, Ne  Furthermore, when only the samples of altitudinal subpopulations were considered, PaZ.800 presented the highest genetic diversity and PaZ.1200 the lowest, although the latter presented the highest presentence of polymorphic loci and a relatively high number of alleles per locus (Table 2).
The Bayesian clustering analyses were performed at three levels, nevertheless results were not explicitly clear.Analysis of all subpopulations, showed the highest Evanno's ∆K index for K = 2 (with relatively high statistical support; ∆K = 55.64), while two more K values demonstrated a trend of grouping, but with low statistical support of ∆K index, K = 5 (∆K = 13.01) a number corresponding to the geographic origins of sampled subpopulations and K = 10 (∆K = 10.05) a number equal to the sampled subpopulations (Figure 2a(i)).Similar outcomes were recorded when the Bayesian analysis was performed including the range-wide subpopulations, by recording the highest Evanno's ∆K index for K = 2 (with high statistical support; ∆K = 81.81)and a further peak on the graphical illustration of ∆K in K = 7 (∆K = 30.62),a number equal to the range-wide subpopulations unit (all sampled subpopulation without the PaZ.400,PaZ.800, and PaZ1200) (Figure 2a(ii)).Bayesian clustering analysis was implemented in the altitudinal subpopulations, revealing three distinct clusters (K = 3; ∆K = 56.25)(Figure 2a(iii)).The AMOVA showed that most of the genetic variation occurs within the sampled subpopulations (Table 3a,b) and the hierarchical analysis based on the subpopulations origin (grouping based on forest origin) revealed significant genetic differentiation (PhiRT = 0.118 ***) (Table 3a).The same analysis did not detect significant genetic differentiation between groups that consist of the core (central) area subpopulations and the peripheral subpopulations (forests).Quantification of genetic differentiation among all sampled subpopulations using AMOVA found a significant PhiST = 0.129 *** (Table 3a), while low but significant genetic differentiation was also detected among the altitudinal populations (PhiST = 0.018 **, Table 3b) and between the two subpopulations with different orientation (PhiST = 0.012 **, Table 3b).The AMOVA showed that most of the genetic variation occurs within the sampled subpopulations (Table 3a,b) and the hierarchical analysis based on the subpopulations origin (grouping based on forest origin) revealed significant genetic differentiation (Phi RT = 0.118 ***) (Table 3a).The same analysis did not detect significant genetic differentiation between groups that consist of the core (central) area subpopulations and the peripheral subpopulations (forests).Quantification of genetic differentiation among all sampled subpopulations using AMOVA found a significant Phi ST = 0.129 *** (Table 3a), while low but significant genetic differentiation was also detected among the altitudinal populations (Phi ST = 0.018 **, Table 3b) and between the two subpopulations with different orientation (Phi ST = 0.012 **, Table 3b).100% § Sampled subpopulations grouped in six "Groups" (based on their geographical origin) for AMOVA -Group #1: Mach; Group #2: Lim; Group #3: Aka; Group #4: Ade; Group #5: Tro; Group #6: PaN, PaS, PaZ.400, PaZ.800, PaZ.1200.† The sampled location grouped in three "Groups" for AMOVA.Group #1: Mach; Group #2: Lim; Group #3: Aka; Group #4: Ade, Tro, PaN, and PaS.} Phi RT : proportion of genetic differentiation due to differences between groups; Phi PR : proportion of genetic differentiation due to different populations within groups; Phi PT : proportion of genetic differentiation among populations among groups.d.f., degrees of Freedom.Significant level of genetic differentiation: n.s., not significant; ***, p < 0.001; **, p < 0.01; *, p < 0.05.
Pairwise genetic distances (Nei's minimum genetic distance-Table S3) among sampled subpopulations were depicted using UPGMA (Figure 3).The genetic similarities among all subpopulations reflect significant subdivisions among two major groups (Figure 3a).One group includes subpopulations Lim, PaZ.400, and PaZ.1200 and the second group the rest.However, the clades in the latter group were shown to be unimportant due to the low values of bootstraps.A similar observation was made in the range-wide subpopulations UPGMA.In this case, Lim formed a separate group (Figure 3b).Contrary to the previous dendrograms, the UPMGA on altitudinal subpopulations reflects significant subdivisions among them, since PaZ.800 seems clearly subdivided from the other two subpopulations (PaZ.400 and PaZ.1200) with the highest bootstrap value (Figure 3c).A PCoA was used to discover and depict the major patterns within a multivariate dataset, by detecting the relationship between the distance matrix elements in a two-dimensional space.When all sampled subpopulations were considered, PCoA revealed four loosely formed groups that do not correspond well to subpopulation geographic origin (Figure 4a).On the other hand, when analysis was performed at the range-wide subpopulations, Ade and Tro were completely isolated from the remaining populations, while the orientation subpopulations (PaN and PaS) grouped with the Mach population (Figure 4b).Contrary to the above analyses that considered the latitudinal and longitudinal sampling of populations, the visualization of the genetic structure of the altitudinal populations showed a clear disjunction among the populations of the three zones sampled (Figure 4c).
Given the significant influence of altitude in subpopulations genetic structure, the Spearman's correlation coefficient analysis examined the relation between allele frequency across loci and altitudinal subpopulations.It revealed a positive significant correlation between the altitudinal subpopulations in a single case, namely regarding one allelic in locus 6PGD.A PCoA was used to discover and depict the major patterns within a multivariate dataset, by detecting the relationship between the distance matrix elements in a two-dimensional space.When all sampled subpopulations were considered, PCoA revealed four loosely formed groups that do not correspond well to subpopulation geographic origin (Figure 4a).On the other hand, when analysis was performed at the range-wide subpopulations, Ade and Tro were completely isolated from the remaining populations, while the orientation subpopulations (PaN and PaS) grouped with the Mach population (Figure 4b).Contrary to the above analyses that considered the latitudinal and longitudinal sampling of populations, the visualization of the genetic structure of the altitudinal populations showed a clear disjunction among the populations of the three zones sampled (Figure 4c).
Given the significant influence of altitude in subpopulations genetic structure, the Spearman's correlation coefficient analysis examined the relation between allele frequency across loci and altitudinal subpopulations.It revealed a positive significant correlation between the altitudinal subpopulations in a single case, namely regarding one allelic in locus 6PGD.The Spearman correlation among the investigated morphological and anatomical traits showed that 52 out of the 190 paired correlations (matrix table of 20 traits) are statistically significant for p-value > 95% (Table S4).Interestingly, significant correlations between needle size and the stomata rows, between cones traits as well as between the cone and seed traits, were found.Investigation on the association (Spearman correlation) between morphoanatomical traits and altitude, detected positive correlations for 12 out of the 20 traits; namely 10 correlations were positively significant and two negatively significant (Table 5).Despite the fact that the correlation coefficient is relatively low, the analysis showed a significant increase of the size of specific traits (e.g., morphological traits for cone and seed; morphoanatomical traits: NShLen, NResTot, and NStom) as the altitudinal subpopulations increased from the PaZ.400 to PaZ.1200.The identification of the traits that contribute more significantly in the overall phenotypic variation observed was investigated by PCA.The use of the eingenvalues (e.g., Kaiser's criterion), reduced the dimension of the 20 morphoanatomical traits to nine axes (for components see Table 6), of which the first six explain 97% of the total variance (Table 6).The first axis presents strong correlations with the initial variable expression of NWid, NThic, and SLen, explaining 40.60% of the total morphological variance.The second axis, accounting for 27.90% of the total variation, was associated with traits NLen and NStom, while the third axis interpreted 14% of the total variation, with NShLen and SWid to be the associated traits for this axis.The next three axes accounted for 14.6% of the overall variation, and were associated with the anatomical traits (stomata and resin ducts) and with the morphological traits of cones and seeds (Table 6).Furthermore, a morphoanatomical Euclidean distance matrix (Table S3) among the sampled subpopulations was produced using the PCA first six axes.The subpopulation analysis showed that the highest morphoanatomical distance was recorded between Aka and PaZ.800 (3.188) and the lowest between PaS and PaZ.800 (0.292).On the other hand, based on the altitudinal subpopulations, the morphoanatomical distance between PaZ.400 and PaZ.1200 was the highest (1.482) distance occurring compared to the other morphoanatomical distances recorded between the PaZ.800 (middle range) and the other altitudinal subpopulations (Table S2).The illustration of these distances in a UPGMA dendrogram indicated a specific grouping pattern (Figure 5).The subpopulations from Pafos forest (PaN, PaS, PaZ.400, PaZ.800, and PaZ.1200) and Tro shaped a clear geographically defined group and formed a common clade.The other clade was built by the Mach and Ade subpopulations.The geographically isolated subpopulations of Lim and Aka were incorporated in two different (separate) clades on the dendrogram.With regards to the altitude gradient subpopulations analysis, and contrary to the UPGMA based on genetic distances, the morphoanatomical Euclidean distances were lower between PaZ.400 and PaZ.800, whereas the PaZ.1200 was grouped in a single clade (Figure 5).Asia Minor where H E = 0.265 [53] using isoenzyme loci.The higher number of rare alleles recorded in the Cyprus population of P. brutia, for specific loci (i.e., loci: AAT, PGI, LAP, 6PGD, IDH, and MNR), compared to the populations from its continuous range [52,54,55], implies that the P. brutia peripheral population in Cyprus contains unique genetic variants due to its distinctive evolutionary history.Thus, this population can be appreciated as a living gene bank for forest genetic resources for this species.The nonsignificant inbreeding coefficient values observed in the two populations (Aka and PaZ.800), in conjunction to the presence of Hardy-Weinberg equilibrium, support the general argument, as in most conifers, of a random mating system.These results indicate that these two subpopulations (which showed reverse ecological and demographic features) are not characterized by a founder effect and have not been affected by strong bottleneck events.Hence, the P. brutia forest in Cyprus may have been somewhat preserved from genetic bottlenecks; especially during displacement and confinement of populations in Pleistocene glaciations, or by raised water levels during intervening warm periods as this event has been observed for several species in Europe [12,56].

The Patterns of Morphoanatomical Diversity and Structure in P. brutia
The above arguments are also supported by the comparison of several morphoanatomical traits from the present study and other studies on P. brutia populations from other islands (i.e., Rodos, Crete, etc.) and/or populations of larger origins (continuous populations).The coefficient of variation (CV) [57,58] and the mean values [57][58][59] for numerous morphoanatomical traits show similarities between the populations of Cyprus and those from other origins, which implies that the morphological and anatomic traits of the former are not affected by restriction of genetic diversity and genetic drift (or strong inbreeding events).A remarkable observation on Cypriot P. brutia is the significant correlation between morphological and anatomical traits (Table S4).This morphoanatomic association in this rear-edge population could be a consequence of macro-environmental differentiation from the continuous distribution range of the species.Such adaptive needle traits have been mentioned for other coniferous species [60].Often the acclimation to their environment leads to the development of specific adaptive phenotypic responses by altering their length, width, number of stomata, their angle towards the shoot, or by forming needle clumping [61][62][63][64].In the current study, the correlation between NLen and Nstom, indicates that long needles are associated with a high number of stomata, an observation which could be a relative advantage for Cypriot P. brutia, providing effective resilience and adaptability for this peripheral population in different environmental conditions.

Patterns of Genetic Diversity of P. brutia in Cyprus
A nonuniform genetic diversity distribution across longitudinal, latitudinal, and altitudinal distributions was detected.Genetic variation among subpopulations is most likely a consequence of different demographic and evolutionary events.The intensive and negative impact of human activities on forests resilience and the deforestations in Cyprus, since the first human presence (11 th millennium B.C.), is well-known [65,66].During the Ottoman period (1570-1878 A.C.) reported goat stocking rates [67] are far above the forest carrying capacity, while the Cypriot forests were repeatedly logged in order to cover the energy needs of Bronze Age copper production (c.3300-1200 BC) [68].These historical facts may reasonably be linked to the fragmentation of P. brutia forest and lead to negative pressures on genetic variability within specific relic subpopulations, and consequently to different genetic variation patterns.Thus, the divergence of genetic variation (H E ) and effective number of alleles (N e ) between subpopulations, seems to be linked to high past anthropogenic pressure, particularly for subpopulations located in peripheral and more isolated forests (Aka, Lim, and Mach) where the pressure was higher (see Figure 1 and Table 2).
The clustering of subpopulations in two groups (Structure analysis; K = 2) and the detection of significant genetic differentiation among the subpopulations (AMOVA, see Table 3), alongside with the unclear geographic clustering of subpopulations from UPMGA and PCoA, could be attributed either to the fragmentation of an earlier larger and uniform population, or to the fact that the present forests are relics of previously differentiated populations of P. brutia.This question could not be directly answered by the present study, as more powerful molecular markers may be needed, for a more clear estimation of the best K-clustering of P. brutia in Cyprus.Meanwhile, the existence of landscape barriers to effective gene flow among the current subpopulations, shape genetic differentiation through generations.
Comparison of the genetic differentiation values (Phi PT = G ST = 0.129) between this and other studies reveals that a notably higher level of differentiation was observed in the case of P. brutia in Cyprus.In particular, populations of P. brutia originating from islands of the north-eastern Aegean sea showed G ST = 0.021 [55], whereas populations from south Asia Minor (G ST = 0.053) support the previous argument that rear-edge populations shape disproportionately high levels of genetic differentiation present even between geographically proximal ones, leading to exceptionally high levels of regional genetic diversity, than the populations of the continuous range [12,[69][70][71].
Fingerprinting of the genetic diversity in Pafos forest allowed assessing the impact of local landscape on genetic patterns.Despite the fact that the comparison of northeast and southwest (aspect) subpopulations detected equal expected heterozygosity (H E = 0.228 and H E = 0.220, respectively), a low but significant genetic differentiation was detected (Phi PT = 0.012).Moreover, the three altitude gradient subpopulations showed significant but low differentiation (Phi PT = 0.018) and structure (K = 3); while the middle altitude gradient subpopulation (PaZ.800)recorded the highest genetic diversity (H E = 0.244) (Table 2).The highest value of H E was found in PaZ.800, being also the highest value among all subpopulations in Cyprus.This result is in concordance with other studies regarding P. brutia, particularly in the Taurus mountains, where the middle altitude zone recorded a higher genetic diversity than other zones [53].The significant differentiation among the three altitude gradient subpopulations could be attributed to a combination of anthropogenic activities and small scale disturbance, or to the processing of different genetic evolutionary factors within each subpopulation (at local microscale), after their last postglaciation separation.An alternative explanation, could be altitudinal movements amplified by local topography (upward and downward movement within a single mountain region) during Pleistocene glaciations and interglaciations.Such an option has been presented for Euro-Mediterranean ecosystems [70], and is supported for instance by the findings on Cedrus brevifolia, a narrow endemic tree in Cyprus [72].Therefore, this shift at different altitudinal gradient zones (i.e., ecological niche) during interglaciation and postglaciation periods is potentially the reason for the formation of an admixture zone.Alternatively, the significant differentiation among the altitude gradient subpopulations could be attributed to the processing of different evolutionary factors at microscale, after their last postglaciation partition.
The hypothesis of local (microenvironmental) adaptation dynamics is supported by relevant literature which demonstrates several examples of wild populations on ecologically marginal sites (i.e., altitudinal gradients, and different ecological aspects) [73][74][75][76] and further reinforced by the presence of altitudinal clinal variation in P. brutia in the various morphoanatomical traits, since significant positive correlation between traits and altitudinal variation was detected (discussed below).

Patterns of Morphoanatomical Traits of P. brutia in Cyprus
The sampled subpopulations exhibit varying degrees of diversity in the 20 morphoanatomical traits studied, in relation to the longitudinal or the altitudinal gradient.In particular, morphological and anatomical traits constituted powerful tools for describing the phenotypic diversity and structure in the present study.The interpretation of the morphoanatomical trait diversity patterns, suggests that the observed diversity and structure, possibly are the result of the manifestation of phenotypic plasticity at different micro-environments.High phenotypic plasticity, leads to the species being distributed in a wider geographical and ecological range.The morphoanatomical traits detected a similar range of values between the populations of Cyprus and those from other origins, which bring the notion of the connection of phenotypic plasticity to genetics (see a previous paper [77]) and to the underlying quantitative trait loci variation that shapes phenotypic patterns.Phenotypic plasticity is reflected in: (i) the clustering of sampled subpopulations in a morphoanatomical dendrogram (Figure 5) where clustering of the subpopulations is concordant to their geographical origin, a relation not seen in the genetic data; and (ii) the back-grouping method results (Figure 6), in which the isolated subpopulations illustrated the highest value of back-grouping (Aka, Lim, and PaZ.1200).Phenotypic plasticity among subpopulations could be a consequence of different reaction norms, where a set of phenotypes can be produced by an individual genotype when exposed to different environmental conditions, such as different soil type and meteorological conditions.Furthermore, phenotypic plasticity is manifested in the positive significant clinal variation of the morphoanatomical traits in the altitudinal range subpopulations (PaZ.400,PaZ.800, PaZ.1200).This clinal variation across the altitudinal range may be due to trade-offs between plasticity and stress tolerance in harsh environments, as indicated for other species [78][79][80].Thus, P. brutia from Cyprus showed an increased number of resin ducts and stomata per needle, as well as the largest size of cones and seeds as the altitude increased within the mountain elevation gradient.
The morphoanatomical trait patterns observed could be the consequence of isolation and phenotypic plasticity seen over different environmental gradients.A similar interpretation was also invoked in analyses of other peripheral populations of the same species in Crete [58], in other tree species, such as Abies cephalonica in Mt.Parnitha [81] and Fagus sp. in western Eurasia [82].Therefore, this study further supports the notion that the evolution of plasticity increases the response to selection, thus reducing maladaptation induced by gene flow (see previous works [83,84]).The same studies support that, in interaction with local population growth, the evolving plasticity allows a species to occupy a larger geographical range, presenting higher plasticity in marginal than in central habitats [84].
The morphological and anatomical diversity of P. brutia in Cyprus is greatly dependent on a number of needle (NWid, NThic, NLen, NStomRow, and NShLen) and seed (SLen and Swid) traits, with high loadings in the first three axes of PCA where 82.5% of the total variance measured is explained (Table 6).In addition, the significant correlation between specific traits, such as needle size (i.e., between length and width), needle morphology, and anatomy (i.e., between length or width of needle and resin ducts or stomata), and cone size with seed size (i.e., between cone width and total length of seed and wing), probably implies that the species developed mechanisms and characters that permit different morphological and anatomical harvesting strategies.Peripheral subpopulations (i.e., Lim, Aka, and PaS) appear to develop different mechanisms to survive in different micro-environments.They are characterized by short, wide, thick needles (especially Aka), which function as an adaptation mechanism (morphological features) in dry sites [85].The Aka and the PaZ.1200 (the subpopulation at the highest altitude) present the highest number of resin ducts.This feature is linked to a plant species mechanism under extreme environment conditions (dry or cold), since resin ducts seem to protect vascular tissues and ground tissues in most species [85,86].Alternatively stomata numbers are low at the driest micro-environments (subpopulations: Aka, Lim, and Tro), which implies a development of adaptation mechanisms against water loss in the summer.Similar adaptation mechanisms to dry environments were recorded in other studies [81,85].Cone and seed size was smaller in peripheral, drier, and ecologically degraded subpopulations (Aka, Lim), while the opposite results were observed in subpopulations of higher altitude, where the same traits present clinal variation.These patterns were also observed in other studies of P. brutia and have been attributed to ecological factors, in particular to the brief sprouting period which is coupled to low temperatures [58,85].Such factors are apparently present in this study, particularly at a higher altitude.

Inferences from the Study of a Mediterranean Oceanic Island's Peripheral Tree Population
The population of P. brutia in Cyprus is an example of an isolated (both geographically and ecologically) island peripheral population, whose patterns of genetic diversity were shaped by past demographic and ecological stochasticity.The detection of high genetic diversity of P. brutia in Cyprus may respond to more stable effective population sizes in eastern Mediterranean forest tree populations as argued for Mediterranean conifers by Fady [8].This study likely constitutes an example of a peripheral forest tree population that does not tend towards genetic erosion, neither towards increased inbreeding or genetic drift, when compared to core populations.These findings are contrary to earlier observations that peripheral tree populations showed genetic erosion and elevated inbreeding leading to genetic drift, in comparison to core populations [87][88][89].The detection of significant genetic structure (Bayesian and AMOVA) among geographically sampled subpopulations shows that the population of P. brutia in Cyprus is not genetically homogeneous, comprising of genetically different subpopulations.Thus, ongoing genetic evolution dynamic processes seem to occur within the study population, in spite of its small geographical distribution within an island.
In the present study, the peripheral population of P. brutia presented a similar value range regarding morphoanatomical traits, to populations from other origins in the Mediterranean Basin (island and mainland).This is in contrast to earlier findings suggesting that in peripheral populations' phenotypic trait descriptive statistics, show different values from those found in core populations, indeed for a limited range of traits usually related to growth [90].The phenotypic plasticity that was observed in the present study is of particular importance for the adaptive potential of the targeted population in an island environment.Such phenotypic plasticity, especially for long-lived plants, is particularly important since it allows wild organisms to accommodate rapid change in local environmental conditions [91,92].

Figure 1 .
Figure 1.The distribution of sampled subpopulations of Pinus brutia in Cyprus.

ForestsFigure 2 .
Figure 2. Bayesian cluster analysis of the optimum K clusters.(a) The optimal number of clusters (K) based on the Evanno method: (i) all sampled subpopulations, (ii) range-wide subpopulations, and (iii) altitudinal subpopulations.(b) Bar plot with individual mosquitoes represented as vertical bars colored in proportion to their assignment to clusters inferred at (i) all sampled subpopulations-K = 2, (ii) range-wide subpopulations-K = 2, and (iii) altitudinal subpopulation-K = 3.

Figure 2 .
Figure 2. Bayesian cluster analysis of the optimum K clusters.(a) The optimal number of clusters (K) based on the Evanno method: (i) all sampled subpopulations, (ii) range-wide subpopulations, and (iii) altitudinal subpopulations.(b) Bar plot with individual mosquitoes represented as vertical bars colored in proportion to their assignment to clusters inferred at (i) all sampled subpopulations-K = 2, (ii) range-wide subpopulations-K = 2, and (iii) altitudinal subpopulation-K = 3.

Figure 3 .
Figure 3.The genetic similarities in Pinus brutia illustrated by UPGMA dendrogram based on Nei's minimum genetic distance.(a) Sampling location level; (b) subpopulation level; (c) altitude level.

Figure 4 .
Figure 4. Two dimensional principal component analysis (PCA) plot of the genetic distance with regard to the first two principal components.(a) All sampled subpopulations; (b) wide-range subpopulation level; (c) altitudinal subpopulations.

Table 4
presents the statistical description of morphoanatomic traits of P. brutia in Cyprus.Aka and Lim presented the lowest values of morphological traits for cones and seeds, while cone and seed length and width increased from low to high elevation in the altitudinal analysis.In addition, needle morphoanatomical traits showed that needle length (NLen) is variable among sampled subpopulations, since Lim presented the lowest values (106 mm), while Ade and PaN presented the largest (131 mm and 128 mm, respectively).Aka showed the highest values of NWid (1.43 mm) and NThic (0.88 mm), while PaZ.800 showed the lowest values NWid (1.25 mm) and NThic (0.77 mm).Concerning needle anatomical traits, Aka presented the highest values for resin ducts (NResIn: 4.23; NResDO: 6.72; NResTot: 10.95) and PaS showed the lowest values (NResIn: 3.15; NResDO: 5.21; NResTot: 8.36).In addition, resin ducts appear to increase with the altitude zonation.Regarding the five anatomical traits of stomata, Aka, Lim, and Tro, presented lower values (measurements relative to stomata traits) than the other subpopulations.

Figure 4 .
Figure 4. Two dimensional principal component analysis (PCA) plot of the genetic distance with regard to the first two principal components.(a) All sampled subpopulations; (b) wide-range subpopulation level; (c) altitudinal subpopulations.

3. 2 .
The Patterns of Morphoanatomical Diversity and Structure in P. brutia Table 4 presents the statistical description of morphoanatomic traits of P. brutia in Cyprus.Aka and Lim presented the lowest values of morphological traits for cones and seeds, while cone and seed length and width increased from low to high elevation in the altitudinal analysis.In addition, needle morphoanatomical traits showed that needle length (NLen) is variable among sampled subpopulations, since Lim presented the lowest values (106 mm), while Ade and PaN presented the largest (131 mm and 128 mm, respectively).Aka showed the highest values of NWid (1.43 mm) and NThic (0.88 mm), while PaZ.800 showed the lowest values NWid (1.25 mm) and NThic (0.77 mm).Concerning needle anatomical traits, Aka presented the highest values for resin ducts (NResIn: 4.23; NResDO: 6.72; NResTot: 10.95) and PaS showed the lowest values (NResIn: 3.15; NResDO: 5.21; NResTot: 8.36).In addition, resin ducts appear to increase with the altitude zonation.Regarding the five anatomical traits of stomata, Aka, Lim, and Tro, presented lower values (measurements relative to stomata traits) than the other subpopulations.

Table 1 .
Geographic location of samples along the Troodos mountain range and sample size per sampled subpopulation.

Table 2 .
Patterns of genetic variation at multilocus level over subpopulations in P. brutia in Cyprus.
PPL: % of polymorphic loci; Na: Number of allelic per locus; Ne: Effective number of allelic; I: Shannon's index; H E : Genetic diversity.

Table 3 .
Summary of the hierarchical analysis of molecular variance (AMOVA).(a) Hierarchical AMOVA based on the sampled subpopulation patterns, (i) all sampled subpopulations and (ii) range-wide subpopulations.(b) Hierarchical AMOVA based on marginality level, (i) altitude gradient subpopulations and (ii) different orientation subpopulations.

Table 4 .
Statistical measures (average value, standard deviation, and coefficients of variance) of morphological and anatomical traits in P. brutia from Cyprus: (a) the values in each sampled subpopulations and (b) the overall value of morphoanatomical traits in P. brutia.

Table 6 .
Principal component analysis (PCA) of the 20 morphoanatomical traits of Pinus brutia in Cyprus.
* Variables that showed strong correlation with each component.