Phenotypic, Genetic, and Epigenetic Variation among Diverse Sweet Cherry Gene Pools

Abstract

Sweet cherry germplasm contains a high variety of phenotypes which are associated with fruit size and shape as well as sugar content, etc. High phenotypic variation can be a result of genetic or epigenetic diversity that may interact through time. Recent studies have provided evidence that besides allelic variation, epiallelic variation can establish new heritable phenotypes. Herein we conducted a genetic and an epigenetic study (using amplified fragment length polymorphism (AFLP) and methylation-sensitive amplified polymorphism (MSAP) markers, respectively), accompanied by phenotypic traits correlation analysis in sweet cherry gene pools. The mean genetic diversity was greater than the epigenetic diversity (h_gen = 0.193; h_epi = 0.185), while no significant relationship was found between genetic and epigenetic distance according to a Mantel test. Furthermore, according to correlation analyses our results provided evidence that epigenetic diversity in predefined populations of sweet cherry had a stronger impact on phenotypic traits than their rich genetic diversity.

1. Introduction

Sweet cherry (Prunus avium L.) consists of phenotypes with high variability in fruit size and shape, sugar content, flowering time, pathogen resistance, and a variety of other traits [1,2,3,4,5,6]. Prunus avium cultivars originated from the Black and Caspian Seas, and the species currently presents a trans-European distribution as a result of the interplay between extant genetic variation and the genetic basis of adaptation [7]. The key to the evolution of sweet cherry breeding is therefore the in-depth study of genetic diversity and genes associated with traits of interest. In Greece, a large number of traditional local landraces are still highly varied, despite the fact that many are considered to have already been lost [7]. Moreover, many modern international varieties have a narrow genetic bottleneck, as reported by [8]. A variety of studies have been conducted in Greece using Simple Sequence Repeat (SSR) molecular markers and concern the analysis of the genetic basis of modern local cherry varieties. The main aim of the sweet cherry breeding program in Greece (‘Cherry4Breed’) was the development of new high-quality sweet cherry cultivars with high consumer appeal and self-compatible traits. Recently, the authors of [9], using whole genome resequencing (WGRS) analysis, found that a wild progenitor of Greek sweet cherry cultivars showed the highest genetic diversity. On the other hand, the mechanisms of epigenetic regulation in cherries have been poorly studied. To date only one such study has been published with regard to Greek wild cherry populations [10], investigating the possible coupling between epigenetic and genetic variation.

In eukaryotic organisms, gene regulation through genetic and epigenetic processes is of utmost importance for determining the phenotype [11,12,13,14]. Epigenetic changes such as DNA methylation [15] histone modifications [16], histone variants [17], and small RNAs [18] can present transgenerational inheritance. The ever-changing environmental conditions necessitate a more extensive study of epigenome variability [19,20,21]. Epigenetic regulation can be inherited, exerting influence on development, adaptation, and phenotype. Therefore, integrating the study of epigenetic variability and epigenome into breeding programs is of paramount importance.

Herein, in order to investigate epigenetic indices of diversity and structure as well as the relationships between epigenetic variants, genetic variation, and leaf traits, we used the amplified fragment length polymorphism (AFLP) and methylation-sensitive amplified polymorphism (MSAP) techniques along with correlation analysis statistics. Most of the leaf traits evaluated are very important for both sweet cherry breeding and cultivar registration in the test guidelines suggested by International Union for the Protection of New Varieties of Plants (UPOV). Our results may help breeders to understand and use epigenetic variation in order to accelerate the improvement of Prunus avium germplasm.

2. Materials and Methods

2.1. Plant Material

Twenty-two cultivated sweet cherry genotypes were selected (

Table 1. Leaf parameters for 24 sweet cherry accessions.

Cultivar	Predefined Population	Common Name	Surface Area	Caliper Length	Compactness	Convex Hull Area
BxS5	Breeding line	-	994,345	2.187.903	0.540	1,841,280
BxS33	Breeding line	-	456,843	1.264.413	0.543	841,147
HGxS11	Breeding line	-	1,293,434	2.314.375	0.614	2,105,830
BxS21	Breeding line	-	1,642,672	2.395.099	0.677	2,424,134
BxS14	Breeding line	-	1,421,852	2.214.082	0.646	2,200,543
BxS19	Breeding line	-	1,169,272	2.341.492	0.580	2,015,571
BxS17	Breeding line	-	1,163,390	2.030.177	0.610	1,904,195
PtrTrAch	Breeding line	Petrokeraso Tragano Achaias	897,861	1.965.073	0.580	1,545,661
HGxS30	Breeding line		990,953	1.929.465	0.679	1,457,352
ChalkAn	Landrace	Chalkidos Anonimo	1,211,613	2.042.699	0.590	2,051,944
Chi	Landrace	Chiou	959,994	1.941.165	0.546	1,756,523
Mz	Landrace	Mieza	1,334,036	2.177.894	0.630	2,116,484
PrKld	Landrace	Proimo Kolindrou	1,235,747	2.180.397	0.598	2,065,893
TrRd	Landrace	Tragana Rodohoriou	1,192,488	2.144.086	0.585	2,036,052
BxS22	Landrace	-	1,443,198	2347.44	0.706	2,042,288
AgLd	Landrace	Agiorgitiko Lilantiou	970,705	2012.02	0.602	1,611,516
Lmnd	Modern cultivar	Lemonidi	1,227,253	2.396.711	0.510	2,404,142
TrEds	Modern cultivar	Tragana Edes-sis	1,564,111	2.282.969	0.666	2,345,886
TrEdsN	Modern cultivar	Tragana Edessis-Naousis	1,166,668	2.049.908	0.634	1,839,019
Vas	Modern cultivar	Vasiliadi	1,595,101	2.245.907	0.710	2,245,480
Tsol	Modern cultivar	Tsolakeika	1,103,878	1.965.568	0.571	1,931,870
Bak	Modern cultivar	Bakirtzeika	1,184,177	2.363.933	0.611	1,935,432
Mhl	Wild	Prunus mahaleb	345,525	1.892.687	0.359	960,569
Wild	Wild	-	1,285,426	2.066.577	0.629	2,041,743

) from the Greek Fruit Gene Bank collection in Naousa (Institute of Plant Breeding and Genetic Resources- H.A.O. ELGO DEMETER), Greece, to represent the total diversity of Greek sweet cherry cultivars. The 2 elite wild cherry genotypes were obtained from the Wild Cherry Gene Bank, located at the Xyloupolis Forest Nursery of the Hellenic Forest Service in Greece. The 22 sweet cherry accessions are traditional Greek cultivars. We predefined 4 groups: ‘Breeding line’, ‘Landrace’, ‘Modern cultivar’, and ‘Wild’ (Table 1). Moreover, in the ‘Wild’ group we analyzed 1 genotype from Prunus mahaleb as a wild accession and also as outgroup in our study. Genomic DNA was isolated from a pool of young leaf samples of 5 plants per genotype using the NucleoSpin Plant II kit (Macherey-Nagel).

2.2. Leaf Phenotyping and Image Analysis

Plant images were obtained using a Scanalyzer PL semi-automated platform imaging system (LemnaTec GmbH, Aachem, Germany). The resolution of the digital camera (Baster AG, Ahresburg, Germany) was 1628 × 1236 pixels, with light in the visible spectrum in RGB (400–700 nm) with a pixel of 4.4 × 4.4 μm being used. Each plant was imaged individually in 3 (0°, 90°, and 180°) plane orientations. Plants were monitored over the mature leaf stage.

Plant pixel-projected leaf areas from the side view images were used to calculate the so-called ‘digital biomass’, which is an estimate of the plant volume and senescence. Data retrieved from the imaging platform were processed through an analysis Software (LemnaGrid) was specifically adjusted to obtain values for the phenotypic traits obtained from the images of each plant (Figure 1). The analysis was performed using a nearest-neighbor color classification, resulting in a binary image. Thereafter, geometric measurements of the object were taken, including caliper length (the longest dimension of the canopy when viewed from above), convex hull area (cm²), minimum enclosing area, and compactness calculated as the square of plant border length divided by the projected side or top area.

The caliper length, convex hull area, minimum enclosing area, and compactness of each population were assessed using the coefficient of variation (CV = 100 × standard deviation/mean).

2.3. AFLP Analysis

For the AFLP procedure, total genomic DNA (200 ng) was digested with 4 U of EcoRI and MseI for 3 h at 37 °C. Digested DNA fragments and EcoRI and MseI adapters were ligated with T4 DNA ligase (New England Biolabs) for 3 h at 26 °C. The resulting DNA was used as the primary template DNA in the AFLP analysis. A primer pair based on the sequences of the EcoRI and MseI adapters (

Table 2. Adapters and primers used for the amplified fragment length polymorphism (AFLP) and methylation-sensitive amplified polymorphism (MSAP) analysis.

	5′ to 3′ Sequence
EcoRI adapter	CTCGTAGACTGCGTACC AATTGGTACGCAGTC
MseI adapter	GACGATGAGTCCTGAG TACTCAGGACTCAT
HpaII/MspI adapter	GACGATGAGTCTCGAT CGATCGAGACTCAT
Pre-selective EcoRI primer	GACTGCGTACCAATTC-A
Pre-selective MseI primer	GATGAGTCCTGAGTAA-C
Pre-selective HpaII/MspI primer	ATGAGTCTCGATCGG-A
Selective EcoRI primers	GACTGCGTACCAATTC+ATG (FAM) GACTGCGTACCAATTC+ACT (HEX) GACTGCGTACCAATTC+AAC (ROX) GACTGCGTACCAATTC+AAG (TAMRA)
Selective MseI primer	GATGAGTCCTGAGTAA-CAA GATGAGTCCTGAGTAA-CAC GATGAGTCCTGAGTAA-CGT GATGAGTCCTGAGTAA-CTC
Selective HpaII/MspI primer	ATGAGTCTCGATCGGATC ATGAGTCTCGATCGGACT ATGAGTCTCGATCGGAAT

) with one additional selective nucleotide at the 3’ end (EcoRI + A and MseI + C) was used for the first PCR step (pre-amplification). Pre-amplification PCR was performed in a total volume of 20 μL containing 1X Kapa Taq Buffer, 0.2 mM of each dNTP, 2.5 mM MgCl₂, 30 ng of each primer EcoRI + A, MseI + C, 1 U Taq DNA polymerase (Kapa Biosystems), and 5 μL of diluted fragments (from the digestion and ligation reaction). Cycling was carried out on a BioRad thermocycler with a 95 °C hold for 30 s followed by 32 cycles of 95 °C for 30 s, 56 °C for 30 s and 72 °C for 1 min and subsequently followed by a final hold at 72 °C for 10 min. A 5-μL aliquot of the reaction was electrophoresed on agarose to verify amplification; the remaining 15 μL were diluted 5-fold with TE. Selective amplifications were carried out in 10-μL total volumes consisting of 3 μL of diluted pre-selective template and using the same reaction conditions as for pre-selective amplification but using 30 ng of an MseI primer and 5 ng of an EcoRI primer per reaction. Sixteen selective amplifications were performed on a BioRad thermocycler with the following program: an initial cycle of 95 °C for 30 s, 65 °C for 30 s, and 72 °C for 1 min; then 12 cycles of 95 °C for 30 s with an annealing temp starting at 65 °C for 30 s, but decreasing by 0.75 °C for each cycle, and then 72 °C for 1 min; and finally, 23 cycles of 94 °C for 30 s, 56 °C for 30 s, and 72 °C for 1 min, with a final hold at 72 °C for 30 min. Pre-selective and selective primers for the AFLP procedure are described in Table 2. Replicate analyses were conducted once, employing the same DNA extractions for AFLP and MSAP analysis.

2.4. MSAP Analysis

MSAP or methylation-sensitive amplified polymorphism is a modification of the standard AFLP technique. EcoRI is used as the rare cutter and the methylation-sensitive restriction enzymes HpaII and MspI as the frequent cutters due to the fact that the latter represent a pair of isoschizomers recognizing the same target sequence 5′- CCGG -3′, but they have differential sensitivity to methylation at the inner or outer cytosine. For the MSAP procedure, from each sample 200 ng of genomic DNA was digested with 4 U of EcoRI and 3 U of HpaII, and 200 ng was treated with 4 U of EcoRI and 3 U of MspI. The digestion was carried out at 37 °C for 3 h. The resulting DNA fragments and the EcoRI and HpaII/MspI adapters (Table 2) were ligated at 25 °C for 3 h using 400 U/μL of T4 DNA ligase (New England Biolabs). After that period the samples were subjected to a heat shock treatment for 10 min at 65 °C to end the ligation reaction. A primer pair based on the sequences of the EcoRI and HpaII/MspI adapters (Table 2) with one additional selective nucleotide at the 3′ end (EcoRI + A and HpaII/MspI + T) was used for the pre-selective PCR step. Pre-amplification PCR was performed in a total volume of 20 μL containing 1 × Kapa Taq Buffer, 0.4 U dNTPmix (10 mM), 2.5 mM MgCl₂, 30 ng of each primer EcoRI + A, HpaII/MspI + C, 1 U Taq DNA polymerase (Kapa Biosystems), and 5 μL of diluted fragments (from the digestion and ligation reaction). The cycling program was the same as for AFLP procedure. Twelve selective amplifications were performed following same procedure as for AFLP above. The pre-selective and selective primers used are given in Table 2.

2.5. Scoring of AFLP Markers

AFLP product mixtures were denatured in formamide at 94 °C for 2 min and electrophoretically separated using an ABI Prism 3730 xl automated fluorescence sequencer (Applied Biosystems). The size of detected fragments was determined by the Genemapper v4.0 program using an internal standard (GS 500 LIZ, Applied Biosystems). Fragments ranging from 150 to 500 bases in size were counted and further analyzed in order to reduce the impact of potential size homoplasy [22].

Analysis of the AFLP data was carried out by employing POPGENE version 1.32 [23] using the option for dominant diploid markers. GenAlEx 6.502b [24] was used to determine the genetic structure at hierarchical levels (analysis of molecular variance, AMOVA), to reveal band patterns and frequencies, and to plot grouping of individuals in the principal coordinates analysis (PCoA).

A cluster analysis using an unweighted pair-group method with arithmetic averaging (UPGMA) [25] was performed using the POPGENE 1.32 software [23]. Tests for statistical significance were based on 9999 random permutations, followed by sequential Bonferroni correction.

2.6. Scoring of MSAP Markers

Only reproducible fragments ranging from 150 to 500 bases were counted and further analyzed in order to reduce the impact of potential size homoplasy [22]. For MSAP analyses, comparison of the banding patterns of EcoRI/HpaII and EcoRI/MspI reactions resulted in 4 conditions of a particular fragment—I: fragments present in both profiles (1/1), indicating an unmethylated state; II: fragments present only in EcoRI/MspI profiles (0/1), indicating hemi- or fully methylated CG sites; III: fragments present only in EcoRI/HpaII profiles (1/0), indicating hemimethylated CHG sites; and IV: absence of fragments in both profiles (0/0), representing an uninformative state caused either by different types of methylation, or due to restriction site polymorphism [26]. To separate unmethylated and methylated fragments and to test for the particular impact of the methylated conditions II and III, we used the ‘Mixed-Scoring 2′ approach [27].

Epigenetic diversity within populations was quantified using the R script MSAP_calc.r [27] as: (1) number of total and private bands (polymorphic subepiloci), (2) percentage of polymorphic subepiloci (Pepi), and (3) mean Shannon’s information index (Iepi). GenAlEx 6.502b [24] was employed to compute haploid gene diversity (h) within populations. GenAlEx was also used to conduct analysis of molecular variance (AMOVA)—separately for each subepilocus class—in order to study the variation of CCGG methylation states (epiloci) among the 4 populations. Separate principal coordinates analyses (PCoA) were employed to assess differentiation among populations based on MSAP data according to different subepilocus classes (h, m, and u epiloci). The similarity between matrices based on different marker systems (MSAP and AFLP) was calculated using the standardized Mantel coefficient [28]. The same test was used to evaluate the similarity between geographic and genetic, and between geographic and epigenetic distances, respectively.

2.7. Correlations Between Inter-Population Genetic, Epigenetic and Phenotypic Distance

Correlation and regression analyses were used to assess the relationships between intra-population epigenetic and genetic diversity, between intra-population epigenetic diversity and phenotypic variation, and between intra-population genetic diversity and phenotypic variation. Particularly, in each case the linear and the polynomial models of second order (quadratic) were evaluated. In order to opt for the best fit, analysis of variance was used to compare the polynomial model of second order (complete) and the respective linear (nested) models, resulting in the model of choice. For each chosen model, 90% confidence intervals were also developed. Correlation and regression analyses were performed in R 3.6.3. The packages that have been used for analysis with R were “made4” [29], “corrplot” [30], and “ggplot2” [31].

3. Results

3.1. AFLP Genetic Diversity

Sixteen primer combinations revealed 542 markers, of which 51.34% were polymorphic (

Table 3. Intra-population (A) genetic and (B) epigenetic variation. PLP (%): percentage of polymorphic loci (PLP_gen) and epiloci (PLP_epi). h: Nei’s genetic diversity (h_gen) and epigenetic diversity (h_epi). I: Shannon’s information index of genetic diversity (I_gen) and epigenetic diversity (I_epi). CV: coefficient of variation.

	A. Genetic Variation			B. Epigenetic Variation
Pre-Defined Population	PLPgen (%)	hgen	Igen	PLPepi (%)	hepi	Iepi
Breeding line	67.53	0.218	0.335	54.06	0.174	0.268
Landrace	42.99	0.156	0.234	35.69	0.159	0.227
Modern cultivar	57.75	0.210	0.315	66.08	0.204	0.317
Wild	37.08	0.185	0.257	40.64	0.230	0.282
C. Phenotypic Variation
	CV Area	CV Caliper Length	CV Compactness	CV Convex Hull Area
Breeding line	30.427	16.681	8.598	26.147
Landrace	20.202	125.841	29.648	22.562
Modern cultivar	16.465	7.816	11.437	11.490
Wild	81.499	6.211	38.578	50.927

). The percentage of polymorphic loci ranged from 37.08 % for the ‘Wild’ predefined population to 67.53% for the ‘Breeding line’ population. Nei’s gene diversity (h_gen) ranged from 0.156 for the ‘Landrace’ population to 0.218 for the ‘Breeding line’ predefined population and Shannon’s diversity index (I_gen) ranged from 0.234 for ‘Landrace’ population to 0.335 for ‘Breeding line’ population (Table 3).

The analysis of molecular variance (AMOVA) showed that approximately 96% of the total genetic variation was partitioned within populations, while among population differentiation was 4% (

Table 4. Analysis of molecular variance (AMOVA) table for genetic (AFLP) and epigenetic (MSAP) data.

AFLP Summary AMOVA Table					MSAP Summary AMOVA Table
Source	df	SS	MS	Est. Var.	df	SS	MS	Est. Var.	%
Among populations	3	256,885	85,628	3776	3	142,444	47,481	2955	9%
Within populations	20	1,286,698	64,335	64,335	20	636,056	31,803	31,803	91%
Total	23	1,543,583		68,111	23	778,500		34,758	100%
Stat	Value	P (rand ≥ data)			Stat	Value	P (rand ≥ data)
PhiPT	0.055	0.037				0.085	0.001

). Furthermore, PCoA explained 28.89% of the total variation in two-dimensional multivariate space (Figure 2). Furthermore, the technical error rate from replicate analysis was 2% for AFLP replicate analysis.

3.2. MSAP Epigenetic Diversity

Fourteen primer combinations produced 283 markers, of which 102 were m-methylated (36%), 98 were h-methylated (35%), and 83 were u-methylated markers (unmethylated/uninformative; 29%). Among populations, the percentage of polymorphic markers (PLPepi) ranged from 35.69% for ‘Landrace’ to 66.08% for the ‘Modern cultivar’ population (Table 3). Nei’s diversity index (hepi) ranged from 0.159 for the ‘Breeding line’ population to 0.204 for the ‘Modern cultivar’ population, and Shannon’s diversity index (Iepi) ranged from 0.227 for the ‘Landrace’ population to 0.282 for the ‘Wild’ population (Table 3).

The analysis of molecular variance (AMOVA) showed that approximately 91% of the total epigenetic variation was partitioned within populations, while among population differentiation was 9% (Table 4). At the multivariate space, principal coordinate analysis of epigenetic distances revealed varying population differentiation patterns among all subepiloci (Figure 2) and explained 23.68% of the total variation in two-dimensional multivariate space. When a separate principal coordinates analysis was conducted for u-, m-, and h- epiloci, the absence of population differentiation was indicated. The first two axes respectively explained the percentages of the total variation: 28.46% for u-, 33.24% for m-, and 29.98% for h-subepiloci. Lastly, the technical error rate for MSAP replicate analysis was also low, presenting an average of 2.2%.

3.3. Correlations Between Intra-Population Genetic, Epigenetic and Phenotypic Variation

We found a strong concave-down relationship between intra-population genetic diversity and epigenetic diversity (Figure 3; for hepi: R² = 0.998, and for Iepi: R² = 0.975, Figure 3b,c), suggesting that in these cases epigenetic diversity assumed its highest values at the intermediate values of genetic diversity, and decreased at lower and higher values. The correlation for PLPepi was positive linear (R² = 0.547, Figure 3a), suggesting that epigenetic diversity was in concordance with genetic diversity.

Furthermore, we found strong correlations between intra-population genetic diversity and phenotypic variation for: (1) CV caliper length with h_gen (R² = 0.996, Figure 4e), (2) CV compactness with the percentage of polymorphic loci (PLP_gen) and hgen (R² = 0.997 and R² = 0.932, Figure 4g,h, respectively), and (3) CV convex hull area with PLP_gen (R² = 0.967, Figure 4j).

Assessing the relationship between intra-population epigenetic diversity and phenotypic variation, we found a strong concave-up correlation for: (1) CV caliper length with I_epi (R² = 0.998, Figure 5f) and (2) CV compactness with hepi (R² = 0.971, Figure 5h), as well as strong linear correlations for: (1) CV area with hepi (R² = 0.548, Figure 5b), (2) CV caliper length with hepi (R² = 0.551, Figure 5e), and (3) CV compactness with PLP_epi (R² = 0.678, Figure 5g), and weaker linear correlations for: (1) CV caliper length with PLP_epi (R² = 0.410, Figure 5d), (2) CV convex hull area with PLP_epi (R² = 0.359, Figure 5j), and (3) CV convex hull area with h_epi (R² = 0.315, Figure 5k).

3.4. Correlations Between Inter-Population Genetic, Epigenetic, and Phenotypic Distance

A Mantel test (Figure 6) revealed no correlation between inter-population phenotypic and epigenetic distance (R² = 9.10–7, P = 0.610). Similarly, inter-population genetic distance had no correlation neither with phenotypic differentiation (R² = 0.067, P = 0.3) or epigenetic differentiation (R² = 0.096, P = 0.02).

4. Discussion

Recent studies in many plant species report that epigenetic effects can be heritable and can affect plant growth, adaptation, and phenotypic plasticity [15,32]. Therefore, it is important to incorporate epigenetic variation into plant breeding programs in order to comprehensively employ the total available variation from the genome and epigenome [33,34,35].

In our study we used the predefined ‘Breeding line’, ‘Landrace’, ‘Modern cultivar’ and ‘Wild’ sweet cherry populations in order to assess the total diversity of Greek sweet cherry cultivars with the main aim of estimating correlations between genetic, epigenetic, and phenotypic traits.

Results showed that genetic diversity was high (mean value h_gen = 0.193) and the highest value was 0.218 for the ‘Breeding line’ population, which is a confirmation of high genetic diversity for sweet cherry cultivars. This was also reported in [9] after whole genome resequencing (WGRS) analysis of the same material. Epigenetic diversity (h_epi) presented a mean value of 0.185 for all populations, while the ‘Wild’ population had the maximum value of 0.282. This value is numerically higher compared to the previous wild populations which were studied in [10], where the respective mean value was 0.108. Furthermore, overall mean genetic diversity was higher compared to epigenetic diversity for Greek sweet cherry cultivars, a result which is in agreement with the Greek wild cherry populations studied with MSAP and ISSR markers in [5]. On the contrary, in a study of 96 accessions of Prunus mume [33], genetic diversity was reported as being lower than epigenetic diversity. Limited studies have been reported for congruent genetic and epigenetic analyses in other species. For example, Wang [34] reported higher epigenetic diversity for Hydrocotyle vulgaris L., which was introduced to China and is cultivated as a clonal herb.

4.1. Correlations Between Epigenetic and Genetic Variation

The relationship between genetic and epigenetic variation in populations is still in debate according to various studies [26,36]. The central issue is “whether epigenetic variation is completely uncoupled from genetic variation”, suggesting that population-specific selection could act on both genetic and epigenetic variation independently. If we accept this scenario, then we have to consider that this is probably a result of random epigenetic drift—dramatic spontaneous alterations of epigenetic marks occurring during mitotic propagation or epimutations [10,37].

The second scenario suggests that epigenetic variation is at least in part a downstream, subsidiary effect of genetic variation, and various studies support this thesis [26,38,39]. In our study a significant concave-down correlation was found between genetic and epigenetic gene diversity (h) and Shannon´s information index (R² = 0.998 and 0.975, respectively), and a linear correlation for PLP (R² = 0.547), indicating that epigenetic and genetic variations strongly interact. Herrera and Bazaga [39] also found that epigenetic differentiation for Viola cazorlensis populations was correlated with adaptive genetic divergence, which suggests that differences at the epigenome could arise by variable selection.

Furthermore, according to the Mantel test no significant relationship was found between the genetic and epigenetic distance of predefined populations of sweet cherry. A similar result of no concordance between genetic and epigenetic distances was also found for Greek Prunus avium wild populations in [10].

4.2. Correlations of Phenotypic Variation with Genetic and Epigenetic Variation

In order to understand and estimate the proportion of phenotypic diversity in predefined populations that is explained by genetic and epigenetic variation, we conducted a series of regression analyses. According to various genetic mapping studies, genetic variation can explain phenotypic variation, but the ‘missing heritability’ problem still exists [40,41]. The problem lies in the fact that not all of the genetic variance in the phenotype can be accounted for by variance in the genome [40]. On the other hand, numerous studies have shown that epimutations are a faster force of evolution; they can have a strong impact on adaptation and can provide higher phenotypic variance [42,43].

In the Greek sweet cherry gene pool, a high degree of phenotypic variation exists [4]. We herein report strong intra-population correlations between genetic and phenotypic traits. Other studies have also shown associations between AFLP genetic markers and phenotypic traits [44].

Interestingly, our results provide evidence that epigenetic variance generally contributes more to phenotypic variation than genetic variance. We found that epigenetic indices (h_epi, PLP_epi) were correlated with a higher number of phenotypic traits than genetic indices. Epigenetic indices either presented strong concave-up (two phenotypic traits), strong linear (three phenotypic traits), or weaker linear correlations (three phenotypic traits). Ma [33] also showed that leaf length, width, and area were positively correlated with methylation levels. In Zea mays cultivars, Roy [45] reported a significant amount of epigenetic variation between them. Moreover, the authors of [42] found that methylation changes rather than DNA sequence changes in a floral symmetry gene explained the phenotypic change in natural populations of Linaria Vulgaris.

5. Conclusions

Data acquisition, in the era of studying genetic/epigenetic variation and explaining phenotypic diversity in populations and cultivars, is of paramount importance. Herein, we showed that epigenetic variance in predefined populations of sweet cherry possibly has a stronger impact on phenotypic traits than their rich genetic diversity. This first insight could provide a significant baseline in order to support further research focused on targeting more epigenetic markers (through bisulfite sequencing) that could serve for breeding programs. This study has elucidated the notable role of epigenetic diversity in shaping phenotypic variance, and therefore it justifies the need to reveal possible epigenetic differentiation associated with different environmental conditions with greater resolution.