Introduction

Phoenix theophrasti, known as Theophrastus’s Date Palm, Cretan Date Palm and Datça Date Palm, is one of the two palm species native to continental Europe. It is a perennial monocotyledon species belonging to the Arecaceae (Palmae) family. It has been known in the Mediterranean since classical times, when it was recorded by Theophrastus in Enquiry into Plants (370–285 BC) (Hort 1916). However, it was not until 1967 that the species was discovered and described formally by Greuter (1967) from a wild population in Crete-Greece. It was considered the same species (a wild or feral type) with date palm (P. dactylifera L.), before Greuter (1968) identified it as a different species (Liolios et al. 2009). Greuter named it in honor of the Greek botanist-philosopher Theophrastus. Until Boydak and Yaka (1983) and Boydak (1985) reported its occurrence in Turkey, occupying rather large areas in the Datça Peninsula, P. theophrasti was considered a relict endemic to Crete. This record was an important addition to the flora of Turkey, representing the natural occurrence of this palm species in Asia, too. In Turkey, it is the only native palm species. In the Datça Peninsula, two major populations of P. theophrasti exist on the northern and southern hillsides of the mountains extending in a west–east direction, respectively. The population at the northern hillside exists along the Eksera Creek, extending from sea level to 225 m of altitude in an area of about 1 km as the crow flies. The population at the southern hillside begins from sea level as well and reaches an altitude of 350 m in an area of 7 km as the crow flies. Hurmalıbük and Dimitri Bay are the main locations of the species at the southern hillside. Another well-distinct population of P. theophrasti was located in Kumluca-Karaöz at the coast of Finike Bay by Boydak in 1985 (Boydak 1986, 1987). This population extends from sea level to 50 m. In both Datça and Kumluca-Karaöz, the species is mainly restricted to the deep, isolated coastal valleys, although some small groups of individual trees are present along the coasts, especially in Hurmalıbük location. The last big location for P. theophrasti was reported in Bodrum-Gölköy and since it has certain characteristics that differ from P. theophrasti it has been classified as P. theophrasti ssp. Gölköy (Boydak and Barrow 1995). This population is found growing on boggy ground on the boundaries of Gölköy village between patches of Pinus brutia Ten. and the sea. The Gölköy palm differs from P. theophrasti in terms of its fruiting stalk length, fruit size and seed shape. Fruit stalks of Gölköy palm are 0.6–2 m long, whereas those of P. theophrasti almost never exceed 0.3 m. Gölköy palms have marginally larger fruits compared to the fruits of P. theophrasti, and these fruits have a slightly sweet taste, reminiscent of the cultivated date P. dactylifera. The fruit of P. theophrasti is an oval yellowish-brown drupe of about 1.5 cm long and 1 cm diameter, holding a single large seed. Since its fruit pulp is too thin, too fibrous and has an acrid taste, it has no agricultural significance, although they are sometimes eaten by the local people. With 4–8 m height, the Gölköy palms are generally shorter than P. theophrasti of Datça and Kumluca-Karaöz, where some palms up to 17 m height can be seen (Boydak and Barrow 1995). The fourth location of the species is in Patara, which is a very small population consisting of only 16 trees (Vardareli 2012). Because of the harsh environmental conditions of the locations of the species, it is very difficult to estimate the total number of individuals in Turkish populations. However, according to our field experiences, we can say that about 2000 individuals would be a too optimistic guess for the total number in all populations. The existence of P. theophrasti always indicates the presence of a water source in its habitat. Salt tolerance of the species enables it to survive combined pressures of exposure to coastal winds and seawater (Barrow 1998). The presence of P. theophrasti in Turkey is a good example of the floristic linkage existing between Anatolia and Crete.

As a more recently described species of the Phoenix genus, P. theophrasti is not well characterized and not well established as a valid species (Krueger 2011). Although it had been identified as a different species by Greuter (1968), it was sometimes thought to be simply a population of distinct P. dactylifera native or naturalized in the region before being validated as a species by Barrow (1998). It is simply due to its similar appearance to P. dactylifera, although apparently somewhat smaller. It is still believed that these two neighbor species are close relatives, but a need for research on detecting whether these two species are relatives or the same species is frequently emphasized in the literature (Barrow 1998; Citrus and Date Crop Germplasm Committee 2004; Krueger 1995, 1998). Due to their high nutritive, economic and social values, varieties of P. dactylifera have been subject to many genetic diversity studies (Sedra et al. 1998; Khierallah et al. 2011; Srivashtav et al. 2013). However, to the best of our knowledge, there is no comprehensive phylogenetic study showing the relation between P. theophrasti and P. dactylifera, as well as the relations of P. theophrasti with other palm species.

Within the genus, a paucity of taxonomically useful gross morphological and anatomical characters makes the Phoenix species difficult to distinguish one from the other, as well as making the clarification of the number of the species in the genus. Moreover, the general tendency of the palms, under the Arecaceae (Palmae) family, forms intragenus and intergenus hybrid species under natural or unnatural (such as human-mediated transport) conditions making it difficult to clarify species status at both the genus and family levels. All species of the Phoenix genus are dioecious, have the identical chromosome number of 2n = 36, and freely hybridize (Röser 1994; Rivera et al. 2008).

Intense human activities have tremendously affected Earth’s ecosystems. Impacts from human activities range in a broad scale from direct management of species to population fragmentation to global climate change. Reid et al. (2005) reported that in the last 100 years, human activities have caused an estimated 50 to 1000 times more extinctions than would have occurred due to natural processes. According to the International Union for Conservation of Nature and Natural Resources (IUCN) Red List out on July 19, 2012, at Rio + 20 Earth Summit, the IUCN being a well-respected authority on the conservation status of a species, a total of 63,837 species were assessed. Of those, 19,817 were described as threatened with extinction, among which 3947 were listed under “critically endangered,” 5766 under “endangered” and more than 10,000 species under “vulnerable.” Almost half of the threatened species were plants.

The southwest of Turkey and southeast of Greece are the ancient areas of civilization and the dense human populations that have inhabited these areas have seriously degraded most of the original habitats. FAO (2010) reports that, as one of the original habitats, coastal forests hosting very unique plant communities will be especially vulnerable to climate change impacts, e.g., sea level rise and changes in the salinity conditions. This is the case of the Liquidambar orientalis and P. theophrasti forests which are the most notable tertiary relict endemic trees of the southwestern coast of Turkey (Médail and Quézel 2003), as a part of the eastern Mediterranean (Vogiatzakis and Rackham 2008; Kozlowski et al. 2014).

Considering the need to distinguish P. theophrasti from P. dactylifera, the main focus of this study was to genetically characterize P. theophrasti populations and its relationship to the other palm species. The Mediterranean and adjacent regions are considered among the most important refugial zones of the tertiary relict flora in western Eurasia (Milne and Abbott 2002; Milne 2006). Most of the relicts of the region have a restricted distribution and are at risk of extinction, with only a few natural populations that survived the last glaciation in isolated areas (Quézel and Médail 2003; Connor 2009; Garfi et al. 2011). At the same time, relict species provides a unique opportunity to understand past and recent biogeographical and evolutionary processes. Therefore, P. theophrasti is of scientific and conservation interest and this study is an important attempt in this respect. We characterized six Turkish P. theophrasti populations on molecular level by using SSR markers for the first time in the literature. We also assessed the relations of P. theophrasti with nine different palm species which are P. dactylifera, P. reclinata, P. rupicola, P. roebelenii, P. canariensis, P. loureiri, P. acaulis, P. sylvestris from Phoenix genus and Chamaerops humilis from Chamaerops genus by the same markers for the first time in the literature. Chamaerops humilis was included in the study as the second natural palm species of Europe, having a close natural distribution area with P. theophrasti, to determine the genetic relations among them.

Knowledge of the extent and structure of the genetic variation in plant populations is necessary to develop effective strategies for conservation, to understand the processes of evolution and for the selection and preservation of wild relatives of important crops (Frankel and Bennett 1970). The development of the molecular marker systems has given geneticists powerful new tools to use in the assessment of genetic diversity in plant populations. Additionally, molecular markers are increasingly used to address species delimitation problems where morphological methods are unreliable and inconclusive (Henderson et al. 2006). Simple sequence repeats (SSRs) or microsatellites are tandem repeats of 1–6 bp in eukaryotic genomes (Tautz and Renz 1984). They are co-dominant, multiallelic markers (Innan et al. 1997; Senior and Heun 1993) characterized by their great abundance (Condit and Hubbell 1991; Röder et al. 1995), high variability due to slippage during DNA replication (Schug et al. 1998) and large distribution throughout different genomes (Liu et al. 1996; Taramino and Tingey 1996; Röder et al. 1998). Also, they are simple to detect, transferable and reproducible (Zietkiewicz et al. 1994). By these properties, they are useful tools for interspecific and intraspecific studies (Henderson et al. 2006).

Materials and methods

Plant material

Five P. theophrasti populations and one P. theophrasti ssp. Gölköy population sampled and used in the study are shown in Fig. 1. A total of 86 plants were genotyped in this study (Table 1). Leaf samples from individual plants in each population were collected in spring 2011 with the transect method. Each transect consisted of a straight line. In each population, the length, numbers and directions of transects were determined depending on the size of the population and other characteristics of the location such as slope, topography and uniformity. Depending on the size of the population, one plant was sampled at every 10–100 m along the transect. Coordinates and the altitude ranges for the first and last samples taken from the locations were recorded in the field by a GPS device (Garmin GPSMAP 62 Model) and given in Table 1. The collected samples were stored at − 20 °C in a freezer until the DNA extraction procedure.

Fig. 1
figure 1

Locations of the six Phoenix theophrasti populations sampled in the study and genetic clustering of these populations assessed from structure analysis. Three clusters were inferred: Cluster 1 is red, cluster 2 is blue, and cluster 3 is green in color. The vertical coordinate of each cluster indicates the membership coefficient for each individual. Bars with 2 or 3 colors represent genotypes that show admixture character. Geographically presented proportions of assignment of each population to each of the three clusters were shown by pie charts

Table 1 Coordinate and the altitude ranges for the first and the last individuals sampled from linear transect designed across the species population for each location and the numbers of the individuals sampled from the populations

Five individuals from each nine different palm species were obtained from Köyceğiz-Palm Center by the direction of the Directorate of Forestry of Muğla Province. For the characterization of the genetic relationships between P. theophrasti and the other palm species, five individuals of P. theophrasti were randomly selected among the individuals of the six populations sampled and used in the analysis together with the individuals of the nine palm species.

Genomic DNA extraction

Genomic DNA was extracted from young leaves using a protocol optimized according to Doyle and Doyle (1987), Cullings (1992) and Sharma et al. (2003). In this optimized protocol, 2 g of leaf tissue was submerged in 5 ml of 96% ethanol for 30 min. After the evaporation of ethanol, tissue was ground with a mortar and pestle. Then, another 5 ml of 96% ethanol was applied to the ground tissue and the sample was left in the hood for evaporation of the ethanol. After grinding the tissue once more, 850 μL of pre-warmed (65 °C) CTAB buffer (0.1 M Tris, pH 8.0, 1.4 M NaCl, 0.02 M EDTA and 2.0% of cetyltrimethyl ammonium bromide) was added to 0.05 g of the fine powdered tissue, vortexed very well and incubated in 65 °C water bath for 60 min by vortexing at the 20th, 40th and 60th minutes. This was followed by addition of 850 μL chloroform/isoamylalcohol (24:1). The samples were mixed well by inverting the tubes and centrifuged for 10 min at 14,000 rpm. Next, the top aqueous layer was transferred to another tube. The volume of the aqueous phase was estimated to be ~ 300–350 μL, and 0.08 volumes of cold 7.5 M ammonium acetate and 0.54 volumes (using the combined volume of aqueous phase and added ammonium acetate) of cold isopropanol were added to precipitate the DNA. After incubating at − 20 °C for 30 min, the sample was centrifuged at 14,000 rpm for 4 min to pellet the DNA. The supernatant was removed, and 700 μL of cold 70% ethanol was added to wash the pellet, followed by incubation at − 20 °C for 15 min prior to centrifugation at 14,000 rpm for 3 min. Then, the ethanol was removed and the pellet was dried and resuspended in 120 μL of TE, pH 8.0 (10 mM Tris, 1 mM EDTA). Next, 10 μg RNase (Fermentas, Lithuania) was added and incubated at 37 °C for 60 min mixing thoroughly but gently. After incubation, 400 μL of chloroform/isoamylalcohol (24:1) was added, mixed thoroughly for 10 min and centrifuged at 14,000 rpm for another 10 min. The top aqueous layer was transferred to a clean tube, and two volumes of 96% ethanol were added, followed by precipitation at − 20 °C for 30 min, and centrifugation at 14,000 rpm for 4 min to pellet the DNA. The supernatant was removed, and 700 μL of 70% ethanol was added to the pellet, followed by incubation at − 20 °C for 15 min prior to centrifugation at 14,000 rpm for 3 min. Finally, the ethanol was removed; the pellet was dried and resuspended in 120 μL of TE, pH 8.0.

Amplification of SSR loci

Twelve variable and polymorphic SSR loci composed of dinucleotide GA repeats previously identified in P. dactylifera, whose transferability was shown across the genus Phoenix including P. theophrasti (Billotte et al. 2004), were amplified for all samples (Online Resource 1). PCR amplifications were carried out in 20 µL reaction volumes, each reaction containing 25 ng of genomic DNA, 1 × PCR buffer, 2 mM MgCl2, 0.16 mM dNTPs (from each), 0.5 U Taq DNA polymerase, 0.2 µM forward primer (F; labeled with one of the FAM, NED, HEX fluorescent dyes) and 0.2 µM reverse primer (R). The PCR program was the initial denaturation at 95 °C for 1 min, then 35 cycles of 94 °C for 30 s, 52 °C for 1 min and 72 °C for 2 min and a final elongation step at 72 °C for 8 min, for all 12 SSR loci screened. One microliter of each reaction was visualized on 2% agarose gel, and the PCR products were analyzed at ABI PRISM 3100 Genetic Analyzer (Applied Biosystems). Electropherograms were manually checked with the Applied Biosystems Peak Scanner program (http://www.appliedbiosystems.com) and recorded.

Data analysis

Genetic variability parameters (average number of alleles Na, effective number of alleles Ne, observed heterozygosity Ho and expected heterozygosity or gene diversity He) for each microsatellite locus, for each P. theophrasti population and for each palm species included in the study were calculated using POPGENE version 1.31 (Yeh et al. 1999). The same program was used to calculate the genetic distance (DN) values among P. theophrasti populations and among palm species according to Nei (1972) as well as to assess F-statistics and gene flow (Nm) values for P. theophrasti populations and for palm species. If Nm > 1, gene exchange among populations can prevent the impact of genetic drift and decrease the genetic variance among populations (Slatkin 1985). Wright (1978) suggested a qualitative guideline for the interpretation of FST in which the range 0.15–0.25 indicates “great genetic differentiation” and values of FST above 0.25 indicate “very great genetic differentiation.” Genotypic frequencies for conformance to the Hardy–Weinberg equilibrium (HWE) were tested with the Chi-square (χ2) test at a significance level of 5% and likelihood ratio G2. Corrections were performed by applying the Bonferroni correction for multiple comparisons (Rice 1989). The degrees of genetic differentiation among populations and among species were measured as pairwise FST values, and the analysis of molecular variance (AMOVA) tests for investigating the genetic structure of P. theophrasti populations and palm species grouped according to their geographical distribution were performed using the Arlequin ver. 3.5 program (Excoffier et al. 2009). AMOVA test estimates and separates the total molecular variance into components and then tests the significance of partitioned variances. In the AMOVA test of P. theophrasti populations, Hurmalı, Dimitri, Eksera and Gölköy populations were included in group 1, while Patara and Karaöz populations were included in group 2 (Fig. 1). In the AMOVA test among palm species, P. acaulis, P. loureiri, P. roebelenii, P. rupicola and P. sylvestris were considered as group 1 and P. theophrasti, P. dactylifera, P. reclinata, P. canariensis and C. humilis were considered as group 2 (Online Resource 2). Genetic relationships among P. theophrasti individuals at the population level and among the palm individuals at the species level were ascertained and graphically visualized to a better degree via factorial correspondence analysis (FCA) using the Genetix 4.05 software (Belkhir et al. 2004). This method allows a spatial representation of P. theophrasti populations and palm species with respect to the defined axis. A model-based clustering approach as implemented in the Structure software version 2.3.1 (Pritchard et al. 2000; Falush et al. 2003) was used to estimate the number of genetic clusters (K) and to fractionally assign individuals sampled from P. theophrasti populations and from palm species. We performed the analysis by setting number of clusters varying from 2 to 10. Each run was implemented with a burn-in period of 10,000 steps followed by 100,000 Monte Carlo–Markov chain replicates (Hubisz et al. 2009). Simulations were run using the admixture model with no prior population information. Method described in Evanno et al. (2005) was used to identify the appropriate number of clusters via the ad hoc statistic, ΔK, which is based on the second-order rate of change in the likelihood function with respect to successive values of K. Based on the Structure output file, the number of true clusters in the data (K) was determined using Structure Harvester v. 0.6.92 software program (Earl and von Holt 2012), which identifies the optimal K based on the ΔK and the posterior probability of the data for a given K (Evanno et al. 2005). The relationship between expected heterozygosity (He) and geographical distance (Km) from the nearest likely source population was tested using SigmaPlot 12.0 software (Systat Software, San Jose, CA, USA). Pearson’s correlation coefficients (rp) were computed to relate the genetic diversity measures (Na, Ne and He) of P. theophrasti populations to climatic (temperature, precipitation and relative humidity; data not shown) and geographic (altitude, latitude and longitude) data of the locations of these populations by using SPSS version 11 for Windows according to Steel and Torrie (1980). The climatic data were taken from the Turkish State Meteorological Service in Ankara as the mean values of 10 year period (2001–2010) recorded at the nearest meteorology stations to the locations of the populations.

Results

Genetic variation in Phoenix theophrasti populations

A total of 112 alleles were detected in 12 SSR loci, and the mean number of alleles per locus was found as 9.33, ranging from 2 (at locus mPdCIR016) to 13 (at loci mPdCIR015 and mPdCIR032) in the six P. theophrasti populations (Online Resource 1). Loci mPdCIR016 and mPdCIR090 did not amplify in Dimitri and Karaöz populations, respectively. While none of the loci were monomorphic in all populations, only the locus mPdCIR016 was found as monomorphic in Hurmalı, Eksera, Patara and Karaöz populations. The observed heterozygosity (Ho) values per SSR locus varied between 0.00 (mPdCIR016) and 0.76 (mPdCIR070), with a mean of 0.38. Higher values than the observed heterozygosities were detected for the expected heterozygosities (He) with values ranging from 0.06 (mPdCIR016) to 0.80 (mPdCIR070), and a mean of 0.64 (Online Resource 1). Of the 112 alleles detected in SSR loci, a total of 48 were revealed as private to populations. Gölköy population was the population having the highest number of alleles (72) and the highest number (23) and percentage (31.94%) of private alleles. Karaöz population was the population having the lowest number of alleles (25), of which three (12%) were private to this population. Among all populations, Dimitri population was the population having the lowest percentage (5.88%) of private alleles (Table 2). The percentage of polymorphic loci (P%) within populations varied between 75% (for Patara and Karaöz populations) and 100% (for Gölköy population) with a mean of 86.11%. While the lowest values of observed number of alleles (Na) (2.27) and effective number of alleles (Ne) (1.50) were detected in Karaöz population, the highest values of Na (6.00) and (Ne) (3.35) were detected in Gölköy population. The mean values of Na and Ne in the populations were 3.71 and 2.19, respectively. In the same way, the lowest and the highest values of observed (Ho) and expected (He) heterozygosities were found in Karaöz population (Ho = 0.24, He = 0.32) and Gölköy population (Ho = 0.46, He = 0.71), respectively. In the populations, the means of these genetic variability parameters were detected as Ho = 0.38 and He = 0.47. Lower Ho values than He values were found for all populations except for Dimitri and Patara for which very slightly higher values of Ho than He were detected (Table 3). Most of the loci were not in accordance with HWE. The departures from HWE were probably due to small (limited) population sizes or to the poor PCR performance of certain alleles (null alleles). A plot of the expected heterozygosity (He) against geographical distances (Km) from the most northwest population (Gölköy) is shown in Fig. 2. Our data indicate a strong correlation of heterozygosity reduction in northwest–southeast direction with a Pearson’s correlation coefficient of − 0.867, P < 0.05. This pattern advocates of northwest–southeast distribution of the species. While strong negative correlations between all three genetic variation parameters (Na, Ne and He) of P. theophrasti populations and relative humidity as a climatic factor were detected, very strong positive correlations between all three genetic variation parameters and latitude as a geographical factor were detected. The rest of the correlations between genetic data and environmental factors were not significant (Online Resource 3). For the 12 SSR loci screened in P. theophrasti populations, F-statistics, gene flow (Nm), within-population variation (HS) and total variation (HT) values were given in Online Resource 4. The mean of total genetic variation (HT) was 0.66. A high proportion of this variation, 0.44 (66.67%), was due to within-population genetic variation (HS). The mean of genetic differentiation coefficient (FST) was 0.34, and the mean gene flow (Nm) within a generation among the six populations was 0.49. The results of hierarchical AMOVA of all 12 SSR data revealed that only a limited degree of total variance (3.40%) was attributed to among groups, while the variation among populations within groups and within populations explained 23.55% and 73.05% of the total variation, respectively (Table 4). According to the pairwise FST analysis, the lowest differentiation was detected between the Dimitri and Hurmalı populations (FST = 0.1128), while the highest differentiations were detected between the Patara and Dimitri (FST = 0.5348), and Patara and Karaöz (FST = 0.5313) populations (Online Resource 5). Nei’s (1972) genetic distance (DN) coefficients ranged from 0.1822 to 1.5148 among the pairs of P. theophrasti populations. The minimum distance was detected between Hurmalı and Eksera populations, and the maximum distance was detected between Patara and Dimitri populations (Online Resource 5). A UPGMA dendrogram is presented in Online Resource 6 as a visual representation of the DN values of P. theophrasti populations based on SSR data. Three-dimensional FCA revealed the existence of three distinct groups of P. theophrasti populations. While Hurmalı, Dimitri, Eksera and Karaöz populations existed in the same group, Gölköy and Patara populations formed the other two groups separately (Online Resource 7). By the Structure Harvester v. 0.6.92 software program (Earl and von Holt 2012), a clear peak in the value of ΔK at K = 3 was determined. At K = 3, the six P. theophrasti populations were divided into three groups as in three-dimensional FCA analysis. Of the 86 P. theophrasti individuals analyzed, 61 (70.93%) (all the individuals sampled from Hurmalı, Dimitri, Eksera and Karaöz) shared > 70% membership with 1 of 3 clusters and were classified as members of that cluster, whereas 15 (17.44%) and 10 (11.63%) individuals were classified as members of the second and the third clusters, respectively. Among the 86 individuals, 41 (47.67%) showed admixture, between the second and third genetic clusters (Fig. 1). No clear geographic structure was detected using these approaches.

Table 2 Number of alleles detected in each SSR loci in Phoenix theophrasti populations
Table 3 Genetic variability estimates based on data from 12 SSR loci in Phoenix theophrasti populations
Fig. 2
figure 2

Correlation between the heterozygosities (He) of the populations and their geographical distances (Km) in northwest–southeast direction

Table 4 Values of AMOVA test partitioning of SSR variability in Phoenix theophrasti populations

Relations of Phoenix theophrasti with other palm species

Analysis of ten palm species in terms of 12 SSR loci resulted in a total of 166 alleles, so that the mean number of alleles per locus was 13.83. The number of alleles in a locus ranged from 5 (at locus mPdCIR016) to 21 (at locus mPdCIR048) (Online Resource 1). Loci mPdCIR085 and mPdCIR093 in C. humilis, locus mPdCIR078 in C. humilis and P. roebelenii, locus mPdCIR090 in P. sylvestris and P. acaulis, locus mPdCIR016 in P. dactylifera, P. canariensis and P. acaulis, locus mPdCIR070 in P. reclinata, P. rupicola, P. roebelenii, P. acaulis and C. humilis did not amplify. None of the loci were detected as monomorphic in all species. However, only locus mPdCIR057 in P. roebelenii, locus mPdCIR016 in P. sylvestris, P. reclinata and P. rupicola and loci mPdCIR078 and mPdCIR093 in P. sylvestris were detected as monomorphic. The observed heterozygosity (Ho) values per SSR locus varied from 0.12 (mPdCIR016) to 0.82 (mPdCIR010), with a mean of 0.52, while the expected heterozygosity (He) values per SSR locus varied between 0.73 (mPdCIR016) and 0.93 (mPdCIR010, mPdCIR078 and mPdCIR085) with a mean of 0.88 (Online Resource 1).

Among the 166 alleles detected in SSR loci, a total of 74 were found as private to species. Phoenix loureiri was the species for which the highest number of alleles (51) was detected. Among these 51 alleles, 16 (31.37%) were private to this species. On the other hand, the lowest number of alleles (22) was detected for P. acaulis, of which five (22.73%) were found private to that species. Among the ten palm species analyzed, P. dactylifera had the lowest percentage (6.45%) of private alleles (Online Resource 8).

The percentage of polymorphic loci (P%) in the species varied between 66.67% (for P. sylvestris and C. humilis) and 100% (for P. loureiri) with a mean of 82.50%. The means of observed (Na) and effective (Ne) number of alleles in the species were 3.20 and 2.61, respectively. Among the species, P. loureiri had the highest Na (4.25) and Ne (3.43) values. On the other hand, the lowest Na (2.44) value was detected for P. acaulis and the lowest Ne (2.21) value was detected for both P. acaulis and C. humilis species. For the species, the mean of Ho was 0.53 and the mean of He was 0.64. The minimum (0.34) and the maximum (0.66) values of Ho were detected for P. theophrasti and P. dactylifera, respectively, while the minimum (0.52) and the maximum (0.78) values of He were detected for P. sylvestris and P. loureiri, respectively. Lower Ho values than He values were found for all species except for P. dactylifera and C. humilis, for which very slightly higher values of Ho than He were detected (Table 5). The mean of total genetic variation (HT) for the species was 0.87, and a higher proportion of this variation, 0.55 (63.22%), was due to within-species genetic variation (HS). The mean of genetic differentiation coefficient (FST) was 0.45, and the mean gene flow (Nm) was 0.38 (Online Resource 9). In the hierarchical AMOVA test, 25.52% of total variation was attributed to among species within groups and 77.15% of total variation was attributed to within species (Table 6). In the pairwise FST analysis, the lowest differentiation was detected between P. reclinata and P. roebelenii species (FST= 0.0825), while the highest differentiation was detected between P. theophrasti and C. humilis species (FST= 0.4261). The lowest differentiation of P. theophrasti was from P. dactylifera (FST = 0.1932) (Online Resource 10). For the pairs of palm species, the Nei’s (1972) genetic distance coefficients (DN) ranged from 0.8995 (between P. reclinata and P. rupicola) to 5.2943 (between P. acaulis and C. humilis). Phoenix roebelenii and C. humilis were found as the least (DN = 0.9740) and the most (DN = 2.5504) genetically distant species from P. theophrasti, respectively (Online Resource 10). As a visual representation of DN values between the species based on SSR data, a UPGMA dendrogram was shown in Online Resource 11.

Table 5 Genetic variability estimates based on data from 12 SSR loci in the palm species studied
Table 6 Values of AMOVA test partitioning of SSR variability in palm species

Three-dimensional FCA resulted in five distinct groups of palm species. While Groups A and E consisted solely of P. theophrasti and C. humilis, respectively, Group B formed from P. dactylifera and P. roebelenii, Group C included P. canariensis, P. reclinata, P. rupicola and P. acaulis, and Group D comprised of P. sylvestris and P. loureiri (Online Resource 12). A clear peak in the value of ΔK at K = 6 was determined in the Structure Harvester v. 0.6.92 software program (Earl and von Holt 2012). At K = 6, the ten palm species formed six clusters. Phoenix theophrasti formed the cluster 1. Cluster 2 included P. dactylifera and P. roebelenii, cluster 3 consisted of P. canariensis and P. rupicola, cluster 4 and 5 formed only from P. reclinata and C. humilis, respectively, and cluster 6 included P. loureiri and P. acaulis. Individuals of P. sylvestris showed admixture (Fig. 3).

Fig. 3
figure 3

Genetic clustering of ten palm species assessed from structure analysis. Six clusters were inferred: Cluster 1 is yellow, cluster 2 is red, cluster 3 is blue, cluster 4 is pink, cluster 5 is green, and cluster 6 is aqua in color. The vertical coordinate of each cluster indicates the membership coefficient for each individual. Bars with different colors represent genotypes that show admixture character. (1: Phoenix theophrasti, 2: P. dactylifera, 3: P. canariensis, 4: P. reclinata, 5: P. rupicola, 6: P. roebelenii, 7: P. sylvestris, 8: Chamaerops humilis, 9: P. loureiri and 10: P. acaulis)

Discussion

Tree species have experienced various environmental changes during their evolutionary and ecological histories. As was the case in glacial periods, some environmental changes may have lasted as long as 100,000 years or they may have lasted only a decade. These changes may have been gradual or sudden (Hamrick 2004). To present day, only a few natural populations of P. theophrasti survived the last glaciation in isolated areas and remained as relict populations (Médail and Quézel 2003). Since these populations have been facing intense human activities due to their specific distribution areas, this tree species has likely experienced both long and short, as well as gradual and sudden, environmental changes during its evolutionary and ecological history. Recently, enlarged consideration has been given to the effects of environmental changes on tree species because of a tremendous increase in harmful human activities. For the future of tree species, the degree and the rate of changes are important. In this respect, we can mention three possible outcomes. First, changes may be of such a degree that tree species will not have enough genetic variation to adapt to the changed environment (Critchfield 1984; Davis and Zabinski 1992). Second, compared to their long generation times, changes may happen so quickly that trees could not adapt rapidly enough to the changes (Davis and Shaw 2001). Third, trees may contain sufficient genetic diversity but may not be able to disperse into new suitable habitats rapidly enough to overcome environmental change (Clark 1998). All three could be applicable for P. theophrasti, especially the second due to its long generation time; P. theophrasti requires indeed 4–5 years to reach the first flowering.

The level and pattern of genetic diversity in a native plant species can be changed by human activities, leading to loss of genetic diversity that may at the end reduce the evolutionary potential of a species to respond to environmental changes (Ellstrand and Elam 1993). It was shown by Godt et al. (1996) and Godt and Hamrick (2001) that in many species, genetic diversity is directly related to population size and levels of genetic diversity may affect individual fitness and potential population persistence (Newman and Pilson 1997; Fischer and Matteis 1998). First suggested by Stebbins (1942), it has been almost universally assumed that endemic and rare taxa show significantly less genetic variability compared with widespread and abundant ones (Hamrick and Godt 1989; Falk and Holsinger 1991; Affre et al. 1997; Garrido et al. 2012).

The current study has provided information about the relative amounts of genetic variation present within and among six Turkish P. theophrasti populations by using SSRs as a nuclear marker system, for the first time in the literature. A higher proportion (66.67%) of genetic variation was due to differences within the populations, while the remaining relatively small portion of variation was distributed among the populations. Gene flow among the six studied populations was consistently low, the rate being Nm = 0.49 migrants per generation, which are below the critical value of Nm = 1.0, indicating genetic differentiation among populations due to genetic drift. In support of this finding, high genetic differentiation coefficient FST (0.34) was assessed among populations. Therefore, genetic differentiation among Turkish P. theophrasti populations appeared critically high with a critically low level of gene flow among them. According to Hamrick and Godt (1996), the mean value of FST for outcrossing monocot plants was 0.16, for long-lived perennial endemic plants, it was 0.15, and for long-lived perennial rare plants, it was 0.13. In the evaluation of combined categories of life forms and breeding systems by the same authors, the mean value of FST for long-lived perennial outcrossing plants was 0.09. Having the same combined life history traits, the FST value determined for P. theophrasti populations in our study is about 2–4 times higher than these mean values, as well as higher than the values determined for P. canariensis (FST = 0.25) and P. dactylifera (FST = 0.20) (González-Pérez et al. 2004). Peñaloza-Ramírez et al. (2016) showed the consequences of habitat fragmentation on the genetic structure of Chamaedorea alternans (Arecaceae) palm populations in the tropical rain forests of Mexico, with FST values above 0.25 and Nm values below 1.0 reported for fragmented populations of the species. Not only palm but also other plant species show the same pattern of high FST and low Nm values for their fragmented populations. For example, the FST value was detected as 0.39 and Nm value was detected as 0.57 for 35 populations of Rhodiola dumulosa (Crassulaceae), which is a perennial diploid species distributed as fragmented populations across northern, central and northwestern China (Hou and Lou 2011). As one of the rare and endangered plant species of Europe, the present range of Pulsatilla patens in Europe is highly fragmented and the size of the populations has been dramatically reduced in the past 50 years. For 29 analyzed populations of this species, pairwise FST values up to 0.60 were detected (Szczecińska et al. 2016). Generally, observed heterozygosity values were less than expected heterozygosity in P. theophrasti populations, as was the overall mean of Ho = 0.38 and the overall mean of He = 0.47. Mean FIS across all SSR loci was 0.13. These results are the indication of the presence of an excess of homozygotes in the populations. With very slightly higher Ho values than He, Dimitri and Patara populations were the exceptions to these findings, indicating the presence of very low levels of selection in the favor of heterozygotes in these populations. Population fragmentation and isolation and lowered population densities may result in modified gene flow patterns due to reduction in pollen and seed movement between populations, increasing the inbreeding within populations (Cruse-Sanders and Hamrick 2004). The level of gene flow among populations is the key element for anticipating the effects of fragmentation. In case of high measures of gene flow, little variation will be lost at the separate fragment. On the other hand, in case of gene flow rates below pre-fragmentation levels, genetic differentiation among fragmented populations will increase and genetic diversity within fragments will decrease. If there are low effective numbers of parents, inbreeding may also increase within fragments (Hamrick 2004).

Although these findings strongly indicate a severe population fragmentation of Turkish P. theophrasti populations, the FCA and structure analyses showed that Hurmalı, Dimitri, Eksera and Karaöz populations are genetically close. In fragmented populations, the lack of genetic differences or minor differences between fragments could be explained by two hypotheses. The first hypothesis asserts that the distribution of genetic diversity within and among populations indicates historical gene flow processes which led to the fragmentation of larger populations. The second hypothesis postulates that geographically close populations are more efficiently linked by gene flow than populations separated by a longer distance (Szczecińska et al. 2016). In case of P. theophrasti, both explanations could be effectual when we consider the ecological history of the species as a relict population that survived the last glaciation in isolated areas, and the geographical proximity of the current populations, though Karaöz is geographically distant to Hurmalı–Dimitri–Eksera populations present in Datça Peninsula as fragmented populations. The pattern observed for the Karaöz population can be due to an introduction of individuals from Datça populations in Karaöz region (maybe human-mediated or linked to some geological/ecological events).

On Crete, there are eight subpopulations of P. theophrasti and the largest subpopulation is protected under Greek Law (IUCN Red List, version 2018-1). Therefore, it has been assessed under Lower Risk/Near Threatened (LR/NT) category in the Red List of IUCN since 1998. However, it is known that some populations of the Cretan P. theophrasti have substantially decreased or disappeared altogether. In the past, the main threat for this species was the clearance of land for cultivation purposes, but today, the greatest threats arise from tourism development. Increased exploitation of underground water is another important threat for the Cretan populations of the species (Johnson 1996). Similarly, considerable pressure to the P. theophrasti areas from the tourism sector, road building activities, insect epidemics due to exotic palm species introduced to the distribution areas of P. theophrasti and climate change to a drier and warmer climate in the Mediterranean region are the main current and most probable future threats for the populations of the species in Turkey (Küçükala et al. 2008). Considering all these threats for both Cretan and Turkish P. theophrasti populations, it is appropriately annotated in the Red List of IUCN (version 2018-1) as a species whose assessment “needs updating.” Supporting this requirement annotated in Red List, the results obtained in this study showed that P. theophrasti must urgently be included in the Red List of IUCN under the critically endangered (CR) category and “in situ” combined with “ex situ” conservation precautions should urgently be taken for Turkish P. theophrasti populations. Because of their known advantages and disadvantages to each other, “in situ” and “ex situ” conservation strategies are deemed complementary (UNCED 1992; FAO 1998). Several approaches combining ex situ and in situ conservation were proposed in the past, but none are satisfying as a general conservation strategy (Volis et al. 2009). Although only few conservation programs employ both approaches together successfully (Maxted et al. 1997), we think that application of “in situ” combined with “ex situ” conservation precautions will be very useful for this valuable tertiary relict species, considering the special ecological and anthropogenic structure of the region.

As one of the outcomes of our study, the date palm (P. dactylifera L.) was found as the genetically closest species of P. theophrasti. With its high nutritive, economic and social values, the date palm is an important cultivated tree crop in North Africa, India, Spain and especially in the Middle East where it has been domesticated for at least 5000 years (Al-Khalifah and Askari 2003; Akkak et al. 2009). Date palms are generally propagated by vegetative means to make it possible to clone individuals with economically important traits although seeds of the species are breeding material with long backcrossing cycles (Al-Khalifah and Askari 2003; Elshibli and Korpelainen 2008). Widespread use of clonal propagation of elite cultivars with desired traits in the growing regions of this crop reduces the genetic diversity of the cultivars, making them open to biotic and abiotic stress factors. Diseases and drought are the main biotic and abiotic factors limiting the date palm production in the growing regions (Oihabi 2000; Ahmed et al. 2001; Elshibli and Korpelainen 2008). Additionally, the usage of a very limited number of elite cultivars put the species in a dangerous situation because of severe genetic erosion (Rhouma 1994; Zehdi et al. 2004). For these reasons, preservation and use of the genetic information of P. theophrasti are very valuable for the breeding purposes of the date palm.

Certain morphological differences detected between the Gölköy palm and P. theophrasti in other locations have resulted in the classification of Gölköy palm as P. theophrasti ssp. Gölköy. However, to what extent these differences may indicate the Gölköy palm indeed as a subspecies of P. theophrasti, or rather as a distinct species, a variety of P. dactylifera or else, remains unclear. Findings of our study showed that the pairwise genetic distance DN and genetic differentiation coefficients FST of the Gölköy population with the other populations were not so different from the pairwise values of the same coefficients between other populations. On the other hand, the Patara population, having about two times higher pairwise values than the Gölköy population for these parameters, has been accepted as true P. theophrasti according to morphological criteria. Therefore, the classification of the Gölköy population as P. theophrasti ssp. Gölköy seems to be not proper on molecular level.

Genetic variation indices Na, Ne and He of P. theophrasti populations revealed significant negative correlations with relative humidity but significant positive correlations with latitude. The absence of significant correlations of these genetic parameters with other environmental factors somewhat contradicts recent studies (Jordano and Godoy 2000; Chen et al. 2009; Goto et al. 2009; Hübner et al. 2009; Thiel-Egenter et al. 2009; Manel et al. 2012).

Although the value of genetic differentiation coefficients (FST = 0.45) indicated that the genetic differentiation among the species is high, the low value of inbreeding coefficient (FIS = 0.05) indicated the presence of random mating among the species. Since their geographic distribution areas are close to each other, we expected a specific genetic closeness between C. humilis and P. theophrasti. However, C. humilis was found as the most genetically distant species from P. theophrasti. Even though C. humilis belongs to a different genus, genetic distances between C. humilis and Phoenix species were found not too different from the genetic distances among Phoenix species themselves. Therefore, we reached the conclusion that there is a requirement for re-evaluation of the taxonomy of palms, especially by using detailed molecular data.