Epigenetic Variability in the Genetically Uniform Forest Tree Species Pinus pinea L

There is an increasing interest in understanding the role of epigenetic variability in forest species and how it may contribute to their rapid adaptation to changing environments. In this study we have conducted a genome-wide analysis of cytosine methylation pattern in Pinus pinea, a species characterized by very low levels of genetic variation and a remarkable degree of phenotypic plasticity. DNA methylation profiles of different vegetatively propagated trees from representative natural Spanish populations of P. pinea were analyzed with the Methylation Sensitive Amplified Polymorphism (MSAP) technique. A high degree of cytosine methylation was detected (64.36% of all scored DNA fragments). Furthermore, high levels of epigenetic variation were observed among the studied individuals. This high epigenetic variation found in P. pinea contrasted with the lack of genetic variation based on Amplified Fragment Length Polymorphism (AFLP) data. In this manner, variable epigenetic markers clearly discriminate individuals and differentiates two well represented populations while the lack of genetic variation revealed with the AFLP markers fail to differentiate at both, individual or population levels. In addition, the use of different replicated trees allowed identifying common polymorphic methylation sensitive MSAP markers among replicates of a given propagated tree. This set of MSAPs allowed discrimination of the 70% of the analyzed trees.


Introduction
DNA Cytosine methylation has been shown to play a determinant role in a variety of molecular processes such as regulation of plant gene expression during development [1], imprinting [2] or genome stability including mobile elements control [3,4] and polyploidization events [5,6].
These functions have important implications not only in fields like developmental biology [1] but also in ecology and evolution.Epigenetic mechanisms have been proposed to contribute to adaptation in plants [7][8][9].Several recent studies have identified correlations between epigenetic variability and adaptive population differentiation of plants in response to environmental stresses such as drought [10,11], salinity [12,13], or damage by herbivores [14,15].
Environmentally-induced epigenetic changes have been shown to mediate phenotypic plasticity by regulation of specific gene expression as well as plant development after a change in environmental conditions [16][17][18].It is also known that epigenetic variability can be independent from genetic variability [19][20][21][22], becoming a source for adaptive potential in itself [18,[23][24][25].Epigenetic changes induced by stress are potentially reversible but some modifications are not only inherited from cell to cell during mitosis but they can also be inherited across generations [26][27][28][29].This so-called ''stress memory'' allows plants to retain active molecular mechanisms after the stress signal disappears, thus responding more efficiently to recurrent stressful conditions [16][17][18].Stress memory can considerably increase the adaptive potential and may help plants to cope with changing environmental conditions [24].
Although the number of studies about epigenetic variation associated with biotic and abiotic stresses in plants is increasing, few studies are focused on forest tree species, with perhaps the exception of poplar [10,11,30,31].Trees, and especially conifers, are key models for the study of stress adaptation due to their longevity and long-life cycles [9,32].Some conifers like Sequoia sempervirens or Pinus longaeva can live for 3,000 and 5,000 years, respectively [33].Therefore, these species must cope with very variable environments through their life spans.Conifer genomes are very large, with genome sizes ranging between 20,000 and 30,000 Mbp, which on average is 200-fold larger than Arabidopsis genome and 6-10-fold larger than the human genome [34,35].Thus, regulatory machinery for gene expression and genome stability must be a key factor for these species in order to survive under changing environmental conditions.From an ecological and economic standpoint, conifers are the most important group of gymnosperms.Altogether, they represent 39% of the world's forest area [36].However, little research about epigenetics, and more specifically cytosine methylation, has been done in conifers.The most explored fields in conifers have focused on epigenetic processes in tree development [37][38][39] and epigenetic memory to environmental factors [40][41][42].
Pinus pinea L. (stone pine) is one of the most ecologically, economically and socially important Mediterranean forest tree species.It is patchily distributed in the North and Southeast Mediterranean area, from Portugal to Syria.Stone pine is characterized by a very low genetic variation [43][44][45] and high adaptive plasticity that increases its global fitness [46][47][48][49].High degree of phenotypic plasticity has been found in response to water availability.The analysis of propagated trees grown under water deficit revealed a significant variation in functional traits [50].This genetically depauperated but plastic species constitutes an optimal system to study natural epigenetic variability and its potential to shape phenotypic plasticity [24].
The main goal of this work was to analyze cytosine methylation in Pinus pinea genome.Despite the lack of genetic variation of this species we expect to identify methylation variability between individuals that might explain the significant variation in functional traits observed in the species.Two different objectives were outlined: 1) Analyze if P. pinea genome is methylated and the extent of methylation.2) Analyze if cytosine methylation is correlated with genetic variability or if cytosine methylation patterns differ among and within individuals.To carry out this study we have analyzed DNA from vegetatively propagated individuals from natural populations of stone pine using two genome wide profiling techniques, Amplified Fragment Length Polymorphism (AFLP) and Methlylation Sensitive Amplified Polymorphism (MSAP) surveying both genetic and epigenetic variability, respectively.

Plant Material
A total of 20 one-year-old seed-grown individuals from five natural populations representing the distribution of Pinus pinea L. in Spain (Tordesillas, Bogarra, Biar, Don ˜ana and Palafrugell; Table S1 and Figure S1) were selected for this study.The two most represented populations, Tordesillas and Bogarra, have contrasting climates; Tordesillas has a colder continental climate while Bogarra has a temperate Mediterranean climate.Vegetative propagation of these individuals was conducted by planting cuttings in a mix of equal amounts of peat and river sand using 1% IBA (Rhizopon AA powder) to promote rooting.A set of 95 rooted cuttings (ramets) were obtained.After two months, the ramets were transplanted into 1.2 l containers with a 3:1 mixture of peat and river sand.Plants were grown in a climatic chamber under controlled conditions (photosynthetic photon flux density (PPFD) of 600-650 mmol m 22 s 21 , temperature of 20uC, relative humidity of 60% and photoperiod of 16/8) placed in a random block design consisting in four blocks with 1-2 ramets of each propagated tree per block.Four months later, needles of similar developmental stage and 2 cm below tip of main apical shoots were collected from every ramet and stored at 280uC for subsequent DNA extraction.Elongating shoots with very young needles were discarded during sampling as well as needles from initial rooted cuttings, which originated from the mother tree.In this respect, needles of same ontogenic state were carefully selected to reduce methylation variation associated to different developmental stages.

AFLP analysis
A total of 59 ramets from the 13 propagated trees belonging to the two most represented Spanish populations were analyzed using Amplified Fragment Length Polymorphism (AFLP) [53].This analysis was performed by digesting 500 ng DNA with EcoRI/ MseI restriction enzymes according to Cervera et al. [54].The number of selective nucleotides for the two consecutive amplification steps was EcoRI + 1/MseI +1 for the pre-amplification and EcoRI +3/MseI +3 for selective amplifications.Two primer combinations (Table S2) were used: EcoRI + ACC/MseI + CCA and EcoRI + ACG/MseI + CCA.EcoRI +3 selective primers were labeled at their 59 ends with fluorescence dye 800 IRDye to allow visualization of the fragments on a Li-Cor 4300 DNA Analyzer (Li-Cor Biosciences, Lincoln, NE).
Scoring of the resulting fragment patterns was based on a presence/absence (1/0) approach.Only markers with an undoubtedly reliable score of at least 95% of the samples (less than 5% of missing data) were considered to estimate genetic variability.

MSAP analysis
Methylation Sensitive Amplified Polymorphism technique [55], modified by Cervera et al. [20], was used to analyze DNA cytosine methylation level and pattern in the genome of P. pinea propagated trees.
Due to the large genome size of conifers [35] several steps of the technique were optimized for pine species.The initial amount of DNA digested with the two restriction enzyme combinations (EcoRI/HpaII and EcoRI/MspI), was increased up to 500 ng [52].A combination of EcoRI +1//HpaII/MspI +1 selective nucleotides in the pre-amplification followed by EcoRI +3// HpaII/MspI +3 selective nucleotides in the selective amplification step provided the best results.EcoRI +3 selective primers were labeled with fluorescent dye as in AFLP.
Initially, nine different EcoRI +3//HpaII/MspI +3 primer combinations (AAC/AAT, ACA/AAT, ACT/AAT, ATC/AAT, AAC/ACT, ACA/ACT, ACG/ACT, ATC/ACT, AAC/ATC; Table S2) were tested on a sample subset to identify the most informative combinations.For this purpose, each propagated tree was represented by a pool made of equimolar amounts of preamplified DNAs from its corresponding ramets.All pools were analyzed using the nine primer combinations to compare their MSAP profiles, selecting the two most informative ones: EcoRI + AAC//HpaII/MspI + AAT (AAC/AAT) and EcoRI + ACA// HpaII/MspI + AAT (ACA/AAT).These two primer combinations were used to analyze cytosine methylation patterns of the 95 ramets individually.Electrophoresis settings were similar to those applied for AFLP analysis.

Scoring and interpretation of MSAP fragment patterns
Comparative analysis between EcoRI/HpaII and EcoRI/MspI profiles for each primer combination allows establishing the methylation status of each targeted restriction site.Methylationsensitive endonucleases HpaII and MspI cleave CCGG sequences with differential sensitivity to methylation at the inner or outer cytosine: HpaII does not cut if one or both cytosines are fullmethylated (methylation occurs in both DNA strands) but cleaves when cytosine methylation occurs in a single strand.MspI does not cut if the outer cytosine is methylated in one or both strands [56,57].
Initially, separated matrices were constructed for EcoRI/HpaII and EcoRI/MspI fingerprints.MSAP fragment presence or absence was visually determined by two independent observations.We detected fragments differing in intensity probably due to different degree of cytosine methylation in different cell types of the analyzed samples.Only markers with an undoubtedly reliable score of at least 95% of the samples (less than 5% of missing data) were considered to estimate epigenetic variability.Rationale for the comparative scoring was based on differential presence/ absence of a particular fragment in HpaII and MspI digestions.Thus, for a given sample, hypomethylation (fragment present in EcoRI/HpaII and EcoRI/MspI fingerprints) and full methylation of both cytosines (fragment absent in EcoRI/HpaII and EcoRI/ MspI fingerprints relative to other samples) were coded as 0. Since P. pinea shows undetectable levels of genetic variation, the loss of the target sequence motif cannot be considered within this class.On the other hand, for a given sample, full methylation of the internal cytosine (fragment only present in EcoRI/MspI fingerprint) and hemi-methylation of the external cytosine (fragment only present in EcoRI/HpaII fingerprint) were codified as 1.The resulting integrated matrix was used for statistical analysis.
MSAP markers were then classified according to their global pattern in all samples (Figure 1-B).Two main groups were identified depending on whether there was at least a difference between EcoRI/HpaII and EcoRI/MspI digestions profiles.Markers were then identified as Methylation Insensitive (MI) when no difference was found between the two digestions profiles for any sample and as Methylation Sensitive (MS) when difference between both profiles was found for one or more samples.These two groups were further split according to whether difference among samples was found.Following this reasoning, MI markers presenting the same pattern among all samples were classified as Monomorphic Methylation Insensitive (MMI) markers and MIs presenting differences among samples were classified as Polymorphic Methylation Insensitive (PMI) fragments.MS fragments were classified as well into Monomorphic Methylation Sensitive (MMS), when they showed different pattern between isoschizomers but not among samples, and Polymorphic Methylation Sensitive (PMS) fragments when at least one sample did not show the same profile.

Statistical analysis
Percentages of cytosine methylation were subjected to analysis of variance (ANOVA; Statistica [58]) to unveil differences in the degree of cytosine methylation among propagated trees.MSAP markers showing the same profile among all individuals derived from the same propagated tree were identified.To analyze its discriminative power, epigenetic similarity (ES) was estimated from the number of shared amplified fragments by using the Dice similarity coefficient [59] [ES(ij) = 2a/(2a+b+c)] where ES(ij) is the measure of ES between the individuals i and j, a is the number of polymorphic fragments that are shared by i and j, b is the number of fragments present in i and absent in j, and c is the number of fragments present in j and absent in i.The resultant matrix was subjected to cluster analysis by the unweighted pair-group method analysis (UPGMA) and a dendrogram was constructed according to the clustering.Clustering was subjected to bootstrapping in order to obtain values for the reliability of the consensus dendrogram.Similarity matrix was obtained using DistAFLP software [60].Using Bootstrap Computation, 1000 matrices were obtained.Cluster analysis and dendrogram construction were performed with PHYLIP phylogeny software package (programs Neighbor and Consense, respectively) [61].Dendrogram was visualized with MEGA software [62].
Analysis of the Molecular Variance (AMOVA) [63] based on polymorphic methylation sensitive markers was performed over the 59 ramets of the two most represented populations, Tordesillas and Bogarra (Arlequin, version 3.5 [64]).Locus by locus AMOVA was performed to identify markers with a significant effect on population differentiation or differentiation of propagated trees (Table S3).These markers were then used to perform a Principal Component Analysis (PCA; Statistica [58]).

Genetic variability in Pinus pinea
A total of 59 ramets from 13 propagated trees of the two most represented Spanish populations (Tordesillas and Bogarra) were analyzed using Amplified Fragment Length Polymorphism (AFLP).A total of 215 AFLPs were identified with confident reliability using two primer combinations (EcoRI + ACC/MseI + CCA and EcoRI + ACG/MseI + CCA).A single AFLP fragment pattern was observed and no variation was found among ramets from each propagated tree as well as among different propagated trees.

Epigenetic variability in Pinus pinea
DNA methylation variability among the 95 ramets from the 20 propagated individuals was analyzed comparing MSAP profiles.The two selected MSAP primer combinations yielded a total of 216 scored markers (Table S4) from which 139 were classified as MS and 77 as MI.Within MS markers, 91 were identified as PMS (42.13% of the total number of MSAPs).The remaining 48 MS markers were identified as MMS MSAPs.Out of the 77 MI markers, 66 were found to be MMI MSAPs.The remaining 11 MSAPs (5.09% of the total number of MSAPs) were identified as PMI.Ten out of these 11 PMI MSAPs showed a different pattern in at least one ramet of the propagated trees.Detailed classification per primer combination is shown in Figure 1-A.
The EcoRI + AAC//HpaII/MspI + AAT (AAC/AAT) primer combination was the most informative with 119 out of the 216 amplified MSAPs analyzed.The main difference between primer combinations was found in the number of PMS markers since MMS, MMI and PMI markers showed similar values for the two primer combinations (Figure 1).
Comparison of EcoRI/HpaII and EcoRI/MspI profiles showed contrasting levels of polymorphism.While EcoRI/MspI provided a higher number of MSAPs than EcoRI/HpaII, 116 versus 91, their fragment patterns were less polymorphic.In particular, 82 out of the 91 PMS markers showed variation only in the EcoRI/ HpaII profiles, 5 only in the EcoRI/MspI profiles and the remaining 4 were associated with polymorphic MSAPs identified in both profiles (Figure 2).Thus, 91.67% of the MMS fragments were only found in EcoRI/MspI profiles (Figure 2).We explored differences in methylation level among the propagated trees estimating the number of MS markers detected for each ramet vs. the total number of MSAPs.The resulting mean value and standard deviation of all ramets corresponding to the same propagated tree were calculated.Values ranged from 42.7360.88%(Pal 27) to 47.9060.42%(Tor 27) (Table 1).DNA methylation significantly varied among the 20 different propagated trees (ANOVA, p,0.0001).Cytosine methylation polymorphism among ramets of each propagated tree, based on PMS detected among ramets vs. total number of MSAP fragments, was also calculated with values ranging between 0.46% (Tor 3) and 9.72% (Don 13).
Similarity analysis among MSAP profiles of the propagated trees allowed the identification of 15 PMS that, while being polymorphic among the analyzed trees, shown the same pattern among their vegetatively propagated ramets (Table S4).These markers were used to calculate an epigenetic similarity matrix based on DICE coefficient and to perform an UPGMA cluster analysis.As a result, the use of the 15 PMS MSAPs allowed to identify 14 out of the 20 studied individual trees (Figure 3) meaning that for these trees, all their ramets clustered together in common branches.Bootstrap values for these clusters were above 50% in 11 of the 14, and above 25% in all cases.All propagated trees from both Bogarra and Don ˜ana populations were clustered as well as 5 out of 8 trees from Tordesillas.Trees from both Bihar and Palafrugell populations were not clustered together.
Although the number of analyzed trees is scarce to approach population epigenetic studies, we roughly estimated variability of DNA cytosine methylation associated with the studied populations.For this purpose, the mean value of the methylation levels obtained for all propagated trees of a given population was calculated.Palafrugell was the population whose individuals showed the lowest level of DNA methylation with 43.8961.64%methylated cytosines, followed by Tordesillas, Biar, Don ˜ana and Bogarra with 44.8961.32%,45.6662.06%,45.7060.44%and 45.7861.35%,respectively.We also found variation in the percentage of PMS MSAPs among populations (for a given population, PMS MSAPs vs total MSAPs).The less polymorphic population was Palafrugell (13.89%) and the most polymorphic one was Tordesillas (27.31%).Biar, Bogarra and Don ˜ana showed intermediate values of 16.67%, 17.13% and 18.52%, respectively.
We tested the discriminative power of MSAP markers for the same two populations analyzed with AFLP technique by performing an Analysis of the Molecular Variance (AMOVA).The fixation index between populations (F ST ) was 0.274 (p, 0.0001).A locus-by-locus AMOVA was performed to determine which markers showed statistically significant epigenetic variation among populations.The resulting 52 markers (69.3% of the total polymorphic MSAPs identified in the two populations; Table S3) were used to differentiate populations and propagated trees using a Principal Components Analysis.First two components accounted for 35.24% of the total variance (comp. 1 = 21.60%;comp. 2 = 13.64%).Scores for these two components of each analyzed ramet were plotted in Figure 4. First component clearly differentiated both populations.In addition, it was possible to identify three propagated trees from Tordesillas population (Tor 3, Tor 7 and Tor 25) whose ramets clustered in separate groups.

Discussion
Stone pine has been described as a genetically depauperated species, showing a very low level of genetic diversity [43][44][45].The almost undetectable genetic variation has made it very difficult to genetically distinguish stone pine trees.In this study, we were unable to identify a single polymorphic marker using the multi-loci technique AFLP, making impossible genetic discrimination of the analyzed trees.This lack of genetic variation together with the relevant phenotypic plasticity displayed by this species [50], supports Pinus pinea as a suitable model for studying epigenetics and its ecological and evolutionary implications [24].In this study we have analyzed DNA cytosine methylation patterns of 20 selected Pinus pinea individuals from natural populations covering the area of distribution of the species in Spain, as a first attempt to characterize this species.Due to the limited number of cDNAs that are publicly available for this nonmodel tree species, MSAP technology was used to analyze the methylation status of anonymous cytosines in P. pinea genome.We have detected that 64.36% of the analyzed cytosines at CCGG motifs were methylated.The observed percentage of stone pine genome-wide cytosine methylated sites was at least 10% higher than those reported for both annual and other perennial plants [6,12,21,[65][66][67].This result is in agreement with genome hypermethylation of conifer genomes proposed by Nystedt et al. [68] as one of the mechanisms underlying conifer genome evolution.It is interesting to mention that most of the markers that were not selected for the analysis (DNA fragments undoubtedly scored in less than 95% of the analyzed samples) were Methylation Sensitive.Additionally, the degree of cytosine methylation is expected to be even higher since fully methylated sites cannot be detected with the MSAP technique [21,66].The broad presence of DNA methylation found in Pinus pinea genome might be related to certain extent with the repetitive nature of conifer genomes [68].DNA methylation is known to control activity of mobile elements, protecting plant genomes against their mobilization [4,69] and transposable elements are abundant in conifer genomes.Pseudogenes, which are also very abundant in conifers, have been described in mammal genomes as elements encoding long noncoding RNAs involved in the epigenetic regulation of gene expression [70].
This study has also revealed a high level of PMS markers meaning that 42.13% of the total MSAPs analyzed (or 65.46% of MS markers) showed cytosine methylation variation.Taking into account that this study comprises ramets belonging to 20 trees from five populations, this result becomes especially significant to picture the cytosine methylation landscape of the species.Several studies suggest that cytosine methylation variability in particular, and epigenetic variability in general, may be associated with phenotypic plasticity in traits with potential for improving local adaptation [12,18,65].It has been shown in Arabidopsis thaliana that epigenetic diversity favors functional diversity associated with productivity and stability of populations [71].The epigenetic variation found in Pinus pinea may play a role in the fitness of the species by acting as an alternative source of variability, different to genetic diversity, with evolutionary consequences.
Better understanding of this high DNA methylation variation can be achieved by comparing isoschizomers profiles associated to EcoRI/HpaII and EcoRI/MspI digestions.Most of the cytosine methylation variation detected (90% of the PMS fragments) was associated with EcoRI/HpaII profiles, as the result of fragment detection due to HpaII digestion of hemimethylated external cytosines at the CCGG sites ( m CNG), that other samples lack because of fully methylation of the inner (C m CGG) or both cytosines ( m C m CGG) at the corresponding sites.On the other hand, the number of MSAPs present in the EcoRI/MspI profile was higher (117 vs 94) and most of the MMS fragments (91.67%) were found associated with MspI digestion.This lower level of variation may be associated with a higher percentage of fully methylated cytosines at CG motif (C m CGG) in Pinus pinea genome.These results are in agreement with the highest levels of cytosine methylation found in the CG motif observed in both plants and animal genomes [72].In Arabidopsis it is mainly found in heterochromatic regions with transposable elements and repeats, as well as in genic regions.Cytosine methylation in CNG sequence motif (where N denotes A, T or C; in our case N is C), which is also frequent in plant genomes, is associated with histone modification and involved in small non-coding RNA biogenesis in Arabidopsis [73].From an adaptive perspective, modification of their methylation status may allow trees to rapidly respond to abrupt changes in environmental conditions as well as to deal with long term responses to more general environmental scenarios [74].
The extent of cytosine methylation at CCGG sites was statistically different among stone pine individuals, ranging from 42.73% to 47.90%.PMS fragments among vegetatively propagated plants obtained from each original tree were used to estimate variability among the 20 trees initially analyzed.Polymorphism levels ranged from 0.46% (Tor 3) to 9.72% (Don 13).This variability may be associated with the developmental stage of the plant or/and differences in their growing environmental condi-tions.Although all clonally propagated trees shared their chronological stage, methylation variability may be in part associated with differences in their ontological stage, since each ramet derived from a different branch of the corresponding oneyear-old mother tree, developing a specific rooting pattern.Additionally, soil heterogeneity among plant pots and microenvironmental variation among plants due to the block design may be associated to some extent with methylation variability.
It is known that MSAP technique could help assessing genetic variability of the analyzed samples since PMI fragments can be associated with mutations on the restriction site in individuals lacking the fragment.However, considering the absence of genetic variation in the species, also supported by the AFLP analysis, PMI-MSAP fragments (5.09% of the total number of MSAPs) detected in trees from five populations of stone pine should be mainly associated with fully methylated m C m CGG restriction sites, which are demethylated in those individuals with fragment presence.All these results indicate that epigenetic variability is independent from genetic variability in this species and therefore underscore the potentially important role of the epigenetic variability as an evolutionary mechanism [24,75].
To acquire a better understanding of the evolutionary implications of this biological process, additional experiments are required to study modification of the cytosine methylation status in response to different environmental conditions (i.e.drought, different atmospheric CO 2 concentrations, etc.) as well as transgenerational inheritance of these epigenetic marks at both genome -and candidate gene-level.
In relation with genomic resources and from an agronomic point of view, the lack of genetic variation has limited the advance of breeding programs of the species.P. pinea is an economic important tree species mainly due to its edible seeds and elite clones for cone production are cultivated in grafted plantations using mainly non-clonal rootstocks [76].Limitation to identify both, elite tree cuttings and rootstocks, makes more difficult a reliably selection and effective discrimination of materials with a superior capacity for pine nuts yield.The MSAP technique opens an alternative that is worth exploring in order to achieve tree discrimination.In this study we detected a total of 15 PMS that were present or absent in all propagated trees obtained from each original tree (i.e.same profile among ramets from a mother tree but different profiles among different mother trees) that allow to distinguish 14 out of the 20 original trees analyzed.Due to the reduced set of markers in comparison with the number of analyzed samples the epigenetic relationships among the 14 trees were not determined, as indicated by the low bootstrap values obtained at most of the nodes.Additional PMS markers with potential discriminant power could be identified using additional primer combinations.A suitable number of markers can potentially be useful for elite tree identification, supporting stone pine breeding programs with a reliable method to identify improved materials.However, DNA methylation status of cytosines at target CCGG restriction sites from a given organ may vary, as mentioned above, due to plant ontogeny or environmental changes.It is therefore critical to determine the stability of any selected PMS-MSAP marker in different developmental stages and contrasting growing conditions.
Population differentiation for conservation purposes is also a major issue in this species.Different provenances have been identified along the Spanish natural distribution based on environmental characteristics (climatic and geographic) and historical human intervention (fires, clear-cuts, reforestations) but without a genetic structure supporting it [77].Recently, a set of nuclear microsatellites with medium-low or low polymorphic information content have been identified and used to analyze, in a broad sense, stone pine population structure [45].Additionally, inter-population variability has been described for growth related phenotypic traits in common garden assays [48,49].MSAP analysis offers the opportunity to study a source of variability unexplored to date [78].Although the low number of individuals per population in this work does not reach the typical approach from population genetic studies, a preliminary analysis showed epigenetic differences among populations.AMOVA and PCA performed over the two Spanish populations represented in this study with a higher number of trees, Tordesillas and Bogarra, showed that MSAP fragments were informative enough to clearly differentiate them, in contrast with the single AFLP pattern shared between all the samples that made it impossible to distinguish both populations.PCA results showed how ramets from each population clustered together along the first component in a twodimensional scatter plot.In addition, it was possible to identify smaller clusters of ramets corresponding to propagated trees.Even though genome-wide methylation levels were similar among populations, a high percentage of polymorphic MSAPs showed significant epigenetic differentiation between these populations.
Epigenetic variability has been suggested to contribute to the phenotypic plasticity and adaptive potential of individuals and populations and therefore to their evolution [16,79].Several studies have suggested that epigenetic variability alone can cause heritable variability in phenotypic traits [5,14,19,27] although its effect on fitness has still to be elucidated.Pinus pinea is a genetically depauperated but plastic species.Our results reveal a high level of cytosine methylation in stone pine genome as well as high levels of variation of methylation between the analyzed trees.These results, together with the high levels of phenotypic plasticity observed in the species, may suggest a potential role of cytosine methylation in the regulation of gene expression and variation in phenotypic traits to improve Pinus pinea fitness under different environmental conditions.Further analysis of methylation pattern evolution in stone pines subjected to different forecasted environmental conditions, whether isolated events or recurrent stresses, associated with different future scenarios (i.e.water availability, different atmospheric CO 2 concentrations, temperature, etc.), should be carried out to confirm this hypothesis.Table S1 Location, climatic characteristics and number of propagated trees and ramets per tree of the studied populations.

(PDF)
Table S2 Sequences of adaptors and primers used in MSAP and AFLP assays.

(PDF)
Table S3 MSAP markers showing statistically significant epigenetic differentiation among populations.

(PDF)
Table S4 Binary matrix codifying the scoring of the MSAP fragment patterns for both primer combinations.First and second numbers on each cell represent the independent scoring of EcoRI/HpaII and EcoRI/MspI patterns, respectively.The scoring codes were 1 for fragment presence, 0 for fragment absence, and 9 for missing data.(XLSX)

Figure 3 .
Figure 3. UPGM tree for genotype identification.Genetic similarity was calculated and bootstrapped UPGMA clustering was performed for genotype discrimination.Bootstrap computation percentages are shown over the different branches.Tree has been condensed a 25% and clones from the same genotype clustering together are labeled under the genotype code name for easier visualization.doi:10.1371/journal.pone.0103145.g003

Figure 4 .
Figure 4. Two dimensional PCA scatter plot for population differentiation.Principal Components was performed to analyze samples belonging to Tordesillas and Bogarra populations using MSAPs with significant epigenetic effect differentiating populations and propagated trees.Bogarra ramets are identified by squares and Tordesillas ramets by triangles.Propagated trees which ramets clustered together are highlighted with circles.First component (X axis) gathers 21.60% of the variation and second component (Y axis) 13.64.doi:10.1371/journal.pone.0103145.g004

Figure
Figure S1 Map showing the natural distribution of Pinus pinea L. and the location of the studied populations.Map source: EUFORGEN (modified).(PDF)

Table 1 .
Quantification of cytosine methylation in all analyzed genotypes.
Percentage of cytosine methylation and number of polymorphic fragments are provided for each genotype.doi:10.1371/journal.pone.0103145.t001