The study was carried out on twenty-seven autochthonous olive accessions (Olea europeae L.) collected from archeological sites and ancient olive groves localized in Central-Eastern Tunisia, which is actually subjected to a high risk of olive genetic erosion due to rapid and uncontrolled urbanization (Fig. 1). They were collected from old olives surviving as isolated trees or groups of trees. Those genotypes were selected for their different and interesting morphological characteristics compared to standard Tunisian cultivars. They were selected also based on trunk and stump size as the tree age indicator. They were picked respectively from Msaken–Boujriwil, Oued Laya, Sidi Bou Ali, Zaouiet Sousse, Hammam Sousse, Chott Mariem-Foukaia and Mahdia–Aliya (Table 1). All the olive groves studied were under rain-fed conditions. Five Tunisian cultivars: Zarrazi Zarsis (ZZ), Chemlali Zarzis (ChZ), Chemlali Sfax (ChS), Zalmati (Zal) and Chemlali Tataouine (ChT) were supplied from the national collection of olive in Boughrara–Sfax. They were taken as references (Table 2).

DNA was extracted from healthy young leaves using the CTAB method of Doyle and Doyle [22] with minor modifications. 0.5 g plant material were ground in liquid nitrogen and incubated in the extraction buffer (CTAB 2%, 20 mM EDTA, 1.4 M NaCl, 100 mM Tris–HCl pH: 8.0, 2% β-mercaptoethanol). The samples were incubated at 65 °C for 30 min, mixed with 750 μl of chloroform–isoamyl alcohol (24:1) and centrifuged at 10,000 rpm/min for 10 min. The aqueous phase was recovered and the DNA precipitated with 0.6 v of isopropanol and 0.08 v of sodium acetate (3 M). The pellet was washed in 70% ethanol, dried and re-suspended in 50 μL of TE buffer (10 mM Tris–HCl pH: 8, 0.1 mM EDTA). The determination of DNA quality and concentration was performed by a Nanodrop apparatus and electrophoresis in a 0.5% agarose gels.

Sixteen morphological and pomological traits were observed on each olive accession. All the measurements were evaluated for 40 samples of leaves and fruits per tree. After fruits characterization, stones were removed and subject of characterization. Leaf parameters were: leaf length (LL, mm), leaf width (LWI, mm), leaf surface (LS, mm²) and leaf length/width ratio (LR). The fruit parameters were: fruit weight (FW, g), fruit length (FL, mm), fruit width (FWI, mm) and fruit length/width ratio (FR). Stone parameters were: stone weight (SW, g), stone length (SL, mm), stone width (SWI, mm), stone length/width ratio (SR) and number of grooves (SG). Pomological traits were: flesh/stone ratio (FSR), oil content in percentage of wet matter (OWM %), and oil content in percentage of dry matter (ODM %).

From 23 RAPD primers tested, ten were selected for this study (Table 4). PCR was performed on an Eppendorf Mastercycler thermal cycler in a total volume of 20 μl. The PCR mix contained 5× PCR buffer, 8 mM dNTPs, 25 mM MgCl₂, 5 μM primer, 1 U Taq DNA polymerase, and 50 ng template DNA. The temperature profile consisted of an initial 2-min denaturation step at 95 °C, followed by 35 cycles at 94 °C for 1 min, primer annealing at 35 °C for 1 min and extension step at 72 °C for 2 min. The final elongation step was at 72 °C for 5 min. Amplification products were separated by electrophoresis through 1.6% modification products were separated by electrophoresis through 1.6% agarose gel stained with ethidium bromide in 0.5 × TAE buffer. The size of the fragments was estimated using the 1 Kb plus DNA ladder. The RAPD bands were visualized under UV light and photographed with a Photo Print (Vilber Lourmat, France) imaging system and analyzed by BioOne D++ software (Vilber Lourmat, France).

Nine SSR primer pairs were used (Table 1). The 20-μl PCR mix contained 10× PCR buffer, 2.5 mM dNTPs, 25 mM MgCl₂, 20 μM of forward and reverse primers, 0.5 U Taq DNA polymerase and 50 ng template DNA. The temperature profile included an initial denaturation step at 94 °C for 5 min, followed by 35 cycles of 94 °C for 20 s, primer annealing at 50 °C for 45 s and extension step 72 °C for 45 s. The final elongation step was at 72 °C for 8 min. Amplification products were separated by electrophoresis through 4% Agarose gel staining was followed by ethidium bromide staining in a 0.5× TAE buffer. The size of the fragments was estimated using 50 pb plus DNA ladder. The SSR bands were treated similarly to the RAPD bands.

The morphological and pomological traits were evaluated for their quantitative variables: mean, minimum value, maximum value, standard deviation (SD), and coefficient of variation (CV %). The collected data was analyzed by principal component analysis (PCA) using XLSTAT version 3.2016 software. Clustering analysis was also performed using NTSYSpc software version 2.1 [23] after standardization.

For RAPD analysis, primer efficiency was calculated by dividing the number of markers produced for each primer by the total number of markers obtained. The discrimination power was calculated by dividing the number of polymorphic markers amplified for each primer by the total number of polymorphic bands obtained [24].

For SSR analysis, the observed heterozygosity (H_o), the expected heterozygosity (H_e), the frequency of null alleles (R), the probability of identity (PI), and the polymorphism information content (PIC) were calculated using IDENTITY 1.0 software (Centre for Applied Genetics, University of Agricultural Sciences, Vienna, Austria). These parameters served to evaluate the information given by the microsatellite markers. The comparative level of polymorphism and informativeness of RAPD and SSR markers were estimated following the procedures previously described by Belaj et al. [25].

RAPD and SSR bands were scored as 1 (present) or 0 (absent) in a binary matrix for each primer. The Jaccard's similarity coefficient [26] was calculated by using the SIMQUAL module on each matrix. The similarity coefficients were then used to construct dendrograms, by using the unweighted pair group method with arithmetic averages (UPGMA) employing the SAHN. The cophenetic correlation coefficient was used to measure the goodness of fit. Pairwise comparison of distance matrices from independent molecular and morphological data was performed in NTSYSpc with the Mantel test [27] with 1000 permutations. Comparison of cophenetic matrices was done between matrices obtained with the same clustering algorithm. The results of the SSR markers were compared to the FAO Olive Germplasm Database (www.oleadb.it) [28].

The mean range, the maximum and minimum values, the standard deviation and the coefficient of variation (CV) for morphological traits are given in Table 3. The coefficient of variation varied from 11.04% (fruit ratio) to 80.25% (fruit weight). A wide degree of variability was obtained for most evaluated traits. The CV was high (>20%) for the fruit's weight and width (FW, FWI), stone weight (SW), flesh-to-stone ratio (FSR), leaf surface (LS), and oil content in percentage of wet and dry matter (OWM, ODM). All other morphological characters presented medium CV (20–15%) except fruit and stone ratios (FR, SR) and stone grooves (SG). The coefficients of variation recorded in this study were similar or even slightly higher than those previously reported in olive cultivars collection [29] and wild olives [30] using morphometric traits; however, the highest one was noted on fruit by the same authors. Fruit and stone weights varied from 0.56 g to 9.17 g and from 0.15 g to 1.23 g, respectively. The oil content in percentage of wet and dry matter varied from 7.2% to 30.70% and 15.1% to 52.1%, respectively. This can be explained by the fact that studied olive accessions included both olive oil and table olive.

Principal component analysis was performed for all morphological and pomological parameters. The first three components (PC1, PC2 and PC3) expressed respectively 42.4%, 21.6% and 12.8% of the total variance, accumulating 76.7% of the total observed variability. The first component was positively correlated with fruit size (FW, FWI, FL), stone size (SWI, SW), and oil content (OWM, ODM). In fact, a high and significant correlation was noted between fruit and stone weight (r = 0.902), length (r = 0.793), and width (r = 0.842). The correlations between quantitative traits among the fruit and stone dimensions were previously reported for wild olive trees [31], seedlings [32,33], and cultivated olive trees [30,31,34]. A significant correlation was found also between fruit weight and oil content (r = 0.618), which corroborates previous studies carried out in wild olive [30]. The second component was positively correlated with leaf size (LWI, LS) and negatively with flesh-to-stone ratio (FSR). The third component was positively correlated with fruit and stone shape (FR, SR). In fact, a significant correlation coefficient was noted between fruit and stone ratios (r = 0.730), as previously reported by Sorkheh et al. [35].

The PCA performed for morphological and pomological traits was useful for identifying the most important traits associated with the variation among the studied olive accessions (Fig. 2a). The most discriminating features were fruit size (FWI, FW, FL) followed by stone size (SWI, SW), oil content (OWM, ODW), leaf size (LWI, LS), stone shape (SR), flesh-to-stone ratio (FSR), and finally fruit shape (FR). Fig. 2b shows a projection of the olive accessions determined by the first two principal components. The accession (F1) from the region of Chott Mariem seems to diverge significantly from the other ones. This result is expected since it has presented the highest average value for weight (9.17 g), length (28.74 mm), and fruit diameter (23.20 mm). Furthermore, this accession presented a medium to low tree vigor that can be interesting for the expansion of intensification and mechanization of table olive trees. Accessions located on the left-hand side of the plan were distinguished by big fruit size (M17, ZB5, ZZ) and high flesh-to-stone ratio and oil content (M1, M2, M3, M6, M7, M8, M10, TW2). Theses accessions can be valued as a dual purpose olive tree; this means that it is simultaneously oil and table olive like the cultivar reference Zarrazi Zarzis (ZZ) considered as secondary cultivars in Tunisia. It is characterized by high oil content and good acidic composition, and it is usually cultivated in mixture with the variety Chemlali Sfax because of his self-incompatibility. Accessions located on the upper and right parts of the plane presented generally big leaf size (ZB4, HS1, SB3), low oil content and pulp percentage (TW3, ZB2, M4, ZB4, TW1, A, ZB1, SB1). Most of those accessions presented feral olive forms that can be issued from hybridization between wild and cultivated olive varieties as (TW3, A, SB1, HS1, ZB4). Accessions in the middle of the plane presented morphological and pomological characteristics similar to those of Chemlali Sfax, the most cultivated olive trees in Tunisia. It is a productive variety, self-fertile, drought tolerant and well adapted to the local tough environmental conditions. Other cultivars like Zalmati (Zal) and different Chemlali clones as Chemlali Tataouine and Chemlali Zarzis are also grouped in the center of the plan. All those cultivars resemble Chemlali Sfax in some morphological traits; however, they are genetically different, and each one of these cultivars has its own culture area and is differentiated by their own descriptive morphological card as described by Trigui and Msallem [11].

Ten decamer oligonucleotide primers were selected among the 23 ones that were screened (Table 4). Thirteen primers did not give scorable bands (OPA08, OPA12, OPAA03, OPAA16, OPE05, OPE06, OPC07, OPD07, OPD10, OPD15, OPR19, OPB01, and OPF03). Primers yielding valuable bands varied. A total of 92 bands were obtained from 32 olive accessions. The most bands were observed in OPA11 primer with 13 bands, whereas the fewest were observed in OPR13 and OPAA18 primers, with 6 bands. Eighty-six bands were polymorphic, with an average of 8.6 polymorphic markers generated per primer. The lower average has been reported by Brake et al. [36] when analyzing 13 Jordanian olive cultivars using 15 RAPD markers (5.73). The higher average has been reported by Martins-Lopes et al. [37] analyzing 38 Portuguese olive cultivars using 20 RAPD markers (13.1). This depends either on the number and discrimination power of the RAPD markers and the number and amount of diversity of plant material. Ten bands were unique, allowing the identification of nine olive accessions. The molecular weight of the amplification products varied from 200 bp (OPZ13) to 2000 bp (OPA11, OPE1). The primer OPAA18 amplified five polymorphic markers, showing the lowest discrimination power (5.8%), whereas OPA11 showed the highest discrimination power (15.1%). The polymorphism percentage varied from 83.3% to 100%, with an average of 93.4%. Previous investigation in RAPD markers reported similar polymorphism information using OPA9 [38,39], OPA10 [38–40], OPA11 [39], and OPE1 [41].

Microsatellites showed a high level of polymorphism throughout the analyzed accessions and revealed the existence of 63 alleles. The number of alleles per locus ranged from six (UDO99-39 and EMO-90) to nine (ssrOeUA-DCA18), with an average of seven alleles per locus, revealing a high level of variability (Table 5). Higher averages have been reported by other workers [42–44], most likely as a consequence of using more diversified olive genotypes or more polymorphic microsatellite markers. The range of allele sizes in base pairs varied from 105 bp at locus ssrOeUA-DCA17 to 255 bp at locus ssrOeUA-DCA3. The observed heterozygosity value (H_o), which represents the number of heterozygous individuals per locus, showed the lowest value at locus UDO99-39 (0.31) and the highest at locus ssrOeUA-DCA09 and ssrOeUA-DCA18 (0.94), with a mean value of 0.76. Expected heterozygosity values (H_e) ranged from 0.71 (ssrOeUA-DCA16) to 0.86 (ssrOeUA-DCA9) with an average of 0.81. Except for the ssrOeUA-DCA4, ssrOeUA-DCA17 loci and UDO99-39, H_e was lower than H_o. This observed heterozygosity deficiency may be related to the presence of null alleles, whose frequency values were positive. Baldoni et al. [45] and Trujillo et al. [46] reported a positive value of null allele frequency at loci ssrOeUA-DCA4. Muzzalupo et al. [47] and Belaj et al. [48] reported the occurrence of null alleles at locus UDO99-39. Positive null allele's frequencies have been also reported for the ssrOeUA-DCA17 locus by Baldoni et al. [45]. Polymorphic information content (PIC) showed an average value of 0.78. This index ranged from 0.66 (ssrOeUA-DCA16) to 0.85 (ssrOeUA-DCA9). All PIC values were higher than 0.5, which demonstrates the powerfulness of the markers used. This result is in accordance with the consensus list of SSRs for worldwide olive DNA typing recommended by Baldoni et al. [45]. Probability of identity (PI) ranged from 0.037 to 0.14 at DCA09 and DCA16, respectively. The accumulated probability of identical genotypes for all loci of 1.70 × 10⁻¹¹, which confirms the high discriminative level of the used microsatellite set. Similar PI was obtained by the analysis of 92 accessions of ten Iranian cultivars (1.2 × 10⁻⁹) analyzed by 13 SSR markers [49], 24 Spanish cultivars (1.70 × 10⁻¹¹) by 10 SSR [50] and 19 cultivars from Slovenian olive collection by 14 SSR markers (1.43 × 10⁻¹¹) [51].

The three markers gave different dendrograms, in which the 32 olive accessions were separated into different groups with different numbers of subgroups. In general, the grouping of the olive cultivars did not reflect a close relationship with their sampling sites. However, a certain tendency to group the olive trees according to their geographical distribution was observed in the dendrograms issued from molecular data.

The UPGMA cluster analysis of the accessions based on the morphological characters revealed three major groups of accessions, mainly according to fruit size. This result corroborates other studies carried out in other classic cultivars based on morphometric characters [29,52,53]. The first group included six olive accessions (TW2, ZB5, M7, M8, M3 and M6) with high fruit weight and flesh to stone ratio and important oil content. The second group enclosed 13 accessions with medium fruit size and oil content. The rest of the accessions clustered together in Group 3. They were characterized by medium to small fruit size, low oil content and relatively big leaf size. This last group included four accessions which are feral olive trees (SB1, SB2, A, TW3).

The UPGMA cluster analysis of the accessions based on RAPD data indicated clear differences among olive accessions, with the similarity coefficients between all possible pairs of genotypes ranging from 0.25 to 0.88 (Fig. 3a). Eight main clusters were defined. Three accessions picked respectively from Aliya–Mahdia, Foukaia–Chott Mariem, and Hammem Sousse, which are geographically distant, were separated and formed individual branches. The five olive references cultivars were clustered together. The five accessions picked from Zaouiet Sousse were grouped also together. Accessions picket from Msaken–Boujriwil were grouped in a same cluster, except M8 and M9. These results agree with those of previous studies performed using RAPD markers, in which a clustering of olive cultivars from different geographic origins had been observed [7,54]. No apparent clustering according to morphological traits or fruit use of olive accessions was found contrary to other study with RAPD markers [55–57]. However, Belaj et al. [40], Sanz-Costés et al. [54], Martins-Lopes et al. [37] and Parra-Lobato et al. [7] found no apparent clustering of olive cultivars by fruit size or other morphological traits.

The UPGMA cluster analysis of the olive accessions based on SSR data demonstrated the presence of seven major groups. The similarity coefficient varied from 0.21 to 1. The lowest similarity was noted on the accession Mahdiya-Aliya, which was separated in the first cluster. This accession is for millennia characterized by multiple trunks, low oil content although interesting oil composition. The second cluster contained the accessions picked from Hammem Sousse, which were characterized by a vigorous tree and considerable trunk indicating an old age, likely corresponding to more than a hundred years, and an important flowering and pollen production (data not shown). This accession can be valued as pollinator (or Dhokar as genetic Tunisian nomination). The third cluster contained accessions picked from Foukaia–Chott Mariem and one accession from Zaouiet Sousse. The last two accessions were characterized by a very big fruit size. Rekik et al. [16] noted that two Tunisian olive table Besbessi and Tounsi were well separated from the other cultivars, which was probably the result of a different genetic background, as might be expected from their very big fruit size. The fourth cluster grouped the olive cultivars taken as references. The rest of the accessions were divided into three clusters. The highest similarity coefficients were noted respectively between M2 and M3, SB3 and SB6, SB1 and SB2, SB4 and SB5, ZB3, and ZB4. All pairs were originated from same geographic origin. Similar clustering tendencies were reported according to the genotype's geographical origin in previous study using SSR markers [33,58,59]. The clustering of the accessions from the same or nearby region suggests a common genetic base and an autochthonous origin for these accessions, which agreed with the hypotheses of autochthonous origin as well as the limited diffusion of olive cultivars from their origin zones of cultivation [40,60]. No clear clustering according to fruit size or olive use was noted, as reported by Delgado et al. [61] and Trujillo et al. [46].

The classic Tunisian olive cultivars taken as references were clustered together on the two dendrograms issued from molecular markers (RAPD, SSR). Present clustering is in overall accordance with previous studies performed using molecular analysis by SSR markers carried out on a large number of Tunisian olive cultivars [16,18,62]. Ben Ayed et al. [62], studying the genetic relationship of 15 cultivars corresponding to the major and widely distributed Tunisian cultivars, noted that Chemlali Tataouine differed from Chemlali Sfax, which suggests that nomination reflects just similarity on fruit and leaves. Rekik et al. [16] found that the maximum genetic similarity was found between Zalmati and Chemlali Sfax. Compared to the allele size reported for seven out of the 14 SSR markers analyzed in the international database of olive cultivars (www.oleadb.it), the two varieties Chemlali Sfax and Zalmati showed a good identity and minor variations with the material of the database. This result confirms the large clonal variability of those two varieties. Moreover, a comparison was done for all studied genotypes with the material of the database to identify synonymy for our genotype. Any synonymy with the cultivars recorded was found. However, the allele sizes were tentatively harmonized according to the list of allele sizes displayed for each marker.

The cophenetic correlation coefficient value between the dendrogram and the original distance matrix estimated from the morphological, RAPD and SSR markers was r_morph = 0.78, r_RAPD = 0.86 and r_SSR = 0.89, respectively. These cophenetic correlation coefficients indicated a good fit for molecular data and modest fit for morphological data according to the phenograms (Fig. 4).

The Mantel matrix correspondence test was used to compare the similarity matrices. No correlation between either the morphological and RAPD similarity matrices (r = −0.04) or the morphological and SSR similarity matrices (r = −0.10) was found. In fact, no clustering similarities were observed between morphological and molecular data phenograms. Similar estimates of the correlation coefficient between morphological and molecular data were obtained in previous study [30,35,55,63]. This could be due to the limited number of characters that were selected and measured. The discordance between morphological and molecular diversity may result from similar selection pressure, which leads to similar forms with different genetic backgrounds [63]. Owing to the very scarce genomic information that is available for the olive, it is very difficult to predict how many genes are involved in the morphological analyzed characters. Qualitative traits are expected to be under monogenic or oligogenic control [64]. Conversely, quantitative traits have a more complex genetic base, as they are governed by multiple genes and their interactions [64].

The correlation coefficient between matrices based on RAPD and SSR markers was significant but relatively low (r = 0.45) (Fig. 5). Belaj et al. [25] reported significant coefficient of correlation between RAPD and SSR data (r = 0.39) using eight SSR markers and 21 RAPD markers on 32 olive cultivars. A high and significant correlation was noted also between the consensus molecular data (RAPD + SSR) and both SSR data (r = 0.76) and RAPD data (r = 0.90). This can be explained by the large number of RAPD bands compared to that of SSR detected bands.

However, in spite of the relationships and concordance between the three used markers, molecular markers (RAPD-SSR) should not be seen as substitutes for the traditional morphological descriptors. All these marker systems should be considered as complementary tools for providing a more complete understanding of the diversity available in olive genotypes. The lack of consistency between the relationship studies performed with molecular and morphological markers could indicate that each marker system measures different aspects of olive genetic variability [30].

All studied molecular markers proved to be highly effective in discriminating the 32 accessions analyzed. The results obtained are summarized in Table 6. Different levels of expected heterozygosity (H_e) were observed for RAPD and SSR, which is consistent with results from previous studies carried out on olive cultivars by means of different molecular markers [25,57]. The co-dominant nature of the SSR marker permits the detection of a high number of alleles per locus and contributes to higher levels of expected heterozygosity being reached than would be possible with RAPD [25]. For the SSR markers, the effective number of alleles per locus (Ne) was higher than that noted for RAPD markers. This was reflected in lower values of the expected heterozygosity for both RAPD and AFLP markers. The very low value of the effective number of alleles per locus for SSR markers in comparison to the average number of alleles per locus may be explained by the presence of many unique or low frequent alleles [25].