1. Introduction

Cytosine DNA methylation represents a major molecular mechanism that provides critical genomic information and contributes to understand the molecular basis of phenotypic variation based on epigenetic modification in eukaryotic organisms, though the levels and patterns of cytosine DNA methylation appear to vary drastically among different organisms (Avramidou et al., 2015; Okano et al., 1999). Cytosine-5-methylation is the most common DNA modification, which transfers a methyl group from S-adenosyl-L-methionine to fifth carbon of cytosine residues, extensively occurring in eukaryote genomes (Finnegan et al., 2000). In plants, DNA methylation generally refers to three different nucleotide sequence contexts: CG, CHG and CHH methylation (where H = C, T or A) (Gruenbaum et al., 1981; Law and Jacobsen, 2010). Many studies have showed that DNA methylation might play indispensable roles in regulating diverse physiological processes during plant growth and development, including coordination of gene expression to adapt environmental changes (Zhang et al., 2016), involvement in chromatin organization (Chinnusamy and Zhu, 2009), and maintenance of genomic integrity (Wang et al., 2016). Recent studies have showed that environmental stimuli (including biotic and abiotic stresses) can induce changes of genomic DNA methylation (Dowen et al., 2012; Garg et al., 2015), and such environmentally induced methylation can be heritable, which may result in a re-shaping of phenotypes in next generations (Herman and Sultan, 2016). Particularly, the combination of genotype-specific environmental sensitivity, DNA methylation and transgenerational stability, and phenotypic impacts might be a primary candidate mechanism for adaptive evolution in plants (Alsdurf et al., 2013; Herman and Sultan, 2016; Zhang et al., 2013). Investigation of the inherent relationships of phenotypic and DNA methylation variations might be a good way to dissect this potential mechanism for plant adaptive evolution.

Castor bean (Ricinus communis L., Euphorbiaceae), a typical nonedible oilseed crop, displays rich phenotypic variations within the intra-species. Owing to its seed oils (made of high ricinoleic acids) are often used in aviation oil, lubricants, nylon, dyes, inks, soaps, adhesive, biodiesel and so on, with its high economic value. Castor bean is widely cultivated in tropical, sub-tropical and temperate countries, particularly India, China and Brazil (Akpan et al., 2006; Alexandrov and Karlov, 2016; Sujatha et al., 2008; Webster, 1994). Because of the limited domestication, wide seed dissemination in the ancient time, and the ease and rapidity with which it becomes established in margin lands, landraces of castor bean germplasm occur in the warm regions worldwide with rich phenotypic variations in its growth habit, color of foliage and stems, seed size and oil content (Wolf, 2000). Most types are large perennials that often develop into small trees in tropical or subtropical areas; however it is usually shorter and smaller and grown annually in areas prone to frost. It is obvious that castor bean exhibits great phenotypic diversity and phenotypic plasticity to environmental factors (Qiu et al., 2010). However, previous studies have showed that the genetic diversity of castor bean appears relatively low, based on AFLP, SSR and SNP studies (Allan et al., 2008; Foster et al., 2010; Qiu et al., 2010). The genetic or molecular basis of phenotypic variation and plasticity remains unknown. Epigenetic variations may create a potential power to drive the phenotypic diversity and plasticity in castor bean.

As mentioned above, DNA methylation is the main type of epigenetic regulation. The methylation-sensitive amplification polymorphism (MSAP) technique, a modified AFLP (amplified fragment length polymorphism) method, has been well developed for determining the global DNA methylation level without any prior knowledge genomic sequences (Reyna-López et al., 1997; Vos et al., 1995; Xu et al., 2000). Owing to the two methylation-sensitive restriction enzymes HpaII and MspI, which replaces restriction enzymes used in AFLP, HpaII and MspI can recognize unmethylated tetranucleotide sequence 5'-CCGG-3' loci, but be sensitive to different methylation contexts. Specifically, HpaII can recognize and cleave the hemi-methylated sequence (^mCCGG, only one DNA strand methylated), whereas MspI can recognize and cleave fully methylated CCGG sequence (C^mCGG, methylation of both DNA strands). Thus, different methylation status at the cytosines of the CCGG sites could be recognized by the different amplification products obtained from EcoRI–HpaII and EcoRI–MspI enzyme digestion files (Lee et al., 1994; Walder et al., 1983; Xu et al., 2000). Because of its advantages including low cost, operation easy, high sensitivity and high polymorphism, MSAP technique has become a powerful and highly informative tool and been extensively applied to detect global DNA methylation level in plants in last two decades. These studies were mainly focused on detecting the DNA methylation inheritance such as in maize and rice (Lauria et al., 2004; Xiong et al., 1999; Zhao et al., 2007); investigating the epigenetic diversity at the population level such as in mangrove plants, cocoa, white polar and violet (Herrera and Bazaga, 2010; López and Wilkinson, 2010; Liramedeiros et al., 2010; Ma et al., 2013); inspecting the global DNA methylation variation with biotic or abiotic stresses such as in rice, wheat, cotton and rapeseed (Dai et al., 2015; Hong et al., 2013; Marconi et al., 2013; Zhong et al., 2009). However, there were relatively few studies on investigating the inherent relationships between phenotypic diversity and DNA methylation variation within intra-species. In this study, the global DNA methylation variations among the sixty castor bean accessions, which exhibited rich phenotypic diversity and were collected worldwide, were investigated by using MSAP method in order to explore relationships between phenotypic diversity and DNA methylation variation. This study provides new data in understanding the formation of phenotypic diversity and its epigenetic basis in castor bean.

2. Material and methods 2.1. Sample collection and DNA extraction

Sixty worldwide accessions of castor bean were collected from the USDA National Plant Germplasm System (http://www.ars-grin.gov/npgs) and our collected landraces in China, which covered the 29 countries. The detailed information for each accession used in this study is listed in Table 1. Meanwhile, we measured the hundred-grain weight and oil content of castor bean seeds for each accession by minispec mq-one seed-analyzer (Bruker Optik GmbH, Germany).

All seeds were germinated at greenhouse and the leaves of three-leaf stage seedlings were collected. High quality genomic DNA of each individual was then extracted using the Plant Genomic DNA Kit (TianGen, Beijing, China). The DNA quality was confirmed by 0.8% agarose gel electrophoresis, and the DNA concentration was estimated by NanoDrop2000 (Thermo, USA).

2.2. Primer screening and MSAP analysis

In this study, MSAP technique was used to investigate epigenetic polymorphism of different castor bean germplasms (Reyna-Lđpez et al., 1997). Briefly, high-quality genomic DNA was digested with two sets of restriction enzyme combinations, EcoRI/MspI and EcoRI/HpaII. The digestion reaction mix contained 1 ml genomic DNA, 10 × buffer, 8 U EcoRI, 16 U MspI or HpaII (New England Biolabs, Beverly, MA, USA). The digestion reaction was carried out at 37 ℃ for 9 h and ended at 70 ℃ for 10 min. Subsequently, the 4 ml digested DNA fragments were ligated to the specific adapters in a reaction mix containing 5 pmol EcoRI adapter, 50 pmol MspI/HpaII adapter, 2 ml ligase buffer and 1 U T4 DNA ligase (New England Biolabs, Beverly, MA, USA) that was incubated at 16 ℃ for 10 h.

The 2 ml ligated DNA was used as template for pre-amplification in a reaction mix containing 2.5 ml 10 × PCR buffer, 2 ml dNTP (2.5 mmol/ml), 1.5 ml Mgcl2 (25 mmol/ml), 1 ml pre-selective primer (see Table 2), and 1 U of Taq Polymerase (Takara Biotech. Inc., Dalian). The PCR conditions were following profile: 94 ℃ for 3 min, 26 cycles at 94 ℃ for 30 s, 30 s annealing at 56 ℃ and 60 s extension at 72 ℃, final extension at 72 ℃ for 5 min. The diluted (40-fold) amplified products were used as the template for selective amplification. A total volume of 25 ml PCRs were carried out with 20 selective primer combinations (see Table 2). The touch-down PCR conditions were performed and reaction conditions were following profile: 94 ℃ for 3 min followed by 12 cycles of 94 ℃ for 30 s, 65 ℃ for 30 s (reduced 0.7 ℃ per cycle), and 72 ℃ for 1 min, and a subsequent 23 cycles of 94 ℃ for 30 s, 56 ℃ for 30 s, and 72 ℃ for 1 min. A fluorescence-based MSAP technique was used for fragment separation and detection by using an ABI3730 × l DNA analyzer (Applied Biosystems, Milano, IT). The fragments were scored and translated to binary data (0 for absence and 1 for presence) in the EcoRI/MspI and EcoRI/HpaII samples by GeneMarker 2.2.0 software. Meanwhile, for polymorphic primers, selective amplification products were separated on 6% polyacrylamide gels at 180 V for 1.5 h, and then the gels were stained with silver nitrate.

2.3. Scoring of MSAP bands and data analysis

According to different sensitivities to methylation of restriction enzyme HpaII and MspI to the CCGG sequence, the methylation states could be classified into following four basic types based on previous reports (Portis et al., 2004): Type Ⅰ pattern represents the non-methylation sites (E-H/E-M, 1/1), Type Ⅱ pattern represents the hemi-methylation sites (E-H/E-M, 1/0), Type Ⅲ pattern represents full-methylation sites (E-H/E-M, 0/1), and Type Ⅳ pattern represents the uninformative sites (E-H/E-M, 0/0).

For each primer, the well-resolved bands were scored as absence (0) and presence (1) with same molecular weight to form a binary matrix. Epigenetic diversity parameters, including the percentage of polymorphic loci (PPB), the expected heterozygosity (He), Shannon's Information Index (I), the effective number of alleles per locus (Ne) were estimated and analyzed by using the GenAlEx version 6.501 software (Rod and Petere, 2006). Based on the binary matrix, the Nei's genetic distance matrix was calculated by DPS (Data Processing System) software (Vision 7.05), and then the genetic distance matrix was transferred to MEGA6 software (Tamura et al., 2013) to construct the phylogenetic tree with neighbor-joining method.

To determine population differentiation, we used the STRUCTURE program which employs a Bayesian algorithm to calculate the number of populations. The number of the MCMC reps was set to 100, 000 and the length of the burn-in period was set to 50, 000 steps. In order to estimate the appropriate number of clusters (K), ΔK was calculated by the method of Evanno et al. (2005) with ten runs for each K ranging from K = 1 to K = 9.

2.4. Sequencing and annotation of polymorphic MSAP fragments

Some bands showing differential DNA methylation level among castor bean accessions were selected and excised from polyacrylamide gel. The bands were rehydrated with 20 ml of sterile distilled water overnight at 4 ℃. The soluble DNA samples (about 2 ml) were used as PCR template for re-amplified by using the corresponding selective amplification primers and PCR conditions. Re-amplified products were checked by agarose gel electrophoresis, and were purified by using the Gel Extraction Kit (GENEray, China). The purified DNA fragments were then cloned into pMD19-T cloning vector (Takara Biotech. Inc., Dalian), and then transferred into Escherichia coli. Trans1-T1 (TransGen Biotech, China). Positive clones that were confirmed by colony PCR, were selected and sequenced. The obtained sequence information was BLASTN searched against the castor bean genome database including nuclear, mitochondrial and chloroplast genome in NCBI (http://www.ncbi.nlm.nih.gov/).

3. Results 3.1. The DNA methylation level and epigenetic diversity of castor bean germplasms

To obtain the efficient and polymorphic primers, we initially screened PCR products for 8 random selected DNA samples using polyacrylamide gel electrophoresis. Initially, seventy-two EcoRI/ HpaII and EcoRI/MspI primer combinations were tested. The primer combinations which can produce arbitrary, clear and reproductive PCR bands were sort out for further analysis. Finally, twenty primer combinations were determined to test global DNA methylation variation of all accessions. Since the EcoRI primers were fluorescent-labeled, the PCR products were accurately analyzed with ABI3730 × l DNA analyzer. To exclude the noise peaks, the fragment size of PCR products out of the ranging from 50 bp to 500 bp was excluded for the further analysis.

In total, the 20 primer sets produced the reproducible PCR product fragments from the 60 accessions tested, ranging from 123 to 292 fragments among different landraces. The total number of bands ranged from small to large in order according to the size of bands, and these polymorphic fragments obtained from PCR products were combined into a binary matrix to indicate the presence (1) or absence (0) of the fragments in EcoRI/HpaII and EcoRI/MspI digestion profiles. According to the MSAP profiles, the total number of non-methylated sites (CCGG, double-strand; 1, 1), hemi-methylated sites (^mCCGG, only one strand methylated; 1, 0) and fully methylated CCGG sites (C^mCGG, double-strand methylated; 0, 1) were counted respectively. The relative total levels of non-methylation or methylation were calculated as the percentages non-methylation or methylation sites (includes the total number of hemi-methylated sites and fully methylated sites) at the total number of all bands. The total level of DNA methylation ranging from 3.80% to 34.31% of sixty castor bean accessions, of which the fully methylation level ranging from 3.80% to 32.31%. The methylation level of castor bean came from Cuba (Cu-1) was 3.80% and castor bean of Uruguay (Ur) was the highest (34.31%) in all castor bean germplasms. This indicated that DNA methylation level varies greatly among different landraces. In addition, our results showed that the fully methylated ratios (Type Ⅲ, methylation of the internal cytosines on both strands) were higher than the corresponding hemi-methylation levels (Type Ⅱ, hemi-methylation on the external cytosine) of all castor bean germplasms that we adopted.

In order to assess the epigenetic diversity of different castor bean germplasms, the binary matrix (zero and one) was subject to the estimation of the basic parameters including Na (No. of Different Alleles), Ne (No. of Effective Alleles), I (Shannon's Information Index), He (Expected Heterozygosity), uHe (Unbiased Expected Heterozygosity) and the PPB (the percentage of polymorphic loci) using the GenAlEx software. As a result, we found that the value of PPB was 69.07%, and the values of Na, Ne and uHe were 1.691, 1.395 and 0.244, respectively. The average I per locus was 0.366 ranging from 0.389 to 0.693, and the average of He value per locus was 0.242 ranging from 0.218 to 0.498.

3.2. Population structure analysis

STRUCTURE is a Bayesian and model-based algorithm which was used for clustering genetic data analysis. Given the number of clusters (K) and assuming Hardy–Weinberg and linkage equilibrium within clusters, STRUCTURE estimates allele frequencies in each cluster and population memberships for every individual (Hubisz et al., 2009). Here, the Structure software was used for estimating the appropriate number of clusters (K) and inferring the population structure. The admixture model was described by Evanno et al. (2005), the most optimal number of cluster for the whole dataset was determined by the highest ΔK value, whose corresponding the value of K was regarded as the optimal number of cluster. The result of structure analysis for the sixty samples yielded a highest ΔK value (ΔK = 333.625) which showed a clear peak at the K value of 2 (Fig. 1B). Thus, sixty castor bean germplasms were properly divided into two clusters, 33 castor bean accessions were included in the first cluster (Fig. 1C, the red color), collected from African, America and Asia across the main collection regions. The second cluster consisted of the 27 castor bean accessions (Fig. 1C, the green color), collected from the countries of America, Asia, Africa and Europe, similarly across the main collection regions.

Further, using the neighbor-joining criteria we constructed a phylogenetic tree based on the Nei's genetic distance matrix data. As showed in Fig. 2, the sixty accessions were clustered as two distinct clades (Clade Ⅰ and Clade Ⅱ). The Clade Ⅰ was comprised of 16 germplasms collected from Africa, South America, North America, and Asia; The Clade Ⅱ included 44 accessions collected across the collection regions. In comparison with the results obtained from the population structure analysis (Fig. 1), the two main groups were consistent though the number of members within each group were slight different. While inspecting the geographical distribution of each accession within the two groups, it was clear that there was no distinct geographic structure among castor bean.

3.3. Sequence analysis of DNA methylated loci

To investigate which loci or genes on genome resulted in the methylation polymorphism, we randomly inspect more than 20 polymorphic fragments by separating the PCR polymorphic bands from the polyacrylamide gels and clone-sequencing to obtain the DNA sequences of polymorphic bands. In total, we successfully obtained 13 DNA sequences, ranging from 48 bp to 582 bp in length (in Table 3). After the BLASTN search in NCBI database, we found that the 13 DNA fragments were successfully annotated as four distinct regions of castor bean genomes. We denoted these DNA fragments as P1–P4 region as shown in Table 3. The P1 region was identified as a fragment located in downstream 600 bp of nad2 gene in the castor bean mitochondrion genome (HQ874649.1); the P2 region including ten DNA fragments was annotated as the ndhB gene (NADH dehydrogenase subunit 2) in chloroplast genome (JF937588.1), functionally involving photomorphogenesis (Catalá et al., 1997; Guéra et al., 2000); the P3 and P4 were identified as uncharacterized mRNAs with unknown function in nuclear genome. These results indicated that the polymorphic methylated loci occur widely in nuclear, chloroplast and mitochondrion genomes.

4. Discussion

DNA methylation is an ancient epigenetic modification that appears in the genomes of organisms throughout the eukaryotic phylogeny (Feng et al., 2010; Suzuki and Bird, 2008). Increasing evidence have shown that genomic DNA methylation contribute to heritable phenotypic variation (Das and Messing, 1994; Liu et al., 2004). As mentioned above, castor bean germplasms displayed rich phenotypic diversity, but previous studies had revealed that the genetic diversity among germplasms is relatively low (Allan et al., 2008; Foster et al., 2010; Qiu et al., 2010). Why castor bean exhibits the rich phenotypic diversity with the low genetic diversity remained uncertain. This study is the first report on natural epigenetic variation occurring in castor bean germplasms. Our results revealed that the genomic DNA methylation of castor bean exhibited a considerable variation among germplasms. About ten times divergence of DNA methylation level was identified between the lowest and the highest germplasms. Although these germplasms were collected from diverse nature environment, since the germplasms were collected worldwide and germinated and grown under the same environment (the greenhouse), the exhibited difference of genomic DNA methylation is supposed to be genetically evolved, not resulting from the environmental heterogeneity. Besides, leave samples used in this study were collected at the same developmental stage so that any developmental variation in methylation, to the extent that it occurs within leaves, would not confound our results to determine genotype-specific variation in methylation patterns. Although based on current study whether the divergent DNA methylation levels occurring among castor bean germplasms is a driving power to shape the rich phenotypic diversity in castor bean remains unclear yet, the considerable variation of genomic DNA methylation might be a potential clue to dissect the molecular basis of phenotypic variation among germplasms. Using the same MSAP methods, previous studies have shown that 30–41% of loci were methylated in Brassica oleraceae accessions (Salmon et al., 2008), 32% of loci were methylated in Gossypium hirsutum accessions (Keyte et al., 2006), and 35–43% of loci were methylated in Arabidopsis thaliana ecotypes (Cervera et al., 2003). The degree of cytosine methylation may be expected to be even higher since some methylated sites cannot be detected with the MSAP technique (Sáez-Laguna et al., 2014), because detection is restricted to the recognition sites of the isoschizomers used (Xiong et al., 1999).

In this study, based on the epigenetic differentiation the two main clusters were identified. Compared to the previous study that investigated the genetic diversity among castor bean germplasms using EST–SSR markers and resulted in five distinct clusters (Qiu et al., 2010), it seems clear that the epigenetic differentiation is not dependent on genetic differentiation among castor bean germplasms. Probably, this divergence of epigenetic and genetic variations might be due to the different molecular mechanisms underlying the epigenetic and genetic variations occur. Interestingly, Rivarola et al. (2011) identified two main genetic clusters based on partial chloroplast genomic SNPs in castor bean germplasms. However, because of the limited samples tested in the study of Rivarola et al. and the different sample strategy from this study we could not figure out the relationships of epigenetic differentiation and genetic variations in castor bean. Combining our current study and with other reports (Kanchanaketu et al., 2012; Shan et al., 2012), it seems to be unexpected to learn that genetic and epigenetic diversity is not always linked together. In addition, the previous studies based on genetic variations revealed that there were no distinct geographic groups in castor bean germplasms (Allan et al., 2008; Qiu et al., 2010). Also, the current results indicated that epigenetic differentiation did not display a distinct geographic structure. This result is understandable because most of castor bean germplasms were multi-introduced across several continents with human activities, although castor bean landraces did not cultivated in many regions.

In addition, thirteen methylated DNA fragments were isolated from polymorphic methylated loci among castor bean accessions, and were mapped into four distinct regions on castor bean genomes (including nucleus, mitochondrion and chloroplast) in this study. This might be related to tissue-specific DNA methylation because our tested samples are young leaf tissues. Interestingly, ten DNA fragments with different length were resided in the same region of chloroplast genome and were identified as the gene ndhB, suggesting that this region or the gene ndhB was sensitive to DNA methylation among natural population. Although there were studies that found the DNA hypomethylation of chloroplast genome in Arabidopsis (Leutwiler et al., 1984), the DNA methylation degrees of chloroplast genomes might be variable among different plants because of some former reports on DNA methylation variations of chloroplast genome in other plants (Kobayashi et al., 1990; Ohta et al., 1991; Takio et al., 1994). The other two methylated loci were identified from genome as uncharacterized mRNA and non-coding RNA. Although limited polymorphic methylated loci were here identified, it seems to imply that polymorphic methylated loci occur widely, involving nuclear, mitochondrion and chloroplast genomes.

Acknowledgments

We thank Dr. Yupeng Geng from Yunnan University for his support and Ms. Ruiwen Li from Yunnan University for her assistance with the polyacrylamide gel electrophoresis experiment. This work was jointly supported by Chinese National Key Technology R & D Program (2015BAD15B02) and National Natural Science Foundation of China (31661143002 and 31501034).

Appendix A. Supplementary data

Supplementary data related to this article can be found at http://dx.doi.org/10.1016/j.pld.2017.05.007.