ORIGINAL RESEARCH article

Front. Plant Sci., 26 April 2016

Sec. Plant Nutrition

Volume 7 - 2016 | https://doi.org/10.3389/fpls.2016.00447

Analyses of Methylomes Derived from Meso-American Common Bean (Phaseolus vulgaris L.) Using MeDIP-Seq and Whole Genome Sodium Bisulfite-Sequencing

  • 1. Molecular Genetics and Epigenomics Laboratory, Delaware State University, Dover DE, USA

  • 2. Division of Science and Mathematics, Mayville State University, Mayville ND, USA

  • 3. Center for Integrated Biological and Environmental Research, Delaware State University, Dover DE, USA

Abstract

Common bean (Phaseolus vulgaris L.) is economically important for its high protein, fiber, and micronutrient contents, with a relatively small genome size of ∼587 Mb. Common bean is genetically diverse with two major gene pools, Meso-American and Andean. The phenotypic variability within common bean is partly attributed to the genetic diversity and epigenetic changes that are largely influenced by environmental factors. It is well established that an important epigenetic regulator of gene expression is DNA methylation. Here, we present results generated from two high-throughput sequencing technologies, methylated DNA immunoprecipitation-sequencing (MeDIP-seq) and whole genome bisulfite-sequencing (BS-Seq). Our analyses revealed that this Meso-American common bean displays similar methylation patterns as other previously published plant methylomes, with CG ∼50%, CHG ∼30%, and CHH ∼2.7% methylation, however, these differ from the common bean reference methylome of Andean origin. We identified higher CG methylation levels in both promoter and genic regions than CHG and CHH contexts. Moreover, we found relatively higher CG methylation levels in genes than in promoters. Conversely, the CHG and CHH methylation levels were highest in promoters than in genes. This is the first genome-wide DNA methylation profiling study in a Meso-American common bean cultivar (“Sierra”) using NGS approaches. Our long-term goal is to generate genome-wide epigenomic maps in common bean focusing on chromatin accessibility, histone modifications, and DNA methylation.

Introduction

Common bean (Phaseolus vulgaris) is the most widely consumed legume, as it is high in protein, fiber, and essential nutrients, while low in glycemic index (Mitchell et al., 2009). The haploid genome size of common bean is about 587 Mb (Schmutz et al., 2014). Bean domestication is assumed to be a complex process that involved two distinct gene pools (Meso-American and Andean) and repeated selection of desirable traits within the gene pool (Bitocchi and Nanni, 2012). Common bean is consumed widely across the world due to its affordability and health benefits (Gepts et al., 2008). Food legumes such as common bean are an inherently rich source of vitamins and minerals, including B vitamins (particularly B1/thiamine and B9/folic acid), iron, zinc, calcium and magnesium (Akond et al., 2011; Beebe et al., 2013; Câmara et al., 2013; Petry et al., 2015). Hence, common bean is considered as a viable target for biofortification (Blair, 2013). In order to effectively modulate genetic architecture and develop micronutrient-enriched crops, understanding the genes and their regulatory mechanisms involved in physiological processes is important. Though several regulatory mechanisms have been identified, DNA methylation, chromatin remodeling, histone acetylation, histone methylation, and gene silencing are of paramount importance.

DNA methylation is a covalent, heritable epigenetic modification that plays a significant role in gene expression, tissue specialization, and transposon inactivation (Zhang et al., 2006; Verhoeven et al., 2010). In plants, cytosine methylation occurs in both symmetrical (CG and CHG) and asymmetrical (CHH) contexts, while descending order of the extent of methylation in these contexts is CG, CHG, and CHH, where H indicates A, C, or T (Lister and Ecker, 2009; Henderson et al., 2010; Saze et al., 2012). DNA methylation occurs at higher levels in heterochromatin when compared to euchromatin (Madlung and Comai, 2004; Song et al., 2013). Genome-wide DNA methylation and gene expression studies revealed that gene-body methylation is conserved across species (Feng et al., 2010), among constitutively expressed genes (Zhang et al., 2006), and more specifically, CG methylation is often linked to increased gene expression (Miura et al., 2009; Wang et al., 2014; Kim et al., 2015). Conversely, promoter DNA methylation in plants is largely associated with transcriptional repression (Li et al., 2008) and is highly tissue-specific in nature (Zhang et al., 2006; Song et al., 2013). In soybean (dicot) and rice (monocot), it has been reported that promoter hypomethylation resulted in increased gene expression of flanking genes (Li et al., 2008; Song et al., 2013).

There are several methods for determining the presence and location of methylated cytosines in the genome, which include methylation-sensitive restriction enzyme sequencing (MRE-seq), affinity purification, and whole genome bisulfite-sequencing (BS-Seq). This study used one of the three affinity purification methods, i.e., methylated DNA immunoprecipitation sequencing (MeDIP-seq) and BS-Seq. MeDIP-Seq utilizes anti-5-mC antibody and generates ∼107 reads while BS-Seq being a “gold standard” of methylation sequencing, generates ∼108 methylated cytosines at single nucleotide resolution. This technique converts unmethylated cytosines to uracil by deamination, while methylated cytosines remain as cytosines (Lister and Ecker, 2009; Laird, 2010; Kim et al., 2014). A genome-wide DNA methylation map has been reported in Andean common bean (Kim et al., 2015) but it is not available in Meso-American common bean. The previous methylome study compared the methylation in Andean common bean leaf to soybean leaf, stripped root, and root hair methylomes (Kim et al., 2015). Screening the methylomes of Meso-American common bean (P. vulgaris) using high-throughput approaches such as BS-Seq and MeDIP-Seq for genome-wide DNA methylation may aid in understanding key regulatory mechanisms linked to important biological processes.

Results

Whole genome bisulfite-sequencing (BS-Seq) is widely used to study DNA methylation at the whole genome level while MeDIP-Seq is used to detect high DNA methylation in low-density CG/CHG/CHH areas or moderate methylation levels in high-density CG/CHG/CHH regions. To understand the DNA methylation levels in the three contexts (CG, CHG, and CHH), we constructed two MeDIP-Seq libraries (one input control that served as no antibody sample and one immunoprecipitated sample) and one BS-Seq library using tissue collected from three sets of bean plants from the Meso-American common bean cultivar, Sierra, grown at three different times (three biological replicates grown in triplicate) which were pooled and sequenced on the Illumina/HiSeq-2500 platform. A similar replicate pooling approach has been used in other epigenomic studies that utilized next generation sequencing technologies (Taylor et al., 2007; Brinkman et al., 2012; Statham et al., 2012). The libraries generated in this study were labeled as Sierra_Input, Sierra_MeDIP and Sierra_BS (Table 1). The input control (Sierra_Input) was used for comparing against immunoprecipitated DNA (Sierra_MeDIP) during MeDIP-Seq analysis to identify differentially methylated regions (DMRs).

Table 1

Sample nameMethodTotal read #Mapped read #Mapping ratio (%)
Sierra_BSBS-Seq214,016,81961,362,99828.67
Sierra_InputINPUT22,509,8738,897,76439.53
Sierra_MeDIPMeDIP-seq23,709,4237,173,10430.25

Reads and mapping ratios for each sample.

Data Collection and Preprocessing

Deep sequencing of three methylome libraries in common bean resulted in ∼260 million-50 bp Illumina reads (Table 1). Of these, ∼214 million reads were from BS-Seq and the rest (∼47 million reads) were from MeDIP-Seq. Among MeDIP-Seq reads, ∼22.5 million reads were from input control and ∼23.7 million were methylated DNA immunoprecipitated reads. The raw reads obtained from both methodologies were trimmed, filtered, and high quality reads collected were aligned to the reference P. vulgaris genome, G19833 (Schmutz et al., 2014) available at Phytozome (V1.0, accessed October 2014). In total, ∼61.4 (28.67%), ∼8.9 (39.53%) and ∼7.2 (30.25%) million BS-Seq, MeDIP-Seq input control, and MeDIP-Seq immunoprecipitated reads were mapped to the reference genome (G19833). Of these, uniquely mapped reads with ≤2 mis-matches were further used in analysis.

Bisulfite-Sequencing (BS-Seq)

DNA methylation ratios from BS-Seq were evaluated for CG, CHG, and CHH contexts. The increasing order of methylation levels in the three contexts were CHH followed by CHG and CG. Though the methylation patterns for CG and CHG contexts were similar, the overall levels of CHG methylation were relatively lower than CG methylation. The genome-wide DNA methylation levels for each context are shown (Figures 1A,B). The least methylated context, CHH was consistently very low across all 11 common bean chromosomes. In addition to the overall methylation across the three contexts, we also investigated the methylated cytosines in promoter and genic regions. The promoter and gene coverage between CG (67%, 71%) and CHG (71%, 73%) contexts varied when compared to CHH coverage (95%, 90%) context. While the majority of the reads covered in all contexts were more than 50-fold at each cytosine site, some sites were not covered (Figure 2). For each context, the overall read coverage (number of reads per cytosine) was determined as either 0 (no coverage), 1–4, 5–10, 11–20, and >20 reads per cytosine (Supplementary Figure S2). The majority of the covered reads was more than 11 reads per cytosine and this trend was consistent across the three contexts of methylation measured by BS-Seq.

FIGURE 1

FIGURE 2

Methylation ratios and their frequencies throughout the genome were analyzed in each context. Within the CG context, ∼2.65 million total sites and ∼1.25 million methylated sites ≥ 0.8 (∼50%) were identified across the genome. In the CHG context, ∼3.3 million total sites and ∼one million methylated sites (∼30%) were identified. The most prominent context was CHH and we identified approximately 18 million CHH sites in the genome, but only ∼500,000 sites were methylated (2.7%) with the vast majority of sites unmethylated. Our results show that most of the methylation and highly methylated regions are found outside of the promoter and genic regions (intragenic) in all three contexts (Figure 3A). Methylation levels for the annotated promoter and genic regions are presented in box plots for each context (Figure 3B). Higher CG methylation was found in genic regions when compared to the promoters. Conversely, comparatively higher CHG and CHH methylation was observed in promoter regions than genic regions. However, the overall CG methylation was highest among the three contexts.

FIGURE 3

Methylated DNA Immunoprecipitation (MeDIP-Seq)

To understand the methylome of common bean, assessing the peak density and peak shape associated with DNA methylation is necessary. Peaks can be called by mapping the reads to the genome to reveal the loci of selectively methylated DNA. Several peak callers are currently available to identify the global pattern of these modified sites by comparing the antibody tagged sample with its background (input, no antibody). The shapes are based on genomic coverage of the input control, the coverage of the DNA that was not subjugated to immunoprecipitation with an anti-5mC antibody, versus the immunoprecipitated DNA. The increased coverage is displayed as a “peak,” which implies that there is more coverage, due to the enrichment of DNA containing methylated cytosines (Supplementary Figure S1). Peak assessment is very important in minimizing the false discovery rate and developing a gold-standard common bean methylome. A correlation between peak density and chromosomal length has been proposed in mammals (Hughes et al., 2010). We estimated the methylation peak density (peak/Mb) based on chromosome length and total number of peaks identified per chromosome. MeDIP-Seq signatures per chromosome, total peaks across the genome, peaks found only in promoter regions, and peaks identified in genic locations are presented (Table 2). The majority (93.2%) of the peaks identified from MeDIP-Seq analysis were in non-promoter and non-genic regions. The peaks obtained from MeDIP-Seq were compared against the methylation sites from BS-Seq data to corroborate the peaks identified in likely methylated regions and to compare the two methylome sequencing work flows. The top 20 most significant (p < 0.05) MeDIP-Seq peaks were compared to CG and CHG weighted methylation determined from BS-Seq to identify that all of these methylated cytosine regions possess methylation ratios over 0.5 in CG and/or CHG contexts (Table 3).

Table 2

ReferenceChrom length (Mbp)ReadsSignaturesTotal peaksPeaks in promoterPeaks in genes
Chr0152.21,889,08135,4363,26113080
Chr02491,413,71827,4212,97710399
Chr0352.31,711,96631,6143,03712877
Chr0445.91,916,91231,4582,9439265
Chr0540.81,548,87126,4042,57210982
Chr0632967,06517,4341,5448172
Chr0751.82,066,30034,4733,59813083
Chr0859.72,877,21234,3593,774137109
Chr0937.51,194,00816,6101,9199771
Chr1043.31,986,80929,9053,091105101
Chr1150.42,473,34534,5683,740153104
Total514.920,045,287319,68232,4561,265943

Common bean chromosomal length, reads, signatures, and total and annotated peaks from MeDIP-Seq data.

Table 3

LengthPeak shape scoreP-value5′ gene5′ distance3′ gene3′ distanceCGCHG
24315.032.30E–51Phvul.008G00010000.897252N/A
23714.092.24E–45Phvul.008G000100957070.9326420.892308
23113.668.72E–43Phvul.008G0001001120680.6862420.679463
20413.383.70E–41Phvul.004G0469002257Phvul.004G04700049360.949129N/A
21013.291.35E–40Phvul.008G29320084760.850638N/A
24612.873.41E–38Phvul.008G0001001009850.8572330.726482
20612.761.29E–37Phvul.008G29320043720.7630930.166279
22612.712.57E–37Phvul.009G14470060451Phvul.009G144800282620.8972140.823301
24212.182.05E–34Phvul.008G29320064710.94471
20912.031.23E–33Phvul.008G29320087240.92825N/A
23511.953.34E–33Phvul.008G0001001087870.8421050.88024
21211.915.27E–33Phvul.008G0001001133320.8074150.734861
22111.431.41E–30Phvul.009G14470050026Phvul.009G14480038692N/A0.775591
21011.363.29E–30Phvul.008G0001001105520.9253730.892
21411.315.68E–30Phvul.008G0001001052110.8328640.814815
22011.278.84E–30Phvul.009G14470051399Phvul.009G144800373200.50.779538
20411.011.69E–28Phvul.008G000100109425N/A0.892105
20610.831.26E–27Phvul.009G14470058862Phvul.009G144800298710.8399010.689939
27410.762.59E–27Phvul.009G14470053308Phvul.009G144800353570.8503940.730228
22310.352.11E–25Phvul.008G29320024790.961965N/A

The top 20 most significant/highest peak shape scores from the MeDIP-Seq data are shown.

Also shown are peak locations in relation to 5′ and 3′ genes. CG and CHG methylation ratios were determined from BS-Seq data in the corresponding peak locations.

Discussion

DNA methylation is a heritable epigenetic mark that controls a wide variety of physiological processes (Boyko et al., 2007). Though DNA methylation has been reported in many eukaryotic organisms, the contexts and levels of DNA methylation vary significantly between plants and animals. In animals, the symmetric context, CG, is most prominent while non-CG (CHG and CHH) methylation is less frequent (Stroud et al., 2013b). In plants, DNA is primarily methylated in the CG context while non-CG contexts (CHG and CHH) are also abundant. The reports on CG and non-CG methylation in plants are increasingly evident. For example, in Arabidopsis both CG and non-CG methylation are functional (Law and Jacobsen, 2010) and the roles of chromomethyltransferase (CMT) and domains rearranged methyltransferase (DRM) proteins in non-CG methylation has been established (Cao et al., 2003; Stroud et al., 2013b).

MeDIP-Seq yields a resolution of ∼100–200 bp while single base resolution methylome maps can be obtained by both conventional Sanger sequencing of short BS-treated fragments and whole genome bisulfite-sequencing (Taiwo et al., 2012; Stevens et al., 2013). In common bean, we generated 10-fold more reads using BS-Seq when compared to MeDIP-Seq, as previously reported (Laird, 2010). The relatively lower number of reads in MeDIP-Seq compared to BS-Seq is due to the inherent difference between the techniques (Laird, 2010). However, MeDIP-Seq analysis identified densely methylated regions in the genome while BS-Seq analysis estimated overall methylation rates in the genome.

Currently, few reports of BS-Seq are available in plants, which are mostly model species and some are presented (Table 4). Epigenome reports in legume crops are scantly available. Recently, single nucleotide resolution reference methylomes in soybean and Andean common bean have been constructed using BS-Seq in order to understand epigenetic variation in legume crops (Kim et al., 2015). Here we generated single nucleotide resolution (BS-Seq) and 150 bp resolution (MeDIP-Seq) maps of DNA methylation to understand genome-wide DNA methylation in Meso-American common bean. Further, we focused our analyses in estimating the overall methylation levels, context-specific methylation patterns, and promoter versus gene body methylation patterns. This study represents the first genome-wide DNA methylation report in the important rust-resistant Meso-American common bean cultivar, Sierra.

Genome-Wide Methylation Levels

We aligned the sequenced cytosines collected from MeDIP-Seq and BS-Seq to the common bean reference genome (G19833) to identify genome-wide methylation, site-specific methylation and annotated the methylated peaks. In total, 3.0–3.5 million cytosines were found to be methylated across the genome covering all three contexts of DNA methylation, thus accounting for <1% of the total DNA bases in the genome. A total number of 32,456 significant (>1.5-fold-enrichment) peaks were identified in common bean. The average DNA methylation peak density identified was ∼63 peaks/1 Mb of DNA. Based on the average peak density, lower methylation levels were identified in chromosomes, Chr02, Chr03, Chr06, and Chr09. Conversely, higher methylation levels were noted in chromosomes, Chr07, Chr10, and Chr11 (Table 2). Similar patterns in methylation saturation between chromosomes were highlighted by using chromosome-based methylation distribution heat maps (Figure 1B). Genome-wide methylation identified in this study was over 37.8% in CG and CHG contexts, which is in accordance with the other reports in Arabidopsis, rice and soybean (Cokus et al., 2008; Li et al., 2012; Song et al., 2013).

Context-Specific Methylation Patterns

In common bean, we generated methylation histograms in order to show the methylation ratios (0.0–1.0) for the cytosine sites in each context within the genome and for the annotated genomic features (Figure 2). The majority of the overall methylation ratios for hemi-methylated or fully methylated sites ranged between 0.8 and 1.0 in the CG context, which is in accordance with CG methylation in Arabidopsis. The distribution of methylation ratios in the CHG context for the hemi-methylated or fully methylated sites ranged between mainly 0.6-1.0 while in Arabidopsis the ratios reported were between 0.2 and1.0 (Cokus et al., 2008). Within bean, the CG and CHG distributions were similar except for fewer fully methylated sites (ratio = 1.0) in the CHG context. The CHH methylation ratios were predominately <10% in Arabidopsis (Cokus et al., 2008; Feng et al., 2010), maize (Eichten et al., 2013), and many other eukaryotic organisms (Feng et al., 2010) which is consistent with the present findings in common bean. As CHH methylation is commonly found in repetitive and transposable elements, it varies significantly more between tissues and even among biological replicates than CG and CHG methylation (Song et al., 2013; Kim et al., 2015).

In addition to the methylation ratios, the overall methylation levels in our study estimated for CG, CHG and CHH contexts were ∼50, ∼30, and ∼2.7%, respectively. In Arabidopsis, the methylation levels specific to each context (∼24% in CG, 6.7% in CHG, and 1.7% in CHH contexts) were reported (Cokus et al., 2008). A similar study which screened eight eukaryotic species identified ∼22% of CGs, ∼6.0% of CHGs, and ∼2% of CHHs as methylated in Arabidopsis (Feng et al., 2010). In poplar, the methylation level, 41.9%, 20.9%, and 3.25% were reported for CG, CHG, and CHH contexts, respectively (Feng et al., 2010). Recently, among legumes, the methylation percentage for CG (74% and 64%), CHG (62% and 48%) and CHH (21% and 4%) contexts have been reported in common bean and soybean leaves, respectively (Kim et al., 2015). Other leaf methylation levels reported in soybean include, 63% in CG, 44% in CHG, and 5.9% in CHH contexts (Song et al., 2013). Our results were more similar to the methylation levels for each context in poplar than in soybean. This is likely attributed to the fact that bean and poplar have a similar genome size when compared to soybean. Moreover, soybean has experienced a whole-genome duplication after diverging from common bean. Interestingly, our results slightly deviated from the earlier reports in common bean for all the three contexts (Kim et al., 2015). The deviation in the methylation ratios in this study may be attributed to the genotype-specific DNA methylation variation (Meso-American vs Andean), differences in leaf stages (primary vs. trifoliate) and collection time-points (14-day-old vs. 18-day-old) used while analyzing the bean methylome.

Promoter and Gene-Body Methylation

Promoter methylation and genic methylation have been implicated in diverse regulatory roles. Generally, cytosine methylation in promoter regions results in transcriptional repression in plants. A single methylated cytosine in a promoter region can significantly affect the expression level of the corresponding gene. About 16.6% of methylated single CGs identified in the promoter adjacent to transcriptional start sites decreased gene expression in humans (Medvedeva et al., 2014). CG and CHG promoter methylation typically causes repression of adjacent genes, as reported in Arabidopsis, maize, and soybean (Zhang et al., 2006; Song et al., 2013; West et al., 2014). The role of gene-body methylation in controlling alternative promoters and gene splicing has been implicated in mammals (Lou et al., 2014). Gene-body methylation is considered as an evolutionary consequence and it is conserved among plant orthologs (Takuno and Gaut, 2012). With this background we attempted to understand the promoter and genic methylation in common bean.

After annotating the methylated sites for genomic features, we identified higher CG methylation levels in both promoter and genic regions than CHG and CHH contexts (Figure 3A). Moreover, we found relatively higher CG methylation levels in genes than in promoters. Conversely, the CHG and CHH methylation levels were highest in promoters than in genes, which is consistent with the findings in Arabidopsis (Cokus et al., 2008). We further confirmed that the CG methylation levels are higher in the annotated genomic features (Figure 3B). The most significant MeDIP-Seq peaks were found to correlate with high levels of both CG and CHG methylation (Table 3). A study in Arabidopsis revealed that important functional genes evolved slowly when they contained methylation in the gene-body (Takuno and Gaut, 2012). Another study in Arabidopsis found high levels of methylation in gene-bodies (Zhang et al., 2006). Further, it was found that single-copy genes in soybean and common bean were methylated more frequently than duplicated genes (Kim et al., 2015). This suggests that gene-body methylation may have an important evolutionary role, possibly in protecting genes from mutation. DNA methylation has been implicated in physiological and developmental processes in plants. A study conducted in cotton suggested the role of CHH methylation in time-of-day and time-of-year memory (Jin et al., 2013). The fluctuation in DNA methylation pattern in response to different time intervals (Jin et al., 2013) and DNA methylation diversity between genotypes have been proposed in cotton (Osabe et al., 2014). Similarly, we presume differences in DNA methylation found between G19833 (Kim et al., 2015) and Sierra (current study) in common bean is partially due to the fact that two genotypes are derived from different gene pools.

Comparison of BS-Seq and MeDIP-Seq

As BS-Seq and MeDIP-Seq employ different mechanisms for determining methylated cytosines, the bioinformatic output of each method is also different. In the current study, BS-Seq (∼214 million reads) generated approximately 10-fold more reads than that of the immunoprecipitated reads via MeDIP-Seq (∼22.5 million reads). The cost per sample was about a twofold difference, BS-Seq being the more expensive option. Due to the expense, targeted bisulfite conversion, and sequencing is often chosen for specific genes or regions. However, this option gives much less information than both high-throughput sequencing technologies MeDIP-Seq and BS-Seq.

Since MeDIP-seq is affinity-based and the output is derived from peak shapes, the context of the methylated cytosine is more complicated to determine and less informative. BS-Seq yields individual nucleotide output, therefore, the specific context that is methylated is much more obvious. Because of this difference, precise trends about methylation in specific contexts cannot be made from MeDIP-Seq data. Methylation trends in annotated regions can be examined; for example Table 2 shows the total peaks, promoter, and genic peaks per chromosome. Genome viewers can be used to compare multiple datasets aligned to the same reference genome, we compared BS-Seq and MeDIP-Seq (Figures 4A,B), but other datasets such as RNA-Seq and ChIP-Seq could also be viewed simultaneously (Thorvaldsdottir et al., 2013).

FIGURE 4

Since researchers are often interested in finding genome-wide differences between a control and treatment, the use of either method would be useful. If only large changes in DNA methylation patterns are of interest, MeDIP-Seq may be the better option, since it is less expensive and uses an affinity-based approach. However, if determining small changes or very specific changes in DNA methylation are desired, BS-Seq is the better option, as it provides much more data. Further, the use of targeted BS-Seq would be the very least expensive per sample. With decreasing sequencing costs, in the future BS-Seq will likely be the best option, as long as there are means to store the data.

Conclusion

In summary, this is the first genome-wide DNA methylation profiling study in the Meso-American common bean cultivar Sierra using NGS approaches such as BS-Seq and MeDIP-Seq. Methylation statuses of each of the three DNA methylation contexts CG (∼50%), CHG (∼30), and CHH (∼2.7%) were determined from BS-Seq. Overall promoter and genic methylation trends were established. CG was the most commonly methylated context, and was more common in promoters than in genes. The opposite trend was found in CHG and CHH methylation contexts. The majority (93.2%) of MeDIP-Seq peaks were found in intragenic regions.

Materials and Methods

Plant Material

The common bean (P. vulgaris) cultivar “Sierra,” a pinto bean derived from the Meso-American gene pool was used in this study. Seeds were germinated for 2–3 days, on wet filter paper in petridishes and then planted in 6″ pots, filled with Promix Bx mycorrhizae soil, and grown in a greenhouse, at Delaware State University, Dover, DE, USA. Three separate times, three plants each were grown in the greenhouse under standard conditions, with a photoperiod of approximately 14 h days/10 h nights at 28/20°C (Ayyappan et al., 2015). At 14 days after planting, leaves from the three plants were harvested separately, high quality DNA was extracted, and equal amounts of DNA were combined to represent a single sample from these three plants, called Pool 1. Similarly, the second set of three plants, and the third set of three plants were processed to obtain Pool 1, Pool 2, and Pool 3. Finally, equal amounts of DNA from pools 1, 2, and 3 were combined for the BS-seq experiment. The same sample that was used for BS-seq was also used for the MeDIP-seq experiment so as to help us obtain a direct comparison of these two methods and their efficiency in identifying methylated regions in the genome. Leaves from each plant were collected after two weeks, flash frozen in liquid nitrogen, and stored in a –80°C freezer until DNA was extracted.

DNA Extraction

DNA was extracted from 2-week old leaves using a CTAB-based protocol as previously described (Doyle, 1990). The quantity and quality of DNA were determined by 1.0% agarose gel electrophoresis and Nanodrop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA). The experiment included a single pooled DNA sample collected from the leaves of nine different plants (biological replicates) and utilized the same sample for generating the MeDIP-Seq and BS-Seq libraries to maintain uniform experimental conditions. Both the libraries were sequenced on the Illumina/HiSeq-2500 platform.

Methylated Cytosine Sequencing

Methyl-MaxiSeqTM Library Construction

Methyl-MaxiSeqTM EpiQuest libraries were prepared from 100 ng of genomic DNA, which was bisulfite-treated using Zymo Research EZ DNA Methylation - LightningTM kit (Cat#: D5030, Zymo Research, Irvine, CA, USA). The bisulfite conversion rate was determined to be >99%. The bisulfite-converted DNA was subjected to a series of amplifications, which include a primer that contained part of the adapter sequence and four random nucleotides, followed by adding the remaining adapter sequence and barcoding the fragments, respectively. All PCR products were purified using the DNA Clean & Concentrator-5TM (Cat#: D4003, Zymo Research, Irvine, CA, USA). Library fragment size and concentration was checked using the Agilent 2200 TapeStation instrument and sequenced on the Illumina HiSeq 2500 platform (Illumina Inc., San Diego, CA, USA).

Methyl-MaxiSeqTM Sequence Alignments and Data Analysis

Sequence reads from bisulfite-treated EpiQuest libraries were identified using standard Illumina base-calling software and then analyzed using a Zymo Research proprietary analysis pipeline written in Python and using Bismark (http://www.bioinformatics.babraham.ac.uk/projects/bismark/) as the alignment software for analysis (Krueger and Andrews, 2011). Index files were constructed by bismark_genome_preparation command using the entire reference genome. Non_directional and all other default parameters were applied while running Bismark. The methylation level of each sampled cytosine was estimated as the number of reads reporting a C, divided by the total number of reads reporting a C or T. Promoter and gene body annotations were added using P. vulgaris genome annotations available at Phytozome (V1.0).

MeDIP-Seq Library Construction

Libraries for MeDIP-Seq were prepared following immunoprecipitation using the Zymo Research DNA Methylation IP Kit (Cat #D5101, Zymo Research, Irvine, CA, USA). Immunoprecipitated DNA was subjected to a series of amplifications as described above for preparing bisulfite-converted libraries. All PCR products were purified using the Zymo Research DNA Clean & Concentrator-5TM (Cat#: D4003, Zymo Research, Irvine, CA, USA). The input DNA library was prepared from pooled sample DNA that was fragmented and denatured. Libraries were quantified using the Agilent 2200 TapeStation and by qPCR. Sample concentrations were normalized to 4 nM, and then sequenced on the Illumina HiSeq 2500 platform (Illumina Inc., San Diego, CA, USA).

MeDIP-Seq Sequence Alignments and Data Analysis

Sequencing reads were aligned to the reference genome (V1.0) by Bowtie using best mode and other default parameters. Peak calling was done by “MACS2 callpeak” using input DNA as a control. BIGWIG files were generated from the coverage for visualization purposes (Zhang et al., 2008). Sequences derived from MeDIP-seq BS-seq libraries were submitted to the short read archives at NCBI BioProject#PRJNA306503.

Statements

Author contributions

MC conceived and performed experiments, interpreted the data, and developed the seminal draft of the manuscript. VS performed bioinformatic analysis, interpreted the data, and contributed significantly in drafting the manuscript. KH supervised and contributed to the drafting of the manuscript. VK conceived, supervised the whole project, and contributed to the drafting and editing of the manuscript.

Funding

We sincerely acknowledge funding from the National Science Foundation EPSCoR Grant No. IIA-1301765 and the State of Delaware (VK).

Acknowledgments

The authors acknowledge the assistance of Rita K. Hayford and other members of the Molecular Genetics and Epigenomics Laboratory at Delaware State University in editing and preparing this manuscript.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer JR and handling Editor declared their shared affiliation, and the handling Editor states that the process nevertheless met the standards of a fair and objective review.

Supplementary material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2016.00447

FIGURE S1 | The peak shapes for the input and immunoprecipitated DNA samples are shown.FIGURE S2 | Read coverage across the genome for each methylation context is shown. The top is CG, middle is CHG, and bottom is CHH contexts. Coverage read ranges from 0 (no coverage), to 1–4 reads, 5–10 reads, 11–20 reads, and 20 or more reads.

References

  • 1

    AkondG. M.CrawfordH.BertholdJ.TaulkderZ. I.HossainK. (2011). Minerals Zn, Fe, Ca, and Mg and antinutirent phytic acid constitutents in common bean.Am. J. Food Technol.6235243. 10.3923/ajft.2011.235.243

  • 2

    AlbertV. A.BarbazukW. B.DePamphilisC. W.DerJ. P.Leebens-MackJ.MaH.et al (2013). The amborella genome and the evolution of flowering plants.Science34212410891241089. 10.1126/science.1241089

  • 3

    AyyappanV.KalavacharlaV.ThimmapuramJ.BhideK. P.SripathiV. R.SmolinskiT. G.et al (2015). Genome-wide profiling of histone modifications (H3K9me2 and H4K12ac) and gene expression in rust (Uromyces appendiculatus) inoculated common bean (Phaseolus vulgaris L.).PLoS ONE10:e0132176. 10.1371/journal.pone.0132176

  • 4

    BeebeS. E.RaoI. M.BlairM. W.Acosta-GallegosJ. A. (2013). Phenotyping common beans for adaptation to drought.Front. Physiol.4:35. 10.3389/fphys.2013.00035

  • 5

    BitocchiE.NanniL. (2012). Mesoamerican origin of the common bean (Phaseolus vulgaris L.) is revealed by sequence data.Proc. Natl. Acad. Sci. U.S.A.109788796.

  • 6

    BlairM. W. (2013). Mineral biofortification strategies for food staples: the example of common bean.J. Agric. Food Chem.6182878294. 10.1021/jf400774y

  • 7

    BoykoA.KathiriaP.ZempF. J.YaoY.PogribnyI.KovalchukI. (2007). Transgenerational changes in the genome stability and methylation in pathogen-infected plants: (Virus-induced plant genome instability).Nucleic Acids Res.3517141725. 10.1093/nar/gkm029

  • 8

    BrinkmanA. B.GuH.BartelsS. J. J.ZhangY.MatareseF.SimmerF.et al (2012). Sequential ChIP-bisulfite sequencing enables direct genome-scale investigation of chromatin and DNA methylation cross-talk.Genome Res.2211281138. 10.1101/gr.133728.111

  • 9

    CâmaraC.UrreaC.SchlegelV. (2013). Pinto beans (Phaseolus vulgaris L.) as a functional food: Implications on human health.Agriculture390111. 10.3390/agriculture3010090

  • 10

    CaoX.AufsatzW.ZilbermanD.MetteM. F.HuangM. S.MatzkeM.et al (2003). Role of the DRM and CMT3 methyltransferases in RNA-directed DNA methylation.Curr. Biol.1322122217. 10.1016/j.cub.2003.11.052

  • 11

    ChodavarapuR. K.FengS.DingB.SimonS. A.LopezD.JiaY.et al (2012). Transcriptome and methylome interactions in rice hybrids.Proc. Natl. Acad. Sci. U.S.A.1091204012045. 10.1073/pnas.1209297109

  • 12

    CokusS. J.FengS.ZhangX.ChenZ.MerrimanB.HaudenschildC. D.et al (2008). Shotgun bisulphite sequencing of the Arabidopsis genome reveals DNA methylation patterning.Nature452215219. 10.1038/nature06745

  • 13

    DoyleJ. J. (1990). Isolation of plant DNA from fresh tissue.Focus121315.

  • 14

    EichtenS. R.BriskineR.SongJ.LiQ.Swanson-WagnerR.HermansonP. J.et al (2013). Epigenetic and genetic influences on DNA methylation variation in maize populations.Plant Cell2527832797. 10.1105/tpc.113.114793

  • 15

    FengS.CokusS. J.ZhangX.ChenP.-Y.BostickM.GollM. G.et al (2010). Conservation and divergence of methylation patterning in plants and animals.Proc. Natl. Acad. Sci. U.S.A.10786898694. 10.1073/pnas.1002720107

  • 16

    GentJ. I.EllisN. A.GuoL.HarkessA. E.YaoY.ZhangX.et al (2013). CHH islands: de novo DNA methylation in near-gene chromatin regulation in maize.Genome Res.23628637. 10.1101/gr.146985.112.as

  • 17

    GeptsP.AragãoF. J. L.BarrosE.De BlairM. W.BrondaniR.BroughtonW.et al (2008). Genomics of Phaseolus beans, a major source of dietary protein and micronutrients in the tropics.Genomics1113143. 10.1007/978-0-387-71219-2_5

  • 18

    GreavesI. K.GroszmannM.YingH.TaylorJ. M.PeacockW. J.DennisE. S. (2012). Trans chromosomal methylation in Arabidopsis hybrids.Proc. Natl. Acad. Sci. U.S.A.10935703575. 10.1073/pnas.1201043109

  • 19

    HendersonI. R.ChanS. R.CaoX.JohnsonL.JacobsenS. E. (2010). Accurate sodium bisulfite sequencing in plants.Epigenetics54749. 10.4161/epi.5.1.10560

  • 20

    HseihT.IbarraC. A.SilvaP.ZemachA.Eshed-WilliamsL.FischerR. L.et al (2009). Genome-wide demethylation of Arabidopsis endosperm.Science32414511454. 10.1126/science.1172417

  • 21

    HughesT.WebbR.FeiY.WrenJ. D.SawalhaA. H. (2010). DNA methylome in human CD4+ T cells identifies transcriptionally repressive and non-repressive methylation peaks.Genes Immun.11554560. 10.1038/gene.2010.24

  • 22

    JinX.PangY.JiaF.XiaoG.LiQ.ZhuY. (2013). A potential role for CHH DNA methylation in cotton fiber growth patterns.PLoS ONE8:e60547. 10.1371/journal.pone.0060547

  • 23

    KimK. D.El BaidouriM.AbernathyB.Iwata-OtsuboA.ChavarroC.GonzalesM.et al (2015). A comparative epigenomic analysis of polyploidy-derived genes in soybean and common bean.Plant Physiol.16814331447. 10.1104/pp.15.00408

  • 24

    KimK. D.El BaidouriM.JacksonS. A. (2014). Accessing epigenetic variation in the plant methylome.Brief. Funct. Genomics13318327. 10.1093/bfgp/elu003

  • 25

    KruegerF.AndrewsS. R. (2011). Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications.Bioinformatics2715711572. 10.1093/bioinformatics/btr167

  • 26

    LairdP. W. (2010). Principles and challenges of genome-wide DNA methylation analysis.Nat. Rev. Genet.11:191. 10.1038/nrg2732

  • 27

    LawJ. A.JacobsenS. E. (2010). Establishing, maintaining and modifying DNA methylation patterns in plants and animals.Nat. Rev. Genet.11204220. 10.1038/nrg2719

  • 28

    LiX.WangX.HeK.MaY.SuN.HeH.et al (2008). High-resolution mapping of epigenetic modifications of the rice genome uncovers interplay between DNA methylation, histone methylation, and gene expression.Plant Cell20259276. 10.1105/tpc.107.056879

  • 29

    LiX.ZhuJ.HuF.GeS.YeM.XiangH.et al (2012). Single-base resolution maps of cultivated and wild rice methylomes and regulatory roles of DNA methylation in plant gene expression.BMC Genomics13:300. 10.1186/1471-2164-13-300

  • 30

    ListerR.EckerJ. R. (2009). Finding the fifth base: genome-wide sequencing of cytosine methylation.Genome Res.19959966. 10.1101/gr.083451.108

  • 31

    ListerR.O’MalleyR. C.Tonti-FilippiniJ.GregoryB. D.BerryC. C.MillarA. H.et al (2009). Highly integrated single-base resolution maps of the epigenome in Arabidopsis.Cell133523536. 10.1016/j.cell.2008.03.029

  • 32

    LouS.LeeH.-M.QinH.LiJ.-W.GaoZ.LiuX.et al (2014). Whole-genome bisulfite sequencing of multiple individuals reveals complementary roles of promoter and gene body methylation in transcriptional regulation.Genome Biol.15:408. 10.1186/PREACCEPT-1031081530108509

  • 33

    MadlungA.ComaiL. (2004). The effect of stress on genome regulation and structure.Ann. Bot.94481495. 10.1093/aob/mch172

  • 34

    MedvedevaY. A.KhamisA. M.KulakovskiyI. V.Ba-AlawiW.BhuyanM. S. I.KawajiH.et al (2014). Effects of cytosine methylation on transcription factor binding sites.BMC Genomics15:119. 10.1186/1471-2164-15-119

  • 35

    MitchellD. C.LawrenceF. R.HartmanT. J.CurranJ. M. (2009). Consumption of dry beans, peas, and lentils could improve diet quality in the US population.J. Am. Diet. Assoc.109909913. 10.1016/j.jada.2009.02.029

  • 36

    MiuraA.NakamuraM.InagakiS.KobayashiA.SazeH.KakutaniT. (2009). An Arabidopsis jmjC domain protein protects transcribed genes from DNA methylation at CHG sites.EMBO J.2810781086. 10.1038/emboj.2009.59

  • 37

    OsabeK.ClementJ. D.BedonF.PettolinoF. A.ZiolkowskiL.LlewellynD. J.et al (2014). Genetic and DNA methylation changes in cotton (Gossypium) genotypes and tissues.PLoS ONE9:e86049. 10.1371/journal.pone.0086049

  • 38

    PetryN.BoyE.WirthJ.HurrellR. (2015). Review: The potential of the common bean (Phaseolus vulgaris) as a vehicle for iron biofortification.Nutrients711441173. 10.3390/nu7021144

  • 39

    RegulskiM.LuZ.KendallJ.DonoghueM. T. A.ReindersJ.LlacaV.et al (2013). The maize methylome influences mRNA splice sites and reveals widespread paramutation-like switches guided by small RNA.Genome Res.2316511662. 10.1101/gr.153510.112

  • 40

    SazeH.TsuganeK.KannoT.NishimuraT. (2012). DNA methylation in plants: relationship to small RNAs and histone modifications, and functions in transposon inactivation.Plant Cell Physiol.53766784. 10.1093/pcp/pcs008

  • 41

    SchmitzR. J.HeY.Valdés-lópezO.ResG.GentJ. I.EllisN. A.et al (2013). Epigenome-wide inheritance of cytosine methylation variants in a recombinant inbred population.Genome Res.2316631674. 10.1101/gr.152538.112

  • 42

    SchmutzJ.McCleanP. E.MamidiS.WuG. A.CannonS. B.GrimwoodJ.et al (2014). A reference genome for common bean and genome-wide analysis of dual domestications.Nat. Genet.46707713. 10.1038/ng.3008

  • 43

    SeymourD. K.KoenigD.HagmannJ.BeckerC.WeigelD. (2014). Evolution of DNA methylation patterns in the Brassicaceae is driven by differences in genome organization.PLoS Genet.10:e1004785. 10.1371/journal.pgen.1004785

  • 44

    ShenH.HeH.LiJ.ChenW.WangX.GuoL.et al (2012). Genome-wide analysis of DNA methylation and gene expression changes in two Arabidopsis ecotypes and their reciprocal hybrids.Plant Cell24875892. 10.1105/tpc.111.094870

  • 45

    SongQ.-X.LuX.LiQ.-T.ChenH.HuX.-Y.MaB.et al (2013). Genome-wide analysis of DNA methylation in soybean.Mol. Plant619611974. 10.1093/mp/sst123

  • 46

    StathamA. L.RobinsonM. D.SongJ. Z.CoolenM. W.StirzakerC.ClarkS. J. (2012). Bisulfite sequencing of chromatin immunoprecipitated DNA (BisChIP-seq) directly informs methylation status of histone-modified DNA.Genome Res.2211201127. 10.1101/gr.132076.111

  • 47

    StevensM.ChengJ. B.LiD.XieM.HongC.MaireC. L.et al (2013). Estimating absolute methylation levels at single-CpG resolution from methylation enrichment and restriction enzyme sequencing methods.Genome Res.2315411553. 10.1101/gr.152231.112

  • 48

    StroudH.DingB.SimonS. A.FengS.BellizziM.PellegriniM.et al (2013a). Plants regenerated from tissue culture contain stable epigenome changes in rice.Elife2114. 10.7554/eLife.00354

  • 49

    StroudH.DoT.DuJ.ZhongX.FengS.JohnsonL.et al (2013b). Non-CG methylation patterns shape the epigenetic landscape in Arabidopsis.Nat. Struct. Mol. Biol.216472. 10.1038/nsmb.2735

  • 50

    TaiwoO.WilsonG. A.MorrisT.SeisenbergerS.ReikW.PearceD.et al (2012). Methylome analysis using MeDIP-seq with low DNA concentrations.Nat. Protoc.7617636. 10.1038/nprot.2012.012

  • 51

    TakunoS.GautB. S. (2012). Body-methylated genes in Arabidopsis thaliana are functionally important and evolve slowly.Mol. Biol. Evol.29219227. 10.1093/molbev/msr188

  • 52

    TaylorK. H.KramerR. S.DavisJ. W.GuoJ.DuffD. J.XuD.et al (2007). Ultradeep bisulfite sequencing analysis of DNA methylation patterns in multiple gene promoters by 454 sequencing.Cancer Res6785118518. 10.1158/0008-5472.CAN-07-1016

  • 53

    ThorvaldsdottirH.RobinsonJ. T.MesirovJ. P. (2013). Integrative genomics viewer (IGV): high-performance genomics data visualization and exploration.Brief. Bioinform.14178192. 10.1093/bib/bbs017

  • 54

    VerhoevenK. J. F.JansenJ. J.van DijkP. J.BiereA. (2010). Stress-induced DNA methylation changes and their heritability in asexual dandelions.New Phytol.18511081118. 10.1111/j.1469-8137.2009.03121.x

  • 55

    WangJ.MarowskyN. C.FanC. (2014). Divergence of gene body DNA methylation and evolution of plant duplicate genes.PLoS ONE9:e110357. 10.1371/journal.pone.0110357

  • 56

    WestP. T.LiQ.JiL.EichtenS. R.SongJ.VaughnM. W.et al (2014). Genomic distribution of H3K9me2 and DNA methylation in a maize genome.PLoS ONE9:e105267. 10.1371/journal.pone.0105267

  • 57

    ZemachA.KimM. Y.SilvaP.RodriguesJ. A.DotsonB.BrooksM. D.et al (2010). Local DNA hypomethylation activates genes in rice endosperm.Proc. Natl. Acad. Sci. U.S.A.1071872918734. 10.1073/pnas.1009695107

  • 58

    ZhangX.YazakiJ.SundaresanA.CokusS.ChanS. W.-L.ChenH.et al (2006). Genome-wide high-resolution mapping and functional analysis of DNA methylation in Arabidopsis.Cell12611891201. 10.1016/j.cell.2006.08.003

  • 59

    ZhangY.LiuT.MeyerC. A.EeckhouteJ.JohnsonD. S.BernsteinB. E.et al (2008). Model-based Analysis of ChIP-Seq (MACS).Genome Biol.9:R137. 10.1186/gb-2008-9-9-r137

Summary

Keywords

Phaseolus vulgaris, common bean, methylome, whole genome bisulfite-sequencing, MeDIP-Seq

Citation

Crampton M, Sripathi VR, Hossain K and Kalavacharla V (2016) Analyses of Methylomes Derived from Meso-American Common Bean (Phaseolus vulgaris L.) Using MeDIP-Seq and Whole Genome Sodium Bisulfite-Sequencing. Front. Plant Sci. 7:447. doi: 10.3389/fpls.2016.00447

Received

29 October 2015

Accepted

22 March 2016

Published

26 April 2016

Volume

7 - 2016

Edited by

Oswaldo Valdes-Lopez, Universidad Nacional Autónoma de México, Mexico

Reviewed by

Norberto Daniel Iusem, Universidad de Buenos Aires, Argentina; Jose Luis Reyes, Universidad Nacional Autónoma de México, Mexico

Updates

Copyright

*Correspondence: Venu Kalavacharla,

This article was submitted to Plant Nutrition, a section of the journal Frontiers in Plant Science

Disclaimer

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Article metrics