Introduction
The genus Pistacia L. belongs to the Anacardiaceae family originating from Central Asia (Parfitt and Badenes 1997). It comprises eleven accepted species, two subspecies and one hybrid taxon, according to Plants of the World Online (POWO 2023). The genus has a wide but disjunct native distribution including Central to North America, the Mediterranean to Central Asia and northeast to east Africa (AL-Saghir and Porter 2012, Xie et al. 2014, POWO 2023). Representatives of the genus are dioecious, deciduous trees or evergreen shrubs, considered to be mainly xerophytes growing in arid and/or semi-arid habitats. Although the genus does not comprise many species, some of them are well known and have been used by mankind since ancient times due to their valuable products (Golan-Goldhirsh 2009). For example, P. vera L. (the pistachio tree) is a highly economically important species due to its edible seeds and for this reason it has been commercially cultivated and used for centuries in many parts of the world (Whitehouse 1957, Zeng et al. 2019, Kafkas et al. 2023). Other examples include the age-old use of the mastic gum and resins of Pistacia trees for various purposes such as medicine, food, aroma and cosmetics (Golan-Goldhirsh 2009).
The Mediterranean is considered to be the contemporary diversification centre of Pistacia (Xie et al. 2014). Six species of the genus are listed for the Euro-Mediterranean area ( P. atlantica Desf., P. eurycarpa Yalt., P. khinjuk Stocks, P. lentiscus L., P. terebinthus L. and P. vera), with only two native species currently present in the eastern Adriatic according to Euro+Med PlantBase (Euro+Med 2006-2023): P. lentiscus and P. terebinthus. Interestingly, although the hybrid taxon P. × saportae Burnat between parental species P. lentiscus and P. terebinthus is listed as native for the area of western Mediterranean, Morocco, Palestine, Algeria, Cyprus and Turkey according to POWO (2023), it is not currently recorded in the Euro+Med PlantBase. Pistacia × saportae was described by Burnat in 1896 as a probable species from the Maritime Alps. Later, it was suggested and confirmed to be an interspecific hybrid taxon ( P. lentiscus × P. terebinthus) in view of its morphological, wood anatomy and genetic data (Zohary 1972, Grundwag and Werker 1976, Werner et al. 2001, Yi et al. 2008, AL-Saghir and Porter 2012), reported in many Mediterranean countries so far. According to the Flora Croatica Database – FCD (Nikolić 2023), in the Croatian flora (eastern Adriatic) the genus Pistacia is represented by four taxa: P. lentiscus, P. terebinthus, their hybrid P. × saportae and sporadically cultivated P. vera. In addition, Radić (1985) described P. calcivora Radić as a new endemic taxon from Mt. Biokovo based on morphology (On-line Suppl. Fig. 1). Since the taxonomic status of this taxon is however doubtful and currently unresolved, it is not listed in the national flora i.e. FCD. Pistacia × saportae was reported and inserted in FCD for the first time in the Croatian flora by Jeričević et al. (2014) from the island of Korčula, the identification being based on morphology, although it was mentioned earlier for the area of Mt. Biokovo by Radić (1985) (On-line Suppl. Fig. 1).
The taxonomy and phylogeny of the genus Pistacia has been assessed on several occasions on the basis of morphological and genetic (RAPD, AFLP, RFLP, nuclear and plastid sequences) data (Parfitt and Badenes 1997, Golan-Goldhirsh et al. 2004, Kafkas 2006, Yi et al. 2008, AL-Saghir and Porter 2012, Xie et al. 2014), however there are still unresolved issues about some taxa, taxa number and species boundaries due to high morphological variation both among and within taxa (Golan-Goldhirsh 2009). This is further complicated by common hybridization between some of the representatives, often forming interspecific hybrids (Zohary 1952, Morgan et al. 1992, Yi et al. 2008, Aznarte-Mellado et al. 2014). There are several different taxonomic classifications of the genus. The first and probably most comprehensive one based on morphology was carried out by Zohary (1952) who recognized 11 species belonging to four sections: Sect. Lentiscella Zoh. ( P. mexicana HBK, P. texana Swingle), Sect. Eu Lentiscus Zoh. ( P. lentiscus, P. saportae, P. weinmannifolia Poiss. ex Franch.), Sect. Butmela Zoh. ( P. atlantica) and Sect. Eu Terebinthus Zoh. ( P. terebinthus, P. palaestina Boiss. P. khinjuk, P. vera, P. chinensis Bunge). Later, Parfitt and Badenes (1997) carried out the first molecular phylogeny of the genus (based on plastid DNA) including ten Pistacia species and proposed the division of the genus into only two sections: Sect. Lentiscus (including all evergreen species with paripinnate leaflets) and Sect. Terebinthus (deciduous species with imparipinnate leaflets). Yi et al. (2008) provided a more recent, updated molecular phylogeny of the whole genus based on five different data sets (sequences of two nuclear and three plastid regions) and including 11 species and one putative hybrid taxon P. × saportae. They confirmed that the genus is monophyletic, as well as the hybrid origin of P. × saportae suggesting that P. lentiscus is the maternal and P. terebinthus paternal taxon. Their phylogeny is largely consistent with the most recent taxonomic revision of the genus based on morphology (up to our knowledge) carried out by AL-Saghir and Porter (2012), which recognized a total of 13 taxa divided into two sections (sect. Pistacia Zoh. and sect. Lentiscella): nine species ( P. atlantica, P. chinensis, P. eurycarpa, P. khinjuk, P. terebinthus, P. vera, P. lentiscus, P. mexicana, P. weinmannifolia), five subspecies ( P. chinensis subsp. chinensis, P. chinensis subsp. falcata (Becc. ex Martelli) Rech. f., P. chinensis subsp. integerrima (J.L. Stew. ex Brandis) Rech. f., P. lentiscus subsp. lentiscus, P. lentiscus subsp. emarginata (Engl.) AL-Saghir) and one hybrid taxon ( P. × saportae) . Finally, the most recent biogeographic history and diversification of Pistacia genus was carried out by Xie et al. (2014), using an expanded set of sequences of multiple nuclear and plastid loci. The latter study better resolved the relationships within the genus compared to Yi et al. (2008) and did not support previous divisions of the genus into four (Zohary 1952) or two sections (Parfitt and Badenes 1997, AL-Saghir and Porter 2012). The authors estimated that the genus has diverged at 37.60 mya, which is much later than some previous estimations (Parfitt and Badenes 1997).
Over the past decades, the genetic relationships among some specific Mediterranean Pistacia taxa at the molecular level were assessed across different geographic areas of interest using various molecular markers such as RFLP (Parfitt and Badenes 1998), RAPD and/or AFLP (Kafkas and Perl-Treves 2002, Katsiotis et al. 2003, Golan-Goldhirsh et al. 2004, Karimi et al. 2009, Shanjani et al. 2009, Iranjo et al. 2016), ITS (Labdelli et al. 2022), SSR (Vendramin et al. 2009), and SAMPL (Karimi and Kafkas 2011) or in specific countries like Syria (Basha et al. 2007), Turkey (Kafkas and Perl-Treves 2001, Guney et al. 2021), Afghanistan (Kafkas et al. 2006), Iran (Mirzaei et al. 2006) etc. However, no such comparative study exists for the native Pistacia taxa in the eastern Adriatic region. Thus, building on this rich legacy to contribute to ongoing discussions about genetic relationships among Pistacia taxa, the aim of our research was to determine relationships among four putative indigenous taxa of the genus Pistacia in the eastern Adriatic ( P. lentiscus, P. terebinthus, P. × saportae and P. calcivora) using AFLP molecular markers and leaf morphology. Specifically, we aimed to determine the presence of the hybrid taxon P. × saportae and to validate the taxonomic status of the putative taxon P. calcivora in the eastern Adriatic at the molecular and morphological levels. We chose AFLP markers because they proved to be efficient in identifying interspecific hybrids between closely related species (Oberprieler et al. 2022, Zalewska-Gałosz et al. 2023) and have been previously used to address genetic relationship between selected Pistacia species (Karimi et al. 2009). In addition, relationships between our taxa of interest were evaluated based on leaf morphology to complement the genetic data and to provide the robustness of the taxonomic treatment because leaf characteristics have been traditionally used as the main taxonomic characters for Pistacia identification (Parfitt and Badenes 1997).
Plant material and study area
We collected plant material from a total of 13 natural Pistacia populations (116 individuals) corresponding to four putative indigenous taxa of the genus Pistacia present in the eastern Adriatic: P. lentiscus, P. terebinthus, P. × saportae and P. calcivora (Tab. 1) . Our study area belongs to the Mediterranean phytogeographic region characterized by the typical Mediterranean climate with hot dry summers and mild wet winters. An an essential part of the Mediterranean vegetation, Pistacia taxa in the eastern Adriatic thrive in dry, warm, sunny, rocky habitats within the Adriatic coastal and/or coastal-montane belt where they are frequently exposed to prolonged droughts and high temperatures.
Tab. 1. Sampled taxa and populations of the genus Pistacia in the eastern Adriatic coast. Abbreviations: n – number of individuals for genetic analyses, n (M) – number of individuals for morphological analyses, x and y – geographical coordinates.
We collected three putative hybrid populations of P. × saportae (a total of 21 individuals) from the localities where the hybrids were previously reported or assumed based on their morphology: islands of Šolta, Korčula and Vis (Fig. 1). We also collected three populations (a total of 30 individuals) of the putative taxon P. calcivora from three different localities from the type locality on Mt. Biokovo (Tab. 1, On-line Suppl. Fig. 2). Both of these putative taxa were identified in the first place in the field based on leaf morphology by an expert botanist. In addition, at each sampling site of P. × saportae we collected individuals of both parental species: a total of 32 individuals of P. lentiscus and 30 individuals of P. terebinthus (Tab. 1, On-line Suppl. Fig. 2). One additional population of P. lentiscus was also collected from Mt. Biokovo for comparison with the putative P. calcivora (Tab. 1). Collected plant material for the DNA extraction was immediately stored in silica-gel. Finally, in each studied population we collected leaves for the morphological analyses from the same individuals sampled for the genetic analyses. From each individual tree or shrub we collected 10 leaves from short shoots which were herbarised in the field until the morphometric measurements. However, some of the collected herbarised plant material decayed prior to the measurements, thus a few individuals had to be excluded from the morphological analyses (Tab. 1). Herbarium specimens are deposited in ZAGR herbarium.
Fig. 1. Leaf variability of Pistacia terebinthus (left), P. × saportae (middle) and P. lentiscus (right) from the island of Vis in the eastern Adriatic (Photo: S. Bogdanović).
DNA extraction and AFLP fingerprinting
Total genomic DNA was extracted from 20 mg of silica-gel dried leaf tissue of each individual in the study using the DNAeasy plant mini kit (QIAGEN). The manufacturer’s protocol was followed and slightly modified by adding 2-mercaptoethanol and 1% PVP (polyvinylpyrrolidone) to increase quality and concentration of DNA extract. DNA quality and concentration of 116 samples was verified using Nanophotometer P330 (Implen®, Munich, Germany). For AFLP protocol we followed Vos et al. (1995) with modifications described in Carović-Stanko et al. (2011).
We used the EcoRI and MseI (= Tru1I) restriction enzymes and five primer combinations for the selective amplification: FAM-EcoRI-ACA, NED-EcoRI-AGA, VIC-EcoRI-ACG, PET-EcoRI-ACC and Tru1I-CGA. For preselective amplification the PCR program was as follows: 2 min at 94 °C followed by 20 cycles: at 94 °C for 20 s, 56 °C for 30 s, 7 °C for 2 min, and 60 °C for 30 min. For selective amplification we used the PCR program from Carović-Stanko et al. (2011). Amplified PCR fragments were separated by capillary electrophoresis on an ABI3730 XL automated sequencer (Applied Biosystems) by Macrogen Europe (Amsterdam, The Netherlands) using GeneScan-500 LIZ size standard (Applied Biosystems) for size calibration. We identified the AFLP alleles using GeneMapper® 5.0 software (Applied Biosystems) within the size range of 50–600 base pairs (bp) which were then exported as binary matrix for subsequent analyses with fragments scored as present (1) or absent (0). In total, we obtained 1502 AFLP alleles for the 116 individuals. Based on the DNA quality and quality of chromatograms in GeneMapper, three individuals were excluded from the further analyses. We then estimated the error rate based on eight replicated samples and selected the final set of AFLP alleles for the subsequent analyses using ScanAFLP ver. 1.2 R script following Herrmann et al. (2010). This resulted in a dataset of 584 AFLP loci for a total of 113 individuals used for all the subsequent statistical analyses, with a mean error rate of 1.79%.
AFLP analyses
We first estimated within-population genetic diversity for each investigated population by calculating the proportion of polymorphic markers (%P), Nei’s gene diversity (HE, Nei 1987) and the frequency down-weighted marker values (DW; Schönswetter and Tribsch 2005) which indicate the genetic rarity of a population using AFLPdat R scripts (Ehrich 2006) in R version 4.1.3. In addition, we used AFLPdat to prepare the input file for the Bayesian model-based clustering method implemented in the STRUCTURE software 2.3.4 (Pritchard et al. 2000), which was used to assess the genetic structure of the investigated populations. We used the following parameters for the STRUCTURE analysis: number of genetic clusters (K) was set from 1 to 10, with 10 replicate runs per each K, each run consisting of a burn-in period of 10,000 Markov Chain Monte Carlo (MCMC) iterations followed by 90,000 MCMC iterations. We used an admixture model of individual ancestry and correlated allele frequencies. To determine the optimal number of K we used the ΔK method of Evanno et al. (2005). The STRUCTURE results were subsequently processed (averaging replicates), summarized and visualized using the STRUCTURE HARVESTER 0.6.94 software (Earl and von Holdt 2012) and CLUMPAK (Kopelman et al. 2015). Hybrid individuals with putative mixed ancestry were defined as those with an assignment probability to any of the genetic clusters below 0.7 (Q ≤ 70%). Next, we calculated pairwise genetic distances between all individuals based on Dice coefficient (Dice 1945; Nei and Li 1979) and conducted a principal coordinate analysis (PCoA) based on the resulting distance matrix using PAST 3.22 (Hammer et al. 2001). Finally, to investigate the phylogenetic relationships between individuals and putative taxa under study, we constructed a NeighborNet diagram of individuals based on Dice coefficient transformed into distance matrix using SplitsTree 4 (Huson and Bryant 2006). The reliability of the phylogenetic network was assessed by constructing a neighbor-joining (NJ) tree based on the same distance matrix using 1000 bootstrap (BS) repetitions, both implemented in PAST.
Morphological analyses
We carried out morphological analyses on a total of 13 Pistacia populations and 109 individuals (Tab. 1). After being herbarised and fully dried, leaves were scanned, using an A3-format scanner (MICROTEK ScanMaker 9800XL), at resolution of 600 dpi (.TIF). Obtained image files were further analysed in the software package WinFolia PRO (Regent Instruments Inc., QC, Canada). Using the interactive measurement option, the length of rachis (RL) and leaf petiole (PL) were measured on each leaf. Afterwards, on each leaf one lateral leaflet in the middle part of the leaf was selected and measured using the leaf morphology option. In total, eight leaflet morphometric traits were measured, with accuracy of 0.1 mm: leaflet area (LA); leaflet length (LL); maximum leaflet width (MLW); leaflet length, measured from the leaflet base to the point of maximum leaflet width (PMLW); leaflet width at 90% of leaflet length (LW90); angle enclosed by the main leaflet vein (the centre of leaflet blade) and the line connecting the leaflet base to a set point on the leaflet margin, at 10% (LA10) of total leaflet length; and form coefficient (FC). In addition, the total number of leaflets per each leaf (NL) was counted.
To assess the possibility of conducting parametric tests, the symmetry and homoscedasticity of data were verified (Sokal and Rohlf 2012). Assumptions of normality were checked using the Shapiro–Wilk test, and the assumption of homogeneity of variance was checked using Levene’s test. Arithmetic mean (M), standard deviation (SD) and coefficient of variation (CV/%) were calculated for the particular trait for each putative taxon, population and obtained genetic group (see results). Differences in leaf morphology between the studied taxa and populations were determined by Kruskal-Wallis test. In addition, statistically significant differences between all pairs of putative taxa were identified using Wilcoxon rank sum test, at P ≤ 0.05. Descriptive statistics and nonparametric tests were carried out using the Statistica software package ver. 13 (StatSoft Inc., Tulsa, OK, USA).
We used two multivariate analyses to examine the morphological variability of the studied taxa. First, we used principal component analysis (PCA) to calculate the principal components across all individuals, and all studied morphological traits. PCA biplot was constructed based on the first two principal components showing all the analysed individuals and traits without any predefined groups. In addition, we performed discriminant analysis (DA) to evaluate the utility and importance of the studied morphological traits by determining which traits were most useful for maximizing differentiation between the identified genetic groups (see results). Individual assignment to the specific genetic group was based on the results of the STRUCTURE analysis. We then used classification discriminant analysis to determine the proportion of individuals that were correctly classified into the identified genetic groups. All multivariate statistical analyses were carried out using the “MorphoTools” R scripts in R v.3.4.3 (R Core Team 2017) by following the manual of Koutecký (2015).
Population genetic diversity
Different genetic diversity parameters of each population showed a similar trend. The proportion of polymorphic markers (%P) per population ranged from 1.37% (population P5) to 26.71% (population P6), with an average of 18.81% (Tab. 2). Population P6 had the highest HE, and population P5 the lowest. The average DW value was 566.72, with the highest value detected in population P7 and the lowest once again in population P5 (Tab. 2). Overall, populations P6 ( P. terebinthus, Korčula), P9 ( P. terebinthus, Vis), and P3 ( P. terebinthus, Šolta) had the highest genetic diversity, while populations P7 (putative hybrid population from Vis), P13 ( P. lentiscus, Biokovo), and P2 (putative hybrid population from Šolta) had the highest frequency of rare alleles (DW). Overall, the population with the lowest genetic diversity was P5 ( P. × saportae, Korčula). However, in the case of P5 and P13, genetic diversity values should be treated with caution due to the low sample size of these populations.
Tab. 2. Genetic diversity parameters of the investigated Pistacia populations based on AFLP markers. %P – Proportion of polymorphic markers, HE – Nei’s gene diversity (Nei 1987), DW – frequency down-weighted marker values (rarity index). Genetic group – population assignment to each of the two genetic clusters identified by STRUCTURE analysis.
Genetic structure and phylogenetic relationships
The results of the PCoA analysis based on the Dice distance matrix between all individuals revealed two distinctly separated groups of individuals (Fig. 2a). One group consisted of P. terebinthus individuals and all individuals presumed to belong to the P. calcivora. The other group consisted of P. lentiscus individuals and most of the individuals initially identified based on morphology as P. × saportae. Closer to the P. lentiscus group, four potentially hybrid individuals (two from population P5 and two from P7) were clearly distinguished (Fig. 2a). The first two PCo axes explained 49.8 percentage of variation (Fig. 2a).
Based on the Bayesian clustering method implemented in STRUCTURE, the optimal number of genetic clusters (source populations) that can explain the genetic structure of the investigated populations was at K = 2 following ΔK (Evanno et al. 2005). The first genetic cluster included all individuals of P. terebinthus and the individuals that were originally assigned to the potential taxon P. calcivora based on leaf morphology. The second cluster corresponded to all individuals of P. lentiscus and included the majority of the putative hybrid individuals initially identified as P. × saportae based on leaf morphology. However, the same four hybrid individuals originating from the islands of Korčula and Vis were clearly identified based on their assignment probabilities to the two parental genetic clusters (Q = 0.66 – 0.71 to the P. lentiscus cluster), confirming the presence of P. × saportae at the molecular level (Fig. 2b).
Finally, the phylogenetic network based on Dice distance matrix among all individuals was largely congruent with the STRUCTURE and PCoA results and confirmed the presence of two main distinctly separated groups of individuals: one group corresponding to P. lentiscus and the other to P. terebinthus which also included all putative P. calcivora individuals. The above-mentioned hybrid individuals of P. × saportae were placed between the two parental taxa and positioned themselves closer to the P. lentiscus (Fig. 2c). The two major groups were also supported with high BS values (BS = 96 and 100), as well as hybrid individuals of P. × saportae which were clearly separated from both parental taxa (BS = 100).
Fig. 2. Genetic relationships and genetic structure of the investigated Pistacia taxa (13 populations and 113 individuals) in the eastern Adriatic based on AFLP data. Hybrid individuals corresponding to P. × saportae are indicated with an asterisk. a – Principal coordinate analysis (PCoA) scatter-plot based on Dice distance matrix visualizing the genetic relationships among investigated Pistacia individuals and taxa. Individual assignment to the specific taxon is based on the initial a priori identification based on leaf morphology. b –population genetic structure based on the Bayesian STRUCTURE analysis at optimal K = 2. Each column represents an individual, and the membership coefficient (Q) indicating the assignment probabilities to each of the two genetic clusters is depicted on the scale from 0 – 100%. Dark green corresponds to genetic cluster of P. terebinthus, orange to the genetic cluster of P. lentiscus. Population numbers correspond to Tab. 1 and Tab. 2. c – neighbour-net diagram based on Dice distance matrix depicting the phylogenetic relationships among individuals. Colours of the points are as in a). Thick outlines correspond to the two genetic clusters inferred by STRUCTURE analysis in panel b). Dark green – genetic cluster of P. terebinthus, orange – genetic cluster of P. lentiscus. Indicated are only bootstrap values for the main groups, as obtained by neighbour-joining tree analysis.
Morphology
The results of the descriptive statistical analysis are shown for all the investigated Pistacia taxa in Tab. 3, and individually for each studied population in On-line Suppl. Tab. 1. Descriptive statistics for analysed leaf morphological traits for the inferred genetic groups are shown in Tab. 4. Based on the Kruskal–Wallis test, statistically significant differences ( P < 0.0001) among the studied taxa were confirmed for all studied morphological traits. However, multiple comparison analysis indicated that putative P. calcivora and P. terebinthus differ only in two traits (LW90 and NL), while a higher number of differences were found between the other pairs of taxa, particularly between the putative P. calcivora and P. lentiscus, which differed in all studied morphological traits (On-line Suppl. Tab. 2). This was further confirmed with descriptive statistical analysis, which pointed out the differences, as well as the similarities between the studied Pistacia taxa. In particular, P. terebinthus and P. calcivora exhibited a high degree of morphological similarity, as the two taxa displayed the largest leaflets and the longest petiole and rachis (Tab. 3). However, significantly lower NL was observed in P. calcivora (6.91), in comparison with P. terebinthus (7.75). Furthermore, P. lentiscus differed the most in its morphology. It had the lowest values in almost all of the studied morphological traits, with the average leaflet area of 1.44 cm2 and petiole length of 3.23 cm (Tab. 3). The only exception was NL, which was the highest among the studied taxa (7.94) in P. lentiscus. As expected, the presumed hybrid individuals between parental P. terebinthus and P. lentiscus, P. × saportae, showed intermediate values in all of the studied morphological traits, except for the NL (Tab. 3). However, when only those individuals confirmed as hybrid P. × saportae by genetic analysis were considered, they were of intermediate character for most traits (except for RL and NL) but were more similar to P. terebinthus than to P. lentiscus, when size of leaflets and rachis was considered (Fig. 1, Tab. 4). Furthermore, excluding the presumed hybrids of P. × saportae, putative P. calcivora stands out with the highest variability in leaf morphology (Tab. 3). In general, of the studied morphological traits, LA showed the highest variability, while parameters related to the leaflet shape (FC and LA10) displayed the lowest variability.
Tab. 3. Descriptive statistics for analysed leaf morphological traits of the studied Pistacia taxa. M – arithmetic mean, SD – standard deviation, CV – coefficient of variation (%). Analysed traits: LA – leaflet area; LL – leaflet length; MLW – maximal leaflet width; PMLW – leaflet length measured from the base to the point of maximum leaflet width; LW90 – leaflet width at 90% of leaflet length; FC – form coefficient; LA10 – angle enclosed by the main leaflet vein (the centre of leaflet blade) and the line connecting the leaflet base to a set point on the leaflet margin, at 10% of total leaflet length; RL – rachis length; PL – petiole length; NL – number of leaflets.
Tab. 4. Descriptive statistics for analysed leaf morphological traits for the three identified genetic groups of Pistacia taxa. M – arithmetic mean, SD – standard deviation, CV – coefficient of variation (%). Analysed traits: LA – leaflet area, LL – leaflet length, MLW – maximal leaflet width, PMLW – leaflet length measured from the base to the point of maximum leaflet width, LW90 – leaflet width at 90% of leaflet length, FC – form coefficient, LA10 – angle enclosed by the main leaflet vein (the centre of leaflet blade) and the line connecting the leaflet base to a set point on the leaflet margin, at 10% of total leaflet length, RL – rachis length, PL – petiole length, NL – number of leaflets.
At the population level, statistically significant differences ( P < 0.0001) were found among the studied populations. The locality Grohote (island of Šolta), containing populations P1, P2 and P3, seems to favour overall smaller but very variable leaves. On the other hand, the island of Korčula, with populations P4, P5 and P6, supported larger and less variable leaves, pronounced the most in P. × saportae (P5). Island of Vis, with populations P7, P8 and P9, displayed predominantly intermediate values in leaf morphology and variability, with the exception of P. terebinthus (P9) which had the largest leaflets at this locality. The locality Mt. Biokovo, containing three populations of P. calcivora and one population of P. lentiscus, was characterised with the largest (P12, P. calcivora) and smallest (P13, P. lentiscus) leaflets in this study.
PCA revealed that the first two PC axes had Eigen-values greater than one (>1). In addition, the cumulative variability explained by the first two PC axes was 84.3%, with the expected dominant percentage of 67.9% ascribed to the first PC axis (Fig. 3). The second PC axes contributed less to the overall variability, with 16.5%. The contribution of the individual measured traits to each of the PC axes can be found in On-line Suppl. Tab. 3. The first PC axis was highly negatively correlated with seven traits (MLW, LA, LW90, LL, PMLW, PL and RL), whereas the second PC axis was highly positively correlated with NL (Fig. 3). In the PCA biplot, a grouping of P. terebinthus and P. calcivora individuals, characterised by larger leaflet blades, is noticeable on the left side, unlike the right side, where P. lentiscus and presumed P. × saportae, with small leaflet blades, have grouped (Fig. 3). Identified hybrid individuals of P. × saportae, as confirmed by AFLP, grouped together with individuals of P. terebinthus and P. calcivora, whereas the remaining putative hybrid individuals showed partial overlap with P. lentiscus, although their intermediate character is clearly visible (Fig. 3).
Fig. 3. Principal component analysis (PCA) of the investigated Pistacia taxa (109 individuals) in the eastern Adriatic based on analysed leaf morphological traits. Each individual shrub/tree is indicated by a small sign, while the taxa barycentres are represented by larger ones. Individual assignment to the specific taxon is based on the initial a priori identification based on leaf morphology. Hybrid individuals corresponding to the confirmed P. × saportae are indicated with an asterisk.
Finally, DA for the three genetic groups of Pistacia individuals was performed. The first large group consisted of all P. terebinthus individuals, as well as all putative P. calcivora individuals. The second large group included P. lentiscus individuals, together with those individuals which were initially classified as P. × saportae in the field but were later identified as P. lentiscus by genetic analysis. The third group included only four hybrid individuals of P. × saportae confirmed by genetic analysis. DA revealed that traits LA10, LA and LL had the strongest discriminative power between the identified genetic groups (On-line Suppl. Tab. 4). Furthermore, classification analysis revealed that P. lentiscus individuals were correctly classified in 100% of cases, P. terebinthus individuals in 96.6% of cases, whereas only 50% of the confirmed hybrid individuals were classified correctly. Two confirmed hybrid individuals were correctly classified as P. × saportae, and two were classified as P. terebinthus.
Discussion
In our study, we aimed to resolve genetic relationships and taxonomy of four putative indigenous Pistacia taxa growing in the eastern Adriatic ( P. lentiscus, P. terebinthus, P. × saportae and P. calcivora) using molecular AFLP markers and morphology. Our results based on the genetic data have shown the presence of two well distinct genetic groups; the first group included all individuals of P. terebinthus and individuals initially identified based on morphology as an endemic taxon P. calcivora. The second group corresponded to P. lentiscus and included the majority of sampled individuals of the potential hybrid taxon P. × saportae. However, four individuals of P. × saportae showed a clear hybrid character at the molecular level and were placed between the two parental species in the phylogenetic tree (Fig. 2). Leaf morphology analyses have shown statistically significant differences between the identified genetic groups where LA10, LA and LL showed the highest discriminatory power. Moreover, we have demonstrated an overlap between individuals of P. terebinthus and putative P. calcivora based on the measured morphological traits, while putative individuals of P. × saportae exhibited intermediate morphological characteristics between the two parental species ( P. terebinthus and P. lentiscus) (Fig. 3). Nevertheless, confirmed P. × saportae hybrids were morphologically more similar to P. terebinthus. Overall, based on the combined obtained results, our study clearly confirms the presence of the hybrid taxon P. × saportae in the eastern Adriatic, while we found no support for the described endemic taxon P. calcivora.
Despite numerous previous studies investigating genetic relationships among different wild Pistacia species, particularly in the Mediterranean (Kafkas and Perl-Treves 2002, Golan-Goldhirsh et al. 2004, Kafkas 2006, Yi et al. 2008, Vendramin et al. 2009, Xie et al. 2014), to the best of our knowledge this is the first study of genetic relationships of indigenous Pistacia taxa in the eastern Adriatic. For comparison, Katsiotis et al. (2003) investigated the genetic relationships among selected Pistacia species and cultivars in Greece using RAPD and AFLP markers, including P. lentiscus and P. terebinthus from the native species in this region. However, the samples of these two taxa were collected only from the Botanical Garden in Athens and the study has primarily focused on introduced species and cultivars present in Greece. Moreover, the vast majority of studies investigating genetic variation and phylogeny in Pistacia have been focused on numerous cultivars and/or species of commercial interest involved in fruit and resin production, among other purposes (Shanjani et al. 2009, Ziya Motalebipour et al. 2016).
On average, we detected higher genetic diversity (HE and %P) in the P. terebinthus genetic group (including all individuals initially identified as putative P. calcivora), while lower genetic diversity was identified in the P. lentiscus group (Tab. 2). However, P. lentiscus exhibited higher frequency of rare alleles (DW). Although other more recent studies on the genetic diversity of wild Pistacia species mainly rely on SSR markers (Vendramin et al. 2009, Ziya Motalebipour et al. 2016, Karcι et al. 2020, Guney et al. 2021), making direct comparisons in terms of absolute values of genetic diversity parameters impossible, Ziya Motalebipour et al. (2016) and Karcι et al. (2020) also found higher genetic diversity in P. terebinthus compared to P. lentiscus, which is in relative terms consistent with our results. Higher genetic diversity of P. terebinthus, also in comparison to other species of the genus (Ziya Motalebipour et al. 2016, Karcι et al. 2020), may be due to the fact that P. terebinthus is likely one of the most recently evolved species of the Pistacia genus (Parfitt and Badenes 1997, Vendramin et al. 2009). On the other hand, P. lentiscus is considered the most distant species from the most ancestral P. vera (Parfitt and Badenes 1997, Kafkas 2006, Ziya Motalebipour et al. 2016, Karcι et al. 2020), which may explain higher genetic rarity observed for this taxon.
Pistaciacalcivora, a taxon presumed endemic to the eastern Adriatic, was initially described by Radić (1985: 40) from the Biokovo Mountain based only on morphological characteristics (On-line Suppl. Fig. 1). Moreover, in his work Radić described two varieties and 14 forms of P. calcivora based on morphology, indicating very high intraspecific morphological variability of this putative taxon, which may be influenced by specific microhabitat and/or soil conditions in the Biokovo mountain range known for its limestone bedrock. Pistacia calcivora was described as occurring on rocky limestone slopes and cliffs, being exposed to harsh environmental conditions and strong winds, occurring mainly in the transitional area between the eu-Mmediterranean and sub-Mediterranean zones along the Biokovo mountain range. Indeed, we have confirmed that the putative taxon P. calcivora has very high variability of leaf morphology compared to other investigated taxa and has significantly lower number of leaflets than P. terebinthus. Nevertheless, P. terebinthus and P. calcivora had overall very similar and overlapping leaf morphology (Fig. 3). Moreover, in the genetic clustering and phylogenetic analyses, all presumed P. calcivora individuals clustered in the P. terebinthus genetic group, confirming that taken together there is no difference at the molecular or morphological level between P. calcivora and P. terebinthus. Thus, they are the same species and therefore the name P. calcivora can be only considered as a synonym of P. terebinthus. We can only assume that the description of P. calcivora is a result of distinct morphological variability and genetic diversity within the P. terebinthus.
In addition to consideration of the putative taxon P. calcivora, the emphasis of our research was placed on the potential hybrid taxon P. × saportae in the eastern Adriatic. Pistacia × saportae has been recorded widely in the Mediterranean, however, to date it was not recorded for the area of the eastern Adriatic in the relevant global databases (Euro+Med 2006-2023, POWO 2023). Pistacia × saportae was officially recorded in Croatia for the first time on the island of Korčula (Jeričević et al. 2014) based on morphological characteristics, although we found an older record, specimens collected by Radić from 1985, in the herbarium MAKAR from the Biokovo area (On-line Suppl. Fig. 1). There are also few additional recent records in FCD of the putative P. × saportae along the eastern Adriatic on the islands Šolta, Vis and Krk (Nikolić 2023). In Spain, it was investigated and described by Werner et al. (2001) using RAPD molecular markers and morphological characteristics, which was also the first molecular study of this taxon. The authors confirmed the hybrid origin of P. × saportae and the proposed parental species. This was later further supported by a more detailed phylogenetic study of the whole genus by Yi et al. (2008) who confirmed the taxonomic status of P. × saportae using two nuclear (ITS and NIA-i3) and three plastid (cpDNA) regions (ndhF, trnC-trnD, and trnL-F). The latter study also identified P. terebinthus as paternal and P. lentiscus as maternal taxon. In our study, identified hybrid P. × saportae individuals had a higher proportion of P. lentiscus gene pool and positioned themselves closer to P. lentiscus in the phylogenetic tree (Fig. 2). However, based on our results we cannot rule out the gender of P. × saportae parental taxa.
There is not much available literature data regarding the morphology of the hybrid P. × saportae. AL-Saghir and Porter (2012) in their latest taxonomic revision of the genus describe the leaves of this species as being evergreen, up to 10 cm long, with seven to nine leaflets, which are on average 5.8 cm long and 2.2 cm wide. Their description was based on the examined material from Cyprus, Palestine, and Spain. Provided dimensions are slightly larger than those obtained in our study (4.05 × 1.66 cm). Werner et al. (2001) described morphological characteristics based only on material from the Iberian Peninsula, where they mention five to 11 leaflets. The authors already indicated that P. × saportae shows intermediate characteristics between the parental species, particularly regarding leaf shape, and displays high morphological heterogeneity. This pattern was largely confirmed in our study where putative P. × saportae individuals, identified a priori solely on morphology, were highly variable in their traits and were located between the parental taxa showing intermediate morphological characteristics (Tab. 3, Fig. 3). However, only four individuals (19%) were confirmed as P. × saportae hybrids at the molecular level, which is much lower compared to identification based on morphology. For example, the population from the island of Šolta, which was according to the initial identification based on leaf morphology identified as P. × saportae, was shown to be within the P. lentiscus genetic group (Fig. 2). On the other hand, four hybrid individuals of P. × saportae confirmed at the molecular level grouped together with individuals of P. terebinthus based on morphology (Fig. 3). Thus, our results emphasise previous observations that the correct identification of P. × saportae exclusively based on morphology is challenging and highly uncertain (Werner et al. 2001).
In conclusion, our results provide no support for treating the described taxon P. calcivora as a separate species; instead, it should be considered as a synonym of P. terebinthus. Furthermore, we have confirmed for the first time the presence of the hybrid taxon P. × saportae in the eastern Adriatic at the molecular level on the islands of Korčula and Vis. Given the wide natural distribution of the hybrid parental taxa P. terebinthus and P. lentiscus across the eastern Adriatic and eastern Mediterranean, we can expect P. × saportae to be also present in Montenegro, Albania, and Greece, although it may not been officially confirmed in those countries yet. Thus, its distribution along the eastern Adriatic should be further investigated.
Acknowledgements
This study was partially financed by the student project "Application of genetic markers for determining the taxonomic relationships of woody species (gentaxo)" financed by the Faculty of Forestry and Wood Technology, University of Zagreb.
We would like to thank Vanja Zorić and Lucija Gajić for their support during laboratory work. Thanks to our colleagues Mirjana Jeričević and Nebojša Jeričević who helped us in the field to collect Pistacia specimens on the island of Korčula and to Dinko Sule from Šolta.