Single-base resolution maps of cultivated and wild rice methylomes and regulatory roles of DNA methylation in plant gene expression

1DNA methylation is an epigenetic modification mechanism that plays essential roles in diverse biological processes. It has also been proposed to be an alternative inheritance system playing an important role in evolution, as many case studies in plants and animals have revealed that differentially methylated alleles could create heritable phenotypic changes across generations, including some agronomically important traits in rice. Recently, single-base resolution methylome maps of a dicot plant (Arabidopsis thaliana), human and silkworm have been successfully generated by whole-genome sequencing bisulfite-treated genomic DNA using next-generation sequencing technology (BS-Seq), which revealed more elaborate patterns and functional effects of DNA methylation at the whole-genome level.

2Rice is not only one of the most important crops as the primary food source for more than half of the world’s population, but also an important model system for the evolutionary study of cereals and the molecular study of monocot plants. DNA methylation serves various important functions and thus has been of great interest to rice geneticists and breeders. Pioneer studies of epigenetic modifications in rice, including DNA and histone methylation, using traditional methylated DNA enrichment method suggested possible functional roles of DNA methylation in rice, but this approach is difficult to discriminate major genomic elements including promoters, gene bodies, transposons, and repeats. Recently two other studies comparing methylation patterns among many species using BS-Seq technology briefly reported genome methylation patterns for japonica rice strain Nipponbare. However, the two Nipponbare methylomes had relatively low sequencing coverage (< 4× per base for each strand). Furthermore, one of the studies used different tissues to obtain methylomic (from leaves) and transcriptomic (from shoots) profiles and the other did not obtain gene expression data at all, which made it difficult to accurately analyze the regulatory effects of DNA methylation in rice. As acknowledged by Feng et al., such a low genome-wide sequencing depth permitted good assessment of the level of methylation of major genomic elements including genes, transposons, and repeats, but was not sufficient for quantifing the methylation level at individual cytosines. Furthermore, to what extent and how the cultivated rice has evolved divergent DNA methylation pattern from its wild relative species still need to be addressed.

3In this study, we generated single-base-resolution DNA methylomic maps as well as transcriptomic profiles for young panicles of the two Asian cultivated rice subspecies, Oryza sativa ssp. japonica and O. sativa ssp. indica, and their wild relatives, Oryza rufipogon and Oryza nivara. The panicle is an important organ showing strong differentiation between cultivated and wild rice and directly affects the major yield components including the number of spikelets and the percentage of filled grains. The high-resolution DNA methylomes of cultivated and domestic rice will not only serve as references for future molecular studies of rice epigenetics but also shed new lights into epigenetic mechanisms of plant domestication.

4To investigate the general methylation patterns of rice as well as the DNA methylation divergence between cultivated and wild rice, we included in our samples both subspecies of Asian cultivated rice, Oryza sativa spp. japonica (represented by Dianjingyou1, bred by the Yunnan Academy of Agricultural Sciences, China, which is a typical japonica rice type mainly suitable to be planted in Yunnan) and indica (IR64, from the International Rice Research Institute, IRRI), and their wild relatives, Oryza rufipogon (Accession 105327, originally collected from Sri Lanka and provided by IRRI) and Oryza nivara (Accession 105426, originally collected from India and provided by IRRI).

5For each of the four rice lines, a young panicle was used to generate a methylomic map with the BS-Seq method. 170–320 million sequencing reads were generated for each of the four samples, respectively. After removing low-quality and clonal reads (artificially generated during amplification of bisulfate-treated DNA in constructing sequencing libraries, see Material and methods for details), 58–176 million uniquely mapped high-quality reads were retained for each of the four rice lines, yielding 2.6-9.2 gigabases (Gb) of data which cover 76-91% of the reference Nipponbare genome (International Rice Genome Sequencing Project, 2005), respectively (Table 1). The read depths range from 4.5× to 13.5× per base for each DNA strand (Table 1). Because japonica has the best reference genome and we obtained the highest genome coverage and sequencing depth for the japonica strain Dianjingyou1, here we use the japonica accession as the representative to describe the general methylation landscape of rice, with the differences between wild and cultivated rice being discussed later.

6We used the unmethylated chloroplast genome to calculate the sum of non-conversion rate and T-C sequencing error rate, which is as low as 0.47% to 1.17% for the four samples respectively (Table 2). Using these values we then conducted binomial tests with false positive rate below 5% to exclude those mCs that may be the results of non-conversion of cytosines during our bisulfite treatment or T to C sequencing errors during the base calling process. After this correction, we identified 35,598,491 methylated cytosines (mCs) accounting for 24.3% of all covered cytosines throughout the reference genome in the japonica Dianjingyou1 (Table 2), which is four times higher than that of Arabidopsis in terms of mC density. The percentages of mCs in CG, CHG (with H being A, C or T) and CHH contexts are 54.7%, 37.3% and 12.0%, respectively while the average methylation level in the three contexts are 44.5%, 24.1% and 4.7%, respectively, with methylation level being defined as the proportion of reads showing mC among all reads covering the same cytosine site. Both of these measures reveal that rice has much higher level of genome-wide DNA methylation than A. thaliana (Table 2). These patterns are more or less the same in indica and the two wild rice species (Table 2), and are consistent with the two previous studies of rice leaf methylomes.

7The genome-wide methylation patterns with respect to the genomic structure in rice young panicles are similar to those observed in A. thaliana and in rice leaves, particularly in the relative prevalence of mCs in the contexts of GG, CHG, and CHH (Figure 1a), the tendency toward hypermethylation in CG context but hypomethylation in CHH (Figure 1b), the high methylation level in transposable elements (TEs) and the relatively low level of methylation in genic regions, the enrichment of CG methylation in gene bodies (Figure 2a-f), and the enrichment of small RNA loci in TEs and depletion in genic regions (Figure 2g and 2h).

8However, in contrast to the significantly enriched methylation in all the three sequence contexts around the centromeric regions in A. thaliana genome, rice shows much less difference between centromeric and non-centromeric regions with only slight CG and CHG methylation enrichment around the centromeric regions and almost uniform distribution of CHH mCs across entire chromosomes (Figure 3). This is probably due to the fact that compared to A. thaliana, rice has significantly larger amount of pericentromeric heterochromatin and much higher proportion of TEs distributed across the genome, both of which are heavily methylated.

9In addition, we found a positive correlation between sequence length and methylation density for genes (Figure 2i), but not for TEs (Figure 2j), which is different from the case in A. thaliana wherein both genes and TEs showed positive correlation between sequence length and methylation. This might be due to the fact that TEs in rice genome are almost saturated with mCs regardless of their length. All the above results are consistently found in all the four rice samples (see Additional files 1 and 2 for results of other samples), showing these are general patterns in both cultivated and wild rice. These differences between rice and Arabidopsis suggest that differences in genomic organization and TE composition among plants could result in different epigenomic landscapes.

10To assess the effects of rice DNA methylation on gene expression, we also generated genome-wide gene expression profiles for all the four samples (Table 3) using Digital Gene Expression tag profiling (DGE) technology, which combines the classic SAGE (Serial Analysis of Gene Expression) with the Illumina sequencing techniques. In total we obtained 7,343,629 to 7,774,577 raw reads with 599,304 to 648,947 unique tag sequences for each of the four samples. After filtration and mapping reads to 24,955 non-redundant full-length cDNAs (FL-cDNAs) of rice, our DGE data cover 79% to 83% of all FL-cDNAs with a one-tag-one-cDNA relationship (Table 3).

11Although some previous studies have also explored the relationship between DNA methylation and gene expression in rice, their low resolution of genome-wide methylated cytokines (mCs) or utilization of different tissues in generating transcriptomic and methylomic profiles created the need for more elaborate methylomic studies to comprehensively unveil high resolution rice methylomes and detailed functional effects of rice DNA methylation. Here we shall mainly use the japonica data to present the detailed patterns in rice, which are similar in all four samples (see Additional file 3 for results of the other three samples).

12Our results show that promoter-unmethylated genes have significantly higher expression level than promoter-methylated genes (p = 1.915e-08, Wilcoxon rank sum test), indicating that promoter methylation represses gene expression. Consistent with this conclusion, the repression effect is weak for slightly and moderately promoter-methylated genes while very strong for heavily methylated ones (Figure 4a, and see Additional file 4 for each sequence context). These results from the high resolution genome-wide data proved false the previous speculation that promoter methylation may not affect rice gene transcription but are consistent with the findings from A. thaliana, human, and the recent rice methylome study, confirming that promoter methylation is a general mechanism suppressing gene expression in eukaryotes.

13In contrast to promoter methylation, gene-body methylation generally appears to be positively correlated to gene expression as body-methylated genes have significantly higher expression level than body-unmethylated genes (p < 2.2e-16, Wilcoxon rank sum test). However, further analysis revealed complicated relationships between gene-body methylation and gene expression. At first gene expression levels increase with methylation levels, but after a certain point, heavy gene-body methylation appears to repress gene expression; and consequently genes with moderate levels of body methylation tend to have the highest expression levels (Figure 4b, see Additional file 4 for each sequence context). These observations are consistent with previous studies in A. thaliana and in rice. It has been proposed that gene-body methylation can prevent transcriptional initiation from cryptic sites within genes but at the cost of impeding transcriptional elongation. This trade-off may have led to the observation that moderately body-methylated genes have the highest level of expression.

14Most interestingly, our study revealed that methylation in the transcriptional termination region (TTR) is also highly correlated with gene expression. In a pattern similar to promoter methylation, the TTR-unmethylated genes are expressed at significantly higher level than TTR-methylated genes (p < 2.2e-16, Wilcoxon rank sum test). Moreover, an approximately monotonic negative correlation exists between TTR methylation and gene expression (Figure 4c, see Additional file 4 for each sequence context). Surprisingly, the correlation coefficient is even higher than that of promoters, especially for CG methylation (Figure 4d), suggesting that TTR methylation may play an even more important role in gene expression regulation than promoter methylation. To exclude the possibility that the negative correlation between TTR methylation and gene expression is an indirect effect caused by promoter methylation if promoter-methylated genes are also prone to have methylated TTRs, we repeated the correlation analysis for both TTR and promoter regions using genes without promoter methylation and those without TTR methylation respectively. The strong negative correlation between TTR methylation and gene expression and the higher correlation coefficients than that for promoter methylation could still be observed (Figure 4f and 4 g).

15DNA methylation, together with associated histone modification, has been suggested to influence the binding of RNA polymerase to DNA and thus affect the initiation, elongation and termination of the gene transcription process. It has been established that promoter methylation can repress gene expression and promoter hypomethylation may be required for genes to express efficiently. Gene-body methylation has also been proposed to inhibit transcriptional noise in actively transcribed genes and consequently body-methylated genes usually show moderate to high level of expression. It is plausible that TTR methylation could have significant effect on gene expression through interfering with transcriptional termination. Consistent with this hypothesis, active DNA demethylation mediated by the DEMETER (DME) family has been found primarily at both the 5’- and 3’- ends of genes in A. thaliana, suggesting a functional role of methylation in both regions.

16To further investigate whether such a regulatory mechanism is shared by plants, we also examined the effects of TTR methylation on gene expression in Arabidopsis using existing data. Consistent with the results from rice, our analysis revealed significant negative correlation between TTR methylation and gene expression in Arabidopsis (Figure 4e), suggesting a general regulatory role of TTR methylation both in monocot and dicot plants.

17In addition, we found that the 5’-end of gene coding region is another important regulatory region showing significant positive correlation between its methylation and gene expression in Arabidopsis, but not in rice (Figure 4d and 4e), consistent with the Arabidopsis-specific CHH methylation enrichment in this region. Whether this is a dicot-specific regulation mechanism needs further studies using more plant species of both monocots and dicots. It is worth noting that the positive correlation between the 5’-end gene coding region methylation and gene expression in Arabidopsis could be revealed based on both the absolute (total methylation level of mCs divided by sequence length) and the relative (total methylation level of mCs divided by total number of cytosine sites in a region) methylation levels, but the correspondent CHH methylation enrichment in Arabidopsis could only be revealed using the absolute methylation level (Additional file 5). This may explain the failure of previous studies on Arabidopsis to reveal the positive effect of 5’-end gene methylation on gene expression and suggests absolute methylation level may be better than relative methylation level to affect gene expression.

18To examine the divergence between cultivated and wild rice at genetic, methylation and gene expression levels, we first constructed trees based on genomic, methylomic, and gene expression data, respectively (see Materials and methods) (Figure 5a-c). The genomic DNA tree shows the same topology with the methylation-based trees except for the one constructed using CHH mCs. The consistency between the genomic and methylomic trees suggests that genetic divergence at DNA sequence level may be the major determinant of methylation patterns in the genome. We further examined the relationship between genetic and methylation divergence among samples using sliding windows across the whole genome. The average number of nucleotide differences per site and the average Spearman correlation coefficient of mCs were used to measure genetic and methylation divergence among samples, respectively (see Materials and methods). Through sliding window analysis, we found that regions with high genetic sequence divergence also have high CG methylation divergence (low correlation coefficient) among samples (Figure 6). However, this pattern was not obvious for non-CG methylation, especially CHG methylation. Taken together, these results suggested that genetic divergence may be the major determinant of CG methylation patterns which usually showed high methylation levels, while non-CG methylation level is usually low and may be easily affected by other internal/external factors besides the DNA sequence.

19However, the topology of the gene expression tree is different from both genomic DNA and methylation-based trees, with the two cultivated rice subspecies tightly clustering together and the two wild rice species being most similar to each other, a pattern consistent with the phenotypic relationships among the four samples. The genome-wide gene expression divergence of panicles among the four rice subspecies/species obviously departs from the expectation under the neutral evolution model that posits that gene expression differentiation positively correlates with species’ genetic divergence. Instead, it suggests rice domestication may have occurred through changes in a limited number of genes that have pleiotropic and/or cascading effects on gene expression at the whole-genome level. With considerable gene expressional and phenotypic divergence between wild and cultivated rice, the panicle might have been an major target for artificial selection during domestication, a hypothesis consistent with the finding that many yield-related traits are associated with panicles.

20The above results, i.e., the CHH methylation tree’s inconsistency with the genomic tree and the low correlation between non-CG methylation divergence and genetic divergence, imply that high-level methylation in CG context might be more conserved than low-level methylation in CHG and CHH contexts among species. To further test this hypothesis, we calculated the coefficient of variation (CV) using a sliding window of the methylation level for each of the CG, CHG and CHH context to measure their conservation levels among species. Our results show that CG methylation has the lowest variation (8.6%), followed by CHG (11.0%) and CHH (13.7%), consistent with the decrease of their methylation level in the same order. In addition, we calculated the correlation coefficient between methylation level and its CV among four samples in gene regions for each sequence context. We found that gene methylation level showed significant negative correlation with its CV regardless of sequence contexts (Figure 5d). These results suggested that high-level methylation states are more stable during evolution.

21To examine the methylation variation among species in different functional elements, we calculated the CVs of methylation level for mCs located in gene, promoter, TTR and TE regions (Additional file 6). We found that TEs’ methylation status is most conserved among species with the smallest CV, followed by genes, promoters and TTRs, consistent with the previous observation in Arabidopsis. The methylation conservation levels of different functional elements also follow their methylation levels, again suggesting that high-level methylation states are more stable among species. Further, we compared the methylation variations among different location and sequence contexts in the same methylation level groups. We found that mCs with similar methylation levels have similar variations regardless of their genomic location or sequence contexts (Additional file 7), demonstrating that methylation level is a consistent indicator of methylation variation among species.

22To identify the DNA methylation changes that may be associated with rice domestication, we identified in cultivated rice 14/24/49 methylation-upregulated and 21/10/46 methylation-downregulated genes in promoter, TTR and gene body regions respectively, leading to a total of 155 non-redundant differentially-methylated genes between cultivated and wild rice (Additional file 8). We also picked 6 promoter or TTR regions in total to validate such methylation differentiation between cultivated and wild rice using the traditional bisulfite sequencing method for single genes, and the results conform to those from our whole genome BS-Seq analyses (Additional files 910111213,and 14), indicating the reliability of our BS-Seq results. Among the 155 genes, 11 (7.1%) show methylation-correlated 2-fold gene expression changes, but the proportion of genes with such expressional changes among differentially methylated genes is not significantly different from the proportion of all genes showing 2-fold gene expression changes between cultivated and wild rice regardless whether there is DNA methylation difference. Interestingly, a similar conclusion has also been drawn among natural accessions of Arabidopsis thaliana in which only 6% of differentially methylated genes have significantly different expression levels between ecotypes, and the proportion of expression-altered genes is the same as that among all genes. The results from both rice and Arabidopsis suggest that a variety of mechanisms, including genetic changes in genes’ cis- or trans-regulators and chromatin modification, together with DNA methylation, regulate gene expression at the genomic level. However, it is also possible that subtle (less than two-fold) changes in gene expression caused by DNA methylation differences could have important consequences, particularly if those genes are of major or multiple effects. In the case of rice domestication, therefore, it still cannot be ruled out that a few key genes’ epigenetic and correlated expressional changes might have played important roles in the production of some important agronomic traits in cultivated rice. Given the small sample size in this study, methylomic and transcriptomic analysis of more representative wild and cultivated rice is needed to further clarify this important issue.

23The high resolution DNA methylation maps for the two cultivated rice subspecies and their wild ancestors obtained in this study and the integrated analysis of genomic, epigenomic, and transcriptomic data have not only broadened our understanding of the mechanisms of gene regulation and the complicated relationships between DNA divergence, DNA methylation, and gene expression variation in plant genomes, but also provided valuable data for future studies on rice epigenetics as well as on epigenetic differentiation between wild and cultivated rice.

24To make our methylomes from four samples comparable, we carefully collected all samples at the same developmental stage of panicle initiation to booting. This stage is an important time point, at which rice start transitioning from vegetative growth to reproductive growth.

25A single young panicle from each of the cultivated rice subspecies and the two wild rice species was ground in liquid nitrogen to fine powder using mortar and pestle. Genomic DNAs were isolated using the Plant Genomic DNA Purification Kit (Tiangen Inc., China) and total RNAs were isolated using the RNeasy Plant Mini Kit (Qiagen Inc., Germany). DNA was fragmented by sonication with the Diagenome sonicator to a mean size of approximately 250 bp, followed by blunting, 3’-end addition of dA, and adaptor ligation according to the manufacturer’s instruction (Illumina). The bisulfite conversion of rice DNA was carried out using a modified (NH4)HSO3-based protocol. Bisulfite-treated DNAs were PCR amplified with 16 cycles. The resultant DNAs were applied to paired-end sequencing with the read length of 44 or 75 nt for each end using the ultrahigh-throughput Illumina Genetic Analyzer (GA 2) as per the manufacturer’s instructions.

26Because japonica rice has high quality reference genome sequence and gene annotation information, all reads from four rice samples were mapped to the Nipponbare IRGSP genome sequence (build 4 assembly), which was downloaded from RAP-DB (http://rapdblegacy.dna.affrc.go.jp/archive/build4/OsGenome_RAP2.tar.gz). Since DNA methylation has strand specificity, the plus strand and the minus strand of Nipponbare genome should be separated and used as different alignment target sequences for BS-Seq reads. That is, each cytosine in reference genome sequences was converted to thymine, termed T-genome which represents the plus strand. Meanwhile, each guanine in reference genome sequences was converted to adenosine, termed A-genome which represents the minus strand. To map the raw 44 or 75 nt pair-ended BS-Seq reads, the original reads were computationally converted to the alignment forms with the following steps: 1) observed cytosines on the forward read of each read pair were in silico replaced by thymines; 2) observed guanines on the reverse read of each read pair were in silico replaced by adenosines. The converted reads were then mapped to both strands of the A- and T- genome sequence using SOAP2 allowing up to two mismatches for 44 nt reads and four mismatches for 75 nt reads. Reads mapped to the same start position for both ends were regarded as clonal duplicates, which might have been generated during PCR process, and only one of them was kept. Only reads uniquely mapped to either of the strands were then retained for further analysis. After above filtration, for methylcytosine (mC) detection, we transformed each aligned read and the two strands of the Nipponbare genome back to their original forms to build an alignment between the original forms. Cytosines in the BS-Seq reads matching the corresponding cytosines in the plus strand of the reference genome, or guanines in the BS-Seq reads matching the corresponding guanines in the minus strand of the reference genome will be regarded as potential mCs. To exclude the false positive caused by base calling process, we removed those potential mCs with Q scores lower than 20, which means that a base is correctly called at more than 99% probability, a highly conservative criterion for calling reliable bases.

27We used the unmethylated chloroplast genome to calculate the sum of non-conversion rate and T-C sequencing error rate, and then conducted binomial tests using these values with false positive rate below 5% to exclude those mCs that may be the results of non-conversion of cytosines during our bisulfite treatment or T to C sequencing errors during base calling process.

28The method of mapping BS-seq reads has been successfully applied in profiling the silk gland methylome of silkworm and the peripheral blood mononuclear cell methylome of humans. In these published studies as well as in this study, experimental validation confirmed the reliability of mC calling. To further compare the performance between our pipeline and other reported softwares, we chose data from one lane of japonica sample (17,406,765 paired reads) and the plus strand of chromosome01 to call mCs using both our method and Bismark with the same parameters (because BS seeker does not support paired-end reads, we could not use BS seeker to do the comparison). The results showed that in general the two methods gave very similar results. For example, 92%, 91% and 95% mCs for CG, CHG and CHH sequence context called by Bismark were also identified by our method. And the methylation levels obtained using two methods show high correlation for each sequence context (see Table 4 for details). Furthermore, our method shows slightly higher sensitivity than Bismark, which could be caused by employing different short reads aligners (our method used SOAP, while Bismark used Bowtie). Therefore, we used mCs called with our pipeline to conduct all following analyses.

29Because DNA methylation has strand specificity, BS-Seq reads mapped to T-genome contain the DNA sequence information of plus strand, while reads mapped to A-genome contain the information of minus strand. Thus we could use BS-Seq reads to call SNPs directly for both plus and minus strands for each sample. After mapping and processing of BS-Seq reads described in above step, SNPs in both strands were called using Samtools. To get high-quality SNPs, we used strict standards for our results—only bases with base quality (Q score) of at least 30 and positions covered by at least 5 reads were used for our SNP calling process. Although bisulfite treatment will convert C to T in both strands which will affect judgment of some kinds of genotypes during SNPs calling process, we can use SNP information from both strands to solve this problem. Some genotypes (A/G, G/A, C/A, C/G, T/A, T/G, A/A, G/G, the former indicates reference base and latter indicates sample’s base) can be correctly called using results from plus strand, while the other genotypes (A/C, A/T, G/T, G/C, C/T, T/C, T/T, C/C) can be correctly called using results from minus strand. Under this principle, a total of 37,141, 50,340, 65,680, and 83,594 high-quality SNPs in japonicaindicaO. rufipogon, and O. nivara were obtained respectively and used for further analysis.

30To verify the BS-Seq results, we picked 6 regions, including 4 promoter regions (Os12g0264800, Os01g0116800, Os01g0543000, and Os04g0431700 genes) and 2 TTR regions (Os11g0111101 and Os08g0150200 genes), to validate the methylation status in all the four samples using the traditional bisulfite sequencing method for single genes.

31DGE-tag libraries were constructed from the total RNAs isolated from the same four rice panicles used for extracting DNA with the DGE-Tag Profiling NlaIII Sample Prep Kit (Illumina) according to the manufacturer’s instructions.

32Full-length cDNA sequences of rice genes were downloaded from RAP-DB (http://rapdblegacy.dna.affrc.go.jp/archive/build4/rep_RAP2.tar.gz). A total of 24,955 genes supported with full length cDNA were used for DGE analysis. All possible CATG + 17nt tag sequences were created from 24,955 full-length cDNAs and used as the reference tag database. Unique tag sequences and their numbers were extracted from our raw DGE tags and these tags were aligned against the reference tag database using SOAP. Only perfect matches were kept for further analysis. Within genes, most of the DGE tags were mapped to the most 3’-end CATG sites of genes (Additional file 15), suggesting that transcriptional termination sites of most rice genes we used were annotated correctly and that performing DGE analysis for those genes was suitable and reliable. The expression level of a gene was represented by the total number of tags that uniquely aligned to that gene. Gene expression levels were normalized to tag number per million tags for gene expression comparisons among different samples.

33Annotation of TEs was downloaded from RAP-DB (http://rapdblegacy.dna.affrc.go.jp/archive/build4/OsNIAS_b4_chromOut.tar.gz) and smRNA sequences were downloaded from rice MPSS database (http://mpss.udel.edu/data-files/rice/small/smallRNA_summary.txt). Sequences of smRNAs were mapped to the rice reference genome using SOAP without allowing mismatch, and uniquely mapped smRNAs were used for further analysis. Methylation level refers to the proportion of reads showing mC among all reads covering the same cytosine site. It can further be classified as absolute methylation level (total methylation level of mCs divided by the total sequence length of the calculated region) and relative methylation level (total methylation level of mCs divided by total number of cytosine sites in the calculated region), both of which were used for our analysis. TE or smRNA density was defined as the ratio of the number of bases belonging to TEs or smRNAs to the total length of the calculated region.

34The single-base methylation profile and corresponding gene expression profile of Arabidopsis were downloaded from NCBI SRA (Sequence Read Archive) database (http://www.ncbi.nlm.nih.gov/sra, the accession numbers are SRA000284 for BS-Seq, also referred as MethylC-Seq in Lister et. al.’s paper, and SRA000286 for mRNA-seq). The genome sequences and gene annotation information (TAIR9) were downloaded from the ftp site of TAIR (The Arabidopsis Information Resource, ftp://ftp.arabidopsis.org/home/tair/). The method of mapping and processing of BS-Seq reads is the same as that used in rice. Mapping and processing of mRNA-seq reads were conducted using TopHat software. Gene expression levels, measured by reads per kilobase of transcript per million reads (RPKM), was also calculated using TopHat. Methylation-expression correlation analyses for Arabidopsis were performed with the same methods as used for rice.

35For methylation tree construction, we first calculated the absolute methylation level in sliding windows of 50 kb with steps of 25 kb across the whole genome for the four samples, and then clustered the samples based on pairwise Spearman correlation coefficients estimated from the above whole genome sliding methylation level matrix. Methylation levels of all mCs with ≥ 5 × coverage of the whole genome were also used to calculate the pairwise Spearman correlation coefficients among samples, which in turn produced similar methylation tree. High-quality SNPs called from BS-Seq reads were used to construct genomic tree among the four samples using p distance and neighbor-joining method implemented in MEGA. Finally, an expression tree was constructed using pairwise Spearman expression correlation coefficients among the four samples based on DGE data. Methylation and expression trees were constructed using the hclust function of R statistical software (http://www.r-project.org/).

36To examine the relationship between genetic and methylation divergence, the average number of nucleotide differences per site among samples for different genomic regions was calculated using SNPs obtained from the above step with 50 kb sliding window and a step of 25 kb through the whole genome, which was used to measure genetic divergence. Then average Spearman correlation coefficients of methylation level of all cytosines with ≥ 5 × coverage among samples were calculated for same sliding windows and used to measure methylation divergence in different genomic regions.

37Because mCs in promoters, TTRs and gene bodies have significant effects on gene expression, we identified genes with different methylation status in these three regions among the four rice samples. Only regions with above 80% sequencing coverage were used for further analysis. Methylation levels of all Cs within promoters/TTRs/gene-bodies were calculated and used to perform two-sample Wilcoxon rank sum tests between any samples. Differentially-methylated genes were identified for each pairwise comparison using a significance level of alpha < 0.05. Genes that are significantly methylation-upregulated or methylation-downregulated consistently in all cultivated vs. wild rice pairwise comparisons but are not significantly different within cultivated or wild rice comparisons were respectively identified as methylation-upregulated or methylation-downregulated genes between cultivated and wild rice.

38The methylome data have been deposited into the NCBI Short Read Archive (SRA, http://www.ncbi.nlm.nih.gov/sra/) under accession number SRA012190 and DGE data have been deposited into the NCBI Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo) under accession number GSE20871.

39TTRs, transcriptional termination regions; BS-Seq, bisulfite-treated genomic DNA sequencing; mC, methylated cytosine; TEs, transposable elements; DGE, Digital Gene Expression tag profiling; SAGE, Serial Analysis of Gene Expression; FL-cDNAs, nonredundant full-length cDNAs; CV, coefficient of variation.

40The authors declare that they have no competing interests.

41WW, JW, JZ and SG designed the study. XL, WW, HX, XX, XC and JW wrote the manuscript. FH, XZ, SG, SZ and QL collected the samples. JZ, JY, XZ and HX constructed the BS-Seq libraries and conducted the BS validation. MY, GZ, RL, CY performed the Illumina sequencing. XL, HX, FH, and XX performed the data analyses. All authors read and approved the final manuscript.