Whole-genome DNA methylation patterns of Oryza sativa (L.) and Oryza glumaepatula (Steud) genotypes associated with aluminum response.

Epigenetic mechanisms in crops have emerged as a fundamental factor in plant adaptation and acclimation to biotic and abiotic stresses. Among described epigenetic mechanisms, DNA methylation has been defined as the most studied epigenetic modification involved in several developmental processes. It has been shown that contrasting methylation marks are associated with gene expression variations between cultivated and wild crop species. In this study, we analyzed single-base resolution methylome maps for Oryza sativa (a cultivated species) and Oryza glumaepatula (a wild species) genotypes grown under control conditions. Our results showed that overall, genome-wide methylation profiles are mainly conserved between both species, nevertheless, there are several differentially methylated regions with species-specific methylation patterns. In addition, we analyzed the association of identified DNA methylation marks in relation with Aluminum-tolerance levels of studied genotypes. We found several differentially methylated regions (DMRs) and DMR-associated genes (DAGs) that are linked with Al tolerance. Some of these DAGs have been previously reported as differentially expressed under Al exposure in O. sativa. Complementarily a Transposable Elements (TE) analysis revealed that specific aluminum related genes have associated-TEs potentially regulated by DNA methylation. Interestingly, the DMRs and DAGs between Al-tolerant and susceptible genotypes were different between O. sativa and O. glumaepatula, suggesting that methylation patterns related to Al responses are unique for each rice species. Our findings provide novel insights into DNA methylation patterns in wild and cultivated rice genotypes and their possible role in the regulation of plant stress responses.

aluminum

DMR‐associated gene

Differentially methylated region

genebody methylated

methylated cytosine

Pearson correlation coefficient

transposable element

TE‐like methylated

unmethylated

INTRODUCTION

Rice is an important crop as it represents the primary food source for over half of the world population. Similarly, it has been established as a biological model for monocots molecular and evolutionary research. Currently, numerous rice studies have shown differential methylation patterns that are related to specific developmental stages (Liu et al., 2017; Stroud et al., 2013; Xing et al., 2015), and with contrasting responses to stress conditions (Feng et al., 2016; Hu et al., 2015; Ou et al., 2012; Schmidt et al., 2018). However, most of these analyses have been restricted to elite Oryza sativa (L.) (cultivated rice) crops, analyzing few stress conditions, mainly salinity and drought (Ferreira et al., 2019; Garg et al., 2015; Wang et al., 2011; Zheng et al., 2013, 2017). Additionally, it has been shown that DNA methylation differences might explain gene expression variations between cultivated and wild species (Li et al., 2012). Nevertheless, the epigenomic divergence between wild and cultivated rice and how the model crop differs from its wild relative has been poorly studied. Among staple food crops, rice is the most aluminum (Al) tolerant cereal (Famoso et al., 2010) and there is genetic variation for Al‐tolerance at the genus/species level, with several genotypes exhibiting highly contrasting tolerance responses (Famoso et al., 2010). Remarkably, greater genetic variability and higher Al‐tolerance levels have been reported for wild rice species (Cao et al., 2011; Mao, 2003). The aluminum (Al) toxicity is one of the most limiting factors for plant growth and crop yield on acid soils (pH < 5), being particularly severe in the tropical and subtropical regions that represent most of the worldwide arable land affected by soil acidity. At low pH, Al turns into the toxic trivalent ion (Al+3) (Kochian et al., 2015), which dramatically affects root elongation and consequently, alters water and nutrients uptake (Rascio & Navari‐Izzo, 2011; Zheng, 2010). The identification of mechanisms that are associated with Al responses have been studied from a genomic and transcriptomic perspective. However, to date, there are no studies on DNA methylation influence over Al‐tolerance in rice, although there is some evidence pointing towards epigenetics as a regulatory factor in Al stress responses in other plant species (Gallo‐Franco et al., 2020).

DNA methylation is the most studied epigenetic modification, and it occurs in plants in three different contexts: CG, CHG, and CHH (where H can be A, C, or T) (He et al., 2011). DNA methylation is a key marker in the regulation of plant responses to environmental conditions, exerting its effects over histone conformation, stress‐responsive gene expression, and altering the regulatory interactions between transcription factors and specific cis‐regulatory elements in the promoter regions (Boyko et al., 2010; Ou et al., 2012; Ueda & Seki, 2020). Similarly, it is well known that DNA methylation can inhibit transposons activity or movement for genome protection (Ou et al., 2012; Zhang et al., 2018). In fact, DNA methylation patterns can generate inherited phenotypic variations, whose differential patterns could be playing a significant role in plant evolution and adaptation (Richards, 2006, 2011).

To address the lack of evidence regarding epigenomic variation between cultivated and wild rice species, we generated single‐base resolution methylome maps for four genotypes of the cultivated rice species O. sativa, and two genotypes of the wild specie O. glumaepatula (Steud) grown under control conditions. Our results showed that the genome‐wide methylation profiles are conserved between both species, but there are several differentially methylated regions (DMRs) with species‐specific methylation patterns. Likewise, we found that overall methylome profiles are not shaped by the response to Al toxicity. However, there exist several DMRs between Al‐tolerant and susceptible species, some of which correlate with differentially expressed genes under Al‐stress conditions previously reported in the literature (Arbelaez et al., 2017; Arenhart et al., 2014). Our findings shed light on the methylation patterns associated with cultivated and wild rice species that will serve as a reference scaffold for future studies of rice genetics and epigenetics.

RESULTS

Genome‐wide DNA methylation patterns in cultivated and wild rice genotypes

In this study, we generated and characterized the methylome of four cultivated and two wild rice genotypes with contrasting responses to Al stress, grown under control conditions. In this way, we were able to analyze the methylation profiles associated with each of these species, and their role as a mechanism associated with Al tolerance.

The whole‐genome bisulfite sequencing (WGBS) was performed in duplicate for four O. sativa genotypes: Azucena and Nipponbare (ssp. Japonica) which are considered as Al‐tolerant genotypes; IR64 (ssp. Indica), reported with intermediate Al‐tolerance, and BGI (ssp. Indica), which is considered a susceptible one. For O. glumaepatula WGBS was performed in duplicate for two genotypes: Og97, an Al‐tolerant genotype, and Og131 a susceptible one. In total, 47–60 million and 55–63 million high‐quality reads were obtained for O. sativa and O. glumaepatula, respectively, whose mapping statistics are presented in Table S1. The genome coverage for all samples ranged between 33‐44X, with an average depth per base between 15‐19X. The cleaned reads for all the samples were aligned to the reference genome: Os‐Nipponbare‐Reference‐IRGSP‐1.0 from the Rice Annotation Project Database (RAP‐DB). Only those sequences that mapped uniquely were considered, therefore, duplicated sequences were removed from all samples. The percentage of unique aligned reads ranged between 49% and 65% for O. sativa and 42% and 43% for O. glumaepatula. The bisulfite conversion rate for all the libraries was above 99.5%. The depth and quality of the sequencing were enough to ensure a high‐quality genome‐wide methylation analysis in all the samples.

The percentage of methylated cytosines (mCs) varied from 10% to 14% in O. sativa and from 8% to 15% in O. glumaepatula samples. The CG context had the highest methylation level for both rice species, that is, number of reads showing mCs for all reads covering the same cytosine site, followed by CHG and CHH contexts. Because DNA methylation has a known role in transposable elements (TEs) silencing and gene regulation, we examined the methylation levels for genes and TEs. This analysis was performed inside genes and TEs bodies, and their 2 kb upstream and downstream regions (Figure 1a). The general methylation patterns for these genomic features were similar for both species, but the methylation levels in the Og131 genotype of O. glumaepatula tended to be lower with respect to the other samples. Consistent with other plant species, for genes and TEs, the methylation level was higher in the CG context, followed by CHG and CHH contexts (Niederhuth et al., 2016). Furthermore, TEs had a higher methylation level compared with genes for all sequence contexts. We also calculated the methylation levels variation for all the genes and TEs (Figure 1b), as well as the number of genebody methylated genes (gbM), TE‐like methylated genes (teM), and unmethylated genes (uM) for all the genotypes (Figure 1c). Overall, no clear differences were observed between species for methylation patterns inside these genomics features.

FIGURE 1

DNA methylation levels in and O. glumaepatula varieties. (a) Average methylation levels for genes and TEs bodies, upstream (−2 kb) and downstream region (+2 kb). Each region was divided into 20 bins and the average methylation level was calculated for each bin. (b) Boxplot showing methylation levels variation located in genes and TEs, for all the sequence context in rice genotypes. Blue bars: ; Orange bars: O. glumaepatula. (c) Number of unmethylated genes (μM), genebody methylated genes (gbM), and TE‐like methylated genes (teM) in rice genotypes.

DNA methylation profiles of TEs in rice

The rice genome was divided into windows of 300 kb and for each one of them we computed (i) the number of genes, Gypsy TEs and Mite TEs; and (ii) the average methylation level per sequence context for each species. To understand the relationship between the DNA methylation patterns and the distribution of genes and TEs throughout the genome, we compare the correlation among the computed characteristics per window (Figure 2a). We used the Gypsy (Class I) and Mite TEs (Class II) because they are the most common and diverse TEs families in rice (Song & Cao, 2017). This approach showed that DNA methylation for CG and CHG context for both species correlated positively with Gypsy TEs density (PCC: .69–.72, Figure 2a), which means that windows with a higher methylation level had a greater number of Gypsy TEs. In contrast, the methylation level correlated negatively with the number of genes and Mite TEs per window. The calculated PCC between methylation level and genes density ranged from −.76 to −.85, whereas for methylation level and Mite density, it varied from −.57 to −.64 (Figure 2a,b). It is worth noting that the CHH context showed a different trend, mainly for O. sativa species. Thus, CHH‐methylation showed a positive correlation with respect to genes and Mite densities, and a negative correlation with respect to the methylation in CG and CHG context, as well as Gypsy density (Figure 2a,b, Figure S1). Even though the rice methylation levels along the genome are differentially related to the Gypsy and Mite TEs density, there is not a clear association between the TEs methylation levels and their impact on the gene structure or regulation. In consequence, we explore TE methylation levels with respect to genes distance, classifying them as close (<2 kb from the gene) and distant (>2 kb from the gene) elements. Our results showed no significant correlations between TEs methylation levels and the distance to the closest gene, but the methylation patterns were clearly different among CG, CHG, and CHH contexts, especially for Gypsy TEs as shown in Figure 3a. In fact, for the CG context, the higher distances to a gene correlate with higher methylation levels, but it is an opposite trend in the CHH context, where higher distances to a gene correspond to lower methylation levels. Likewise, for Gypsy TEs there is a clear difference regarding methylation levels for close and distant TEs in all methylation contexts. For Mite TEs similar methylation levels were observed between these two categories (Figure 3b). In addition, we explored the methylation levels distribution of TEs located in the centromeres or chromosome arms, but we do not find clear differences between them (Supplementary Figure S2).

FIGURE 2

(a) Pearson correlation coefficients (PCC) calculated comparing the methylation levels for and O. glumaepatula, associated with the number of genes, Mite TEs, and Gypsy TEs along the genome. (b) Circos plot showing in a genome‐wide context, genes, Mite TEs and Gypsy TEs distribution in association with average methylation levels in the CG, CHG and CHH context for and O. glumaepatula.

FIGURE 3

(a) Gypsy and Mite TEs methylation levels in relation to the distance to the nearest gene for all the sequence contexts. (b) Boxplot showing methylation level variation between Gypsy and Mite TEs close and distant from genes for all the sequence contexts.

Differential methylation patterns between cultivated and wild rice genotypes

Although methylation patterns were similar between rice species, we carried out a detailed comparison between their methylomes by performing a pairwise comparison between species per each sequence context. Here, the Pearson Correlation Coefficient (PCC) was used to contrast cytosine methylation levels. For this analysis, only the cytosines covered in all samples were considered. The hierarchical clustering of all samples using 1‐PCC distance, showed clustering by species but not grouping by subspecies inside O. sativa for all sequence contexts (Figure S3). In fact, the most Al‐tolerant (Azucena) and the most Al‐susceptible (BGI) genotypes were clustered together. Additionally, we repeated the comparisons considering only the cytosines inside genes and TEs for all the sequence contexts (Figure S4, Figure 4a). Our results showed the same clustering profile as the one obtained using whole genome methylation data.

FIGURE 4

Differential methylation between wild and cultivated rice species. (a) Clustering analysis of rice genotypes according to the methylation levels of mC inside genes and TEs for CG context. (b) Number and location of hyper and hypomethylated DMRs between and O. glumaepatula varieties for each sequence context. (c) Number of DMR‐associated genes found between and O. glumaepatula species. Genes were considered as DMR‐associated genes if DMRs were located at gene‐body (GB), upstream (US) (−2 kb), or downstream (DS) (+2 kb) regions.

Because the rice species are separated according to the methylome correlation analysis (Figure 4a), we decided to explore which regions differ between the analyzed species. Therefore, differentially methylated regions (DMRs) were calculated by a tiling window approach (200 bp window size) between O. sativa and O. glumaepatula. O. sativa genotypes were used as reference for all the pairwise comparisons. In total, we found 1601 DMRs of which 60% were hypermethylated and 40% were hypomethylated in O. glumaepatula with respect to O. sativa (Figure 4b). Notably, hypermethylation was mainly found for the CHH context with most DMRs located in TEs, whereas hypomethylation was mainly for the CG context with most DMRs located in the genebody region. Additionally, for the identification of DMR‐associated genes (DAGs), we tagged the genes overlapping with at least 1 bp of distance to a DMR on its annotated genebody, upstream (−2 kb) or downstream regions (+2 kb). We found a total of 1324 DAGs, of which, 799 genes were hypermethylated and 547 were hypomethylated (Figure 4c).

Differential methylation patterns associated with Al‐tolerance: A comparison between wild and cultivated rice

It is well known that the methylome could be related to stress tolerance phenotypes in plants, so we decided to evaluate whether there are specific differentially methylated regions in relation to established Al tolerance in different rice genotypes. We calculated DMRs between Al‐tolerant and susceptible genotypes within each rice species, using the susceptible genotypes as control for all the pairwise comparisons. For O. sativa we selected the overlapping DMRs between NIP‐BGI and AZU‐BGI pairwise comparisons for further analysis. Our results showed 4633 DMRs for O. sativa of which 38% were hypermethylated and 62% were hypomethylated, with most hypomethylated DMRs in the CHH context (mainly inside TEs) (Figure 5a). Complementarily, 8048 DMRs were found for O. glumaepatula of which 72% were hypermethylated and 38% were hypomethylated (Figure 5b). It is worth noting that most DMRs were hypermethylated at the CHH context contrary to the DMRs found for O. sativa. Lastly, we identified 3024 DAGs for O. sativa and 5484 DAGs for O. glumaepatula (Figure 5c). These DAGs are considered as potentially regulated by DNA methylation patterns between Al‐tolerant and susceptible genotypes within each rice species. Besides, we compared obtained DMRs between species to evaluate whether there are similarities in DNA methylation patterns associated with Al‐tolerance. Here, we found 91 shared DMRs between these two species, 63% hypermethylated and 37% hypomethylated of which 58 were associated with 83 genes.

FIGURE 5

DNA methylation patterns associated to aluminum tolerance in wild and cultivated rice species. (a) Hypo and hyper DMRs detected between Al‐ tolerant and susceptible genotypes for . (b) Hypo and hyper DMRs detected between Al‐ tolerant and susceptible genotypes for O. glumaepatula. (c) Number of DMR‐associated genes in and O. glumaepatula genotypes. (d) Clustering analysis of rice genotypes built by using the average methylation levels of TEs close to Al‐responsive genes for CG context. (e) Clustering analysis of rice genotypes built by using the average methylation levels of TEs close to Al‐responsive genes for CHG context. (f) Clustering analysis of rice genotypes built by using the average methylation levels of TEs close to Al‐responsive genes for CHH context.

To gain more insights into the methylation profiles associated with Al tolerance in rice, we analyzed the methylation patterns of 250 genes that were previously reported as differentially expressed under Al‐exposure (upregulated genes Log2FC ≥ 1 and downregulated genes Log2FC ≤ −1), reported by Arbelaez et al. (2017) and Arenhart et al. (2014). Here we reported 21 DAGs for O. sativa (Table S2) and 37 DAGs for O. glumaepatula (Table S3) that were associated with Al‐responsive genes of which, 11 were shared between both species (Table 1). Additionally, for each genotype, the average methylation level for the upstream, genebody, and downstream regions of these 250 genes were calculated (Figure S5). Then, we compared the average methylation level between genotypes. Our results showed a similar clustering among genotypes for all sequence contexts as previously obtained using whole‐genome methylation data (Figure 1d), where O. sativa and O. glumaepatula belong to different clusters, but inside O. sativa the most Al‐tolerant (Azucena) and the most Al‐susceptible genotype (BGI) are grouped together (Figure S5). These results suggest that there are no differential methylation patterns for all these Al‐responsive genes between Al‐tolerant and susceptible species, grown under control conditions. But a few genes are being epigenetically regulated in relation to the Al‐stress response.

TABLE 1

DMR‐associated genes found in both and O. glumaepatula between aluminum‐susceptible and tolerant genotypes

Gene ID	O. sativa			O. glumaepatula
	Location	Status	Context	Location	Status	Context
Os01g0949900	Upstream	Hypo	CHH	Upstream	Hyper	CHH
Os01g0639600	Upstream	Hypo	CHH	Upstream	Hyper	CHH
Os12g0210500	Upstream	Hypo	CHH	Genebody	Hypo	CHH
	Genebody	Hypo	CHH
Os02g0186800	Upstream	Hypo	CHH	Upstream	Hyper	CHH
Os05g0472400	Upstream	Hypo	CHH	Genebody	Hypo	CHH
	Genebody	Hypo	CHH	Downstream	Hyper	CHH
Os08g0158200	Genebody	Hyper	CHG	Genebody	Hyper	CHG
Os07g0509800	Genebody	Hyper	CHH	Genebody	Hyper	CHH
Os01g0609300	Genebody	Hypo	CG	Genebody	Hypo	CG
				Upstream	Hyper	CHH
Os11g0134900	Genebody	Hypo	CHG	Genebody	Hypo	CHG
Os01g0597800	Downstream	Hyper	CG	Upstream	Hyper	CHH
Os10g0459300	Downstream	Hypo	CHH	Downstream	Hyper	CHH

Note: These DMR‐associated genes were reported as differentially expressed under aluminum stress conditions by Arbelaez et al., 2017 and Arenhart et al., 2014.

Finally, we analyzed the methylation status of 127 TEs that we found close to Al‐responsive genes. Seventy‐nine of these TEs were associated with a DMR in O. glumaepatula and 70 were associated with a DMR in O. sativa. In addition, implementing a PCC analysis, we compared the methylation levels of TEs in proximity to Al‐responsive genes between varieties. Genotypes were clustered according to Al‐tolerance level, at least, for the CG methylation context (Figure 5d–f). This result is especially interesting within O. sativa species where the two most tolerant genotypes (AZU and NIP) were grouped together whereas IR64 (intermediate genotype) and BGI (susceptible genotype) were clustered apart. These results suggest a possible role of TEs in the epigenetic regulation of Al‐genes.

DISCUSSION

DNA methylation is the most studied epigenetic variation having an important role in the regulation of many genomic features such as chromatin structure, transposon silencing, gene expression regulation and recombination. From this perspective, DNA methylation patterns are considered an evolutionary result of life histories in plant species (Niederhuth et al., 2016). However, most epigenetic variation studies have focused on model species.

In this work, we compared DNA methylation patterns between O. sativa, a rice cultivated species, and O. glumaepatula, a South American endemic wild rice. Two genetically and ecologically different species (Stein et al., 2018). Our results showed that, despite the genetic and ecological differences in the studied models, their general methylation profiles are very alike. Here we demonstrate that gene methylation showed the typical behavior reported previously for each context with a similar number of gbM, teM, and uM genes between species. Likewise, for all the methylation contexts, TEs showed a strong methylation level compared with genes. This result agrees with the important role of epigenetics in the TEs silencing, preventing their mobilization along the genome (Jones, 2012; Rabinowicz, 2003) and are in accordance with the global genome‐wide methylation patterns reported in several rice species and genotypes (Feng et al., 2016; Garg et al., 2015; Li et al., 2012; Li et al., 2020; Stroud et al., 2013; Zheng et al., 2017). In fact, previous studies have shown that general methylation trends are conserved in phylogenetically distant plant species (Niederhuth et al., 2016; Zemach et al., 2010).

It has been reported that the methylation status of TEs can regulate nearby genes expression, affecting plant responses to environmental conditions. In this regard, a novel finding from this study is that methylation levels throughout O. sativa and O. glumaepatula genomes were positively correlated with Gypsy TEs density. Conversely, these levels were negatively correlated with genes and Mite TEs density. Previous studies have reported the relationship that exists between methylation level, genes and TEs densities (Bhatia et al., 2018; Zhang et al., 2006). For instance, Bhatia et al. (2018) evaluated the methylation levels of different organs of chickpea plants and reported the overlapping of low methylation levels with high gene density regions, whereas high methylation levels positively correlate with high TEs density. These reports agree with the high methylation level found in the TEs body region, possibly associated with TEs silencing through methylation mechanisms (Jones, 2012; Rabinowicz, 2003). However, the mentioned studies reported a generalized trend of high methylation levels for all TEs. In contrast, our results are showing that methylation patterns vary according to the TEs family and the DNA methylation sequence context.

Our findings showed that Gypsy TEs have lower average methylation levels when they are located close to a gene, whereas for Mite TEs this trend is not clear. It has been previously reported that methylation levels of TEs can affect the transcriptional activity of neighboring genes (Hollister & Gaut, 2009; J. Zhang et al., 2015). For instance, Choi and Purugganan (2018) showed that Class I TEs near to genes have reduced methylation levels in comparison with TEs located far from genetic regions. Likewise, they showed that Gypsy TEs located far from genes heavily methylate as a defensive mechanism against TEs transposition. Therefore, the clear methylation level differences found for Gypsy TEs close and distant to genes, in our study, suggests a possible interaction between TEs and surrounding genes regulation. Regarding methylation levels for Mite TEs, there is not a clear difference between mobile elements located close or distant from genes. Class II TEs (DNA TEs) preferentially insert near or inside genes, contributing to allelic diversity and acting as transcriptional regulators (Dubin et al., 2018). In fact, Class II TEs might have strong regulation effects on gene expression associated with methylation spread (Choi & Purugganan, 2018). Given that Mite family is one of the most diverse TEs families in rice, with almost four times more copy numbers in comparison with other sequenced plant genomes (Song & Cao, 2017; Wang et al., 2014; Zhang et al., 2014), it would be interesting to further study, from a functional perspective, the potential relation between Mite TEs distribution and transcriptional control.

Another result derived from TEs methylation patterns was the high methylation levels of Gypsy and Mite TEs close to genes in the CHH context. The genome‐wide methylation levels for CHH context in several plant species, including our study, are lower than CG and CHG methylation (Niederhuth et al., 2016). However, we found higher CHH methylation levels of TEs close to genes. This behavior has been previously related to the presence of mCHH islands, which are defined as short and highly methylated regions, typically found upstream and downstream genes (Martin et al., 2021). In fact, mCHH islands were reported by the first time associated with Mite TEs in rice plants (Zemach et al., 2010). Likewise, Martin et al. (2021) showed a high level of genes associated with CHH islands (78%) in rice plants compared with other Poaceae species. Although the function of this CHH island is not clear, it has been proposed that they modify chromatin states, to prevent the spread of epigenetic modifications into genes or vice versa (Gent et al., 2013).

So far, we have analyzed the general patterns of methylation shared between wild and cultivated rice genotypes. Samples clustering using whole‐genome methylation as well as specific genes and TEs methylation patterns allowed us to group samples to the species level. We found specific genome regions that differ in their methylation status between O. sativa and O. glumaepatula, representing species‐specific epigenetic marks. Then, we contrasted DAGs found in our study with genes, previously reported as differentially methylated when comparing wild and cultivated rice species (Li et al., 2012). We found 11 common genes differentially methylated in both studies: (Os08g0236800, Os03g0170200, Os04g0460900, Os03g0760000, Os06g0608401, Os06g0690600, Os05g0525800, Os06g0499900, Os09g0434600, Os11g0275500, and Os04g0277400). Likewise, we found two DAGs (Os07g0669500 and Os08g0424500) that have been previously reported as key domestication‐related genes in in rice (Chen et al., 2021; Kovach et al., 2007). Os07g0669500 (FRIZZI PANICLE‐FZP) is a gene that increased the number of secondary branch and grains per panicle (Chen et al., 2021) and Os08g0424500 (SCENTED KERNEL – SK2) gene has been related to the grain fragrance (Kovach et al., 2007). Although more studies are needed, our results provide evidence for the role of DNA methylation during the domestication process in rice plants.

Interestingly, the wild rice O. glumaepatula showed a high number of hypermethylated TEs with respect to O. sativa, specifically in the CHH context. Even though epigenetic variations can arise spontaneously in different organisms (Becker et al., 2011; Graaf et al., 2015; Schmitz et al., 2011), genetic background and environmental conditions are considered the most important factors that structure epigenetic patterns (Kawakatsu et al., 2016). Given that O. sativa and O. glumaepatula are phylogenetically distant and have been subjected to their own evolutionary histories the different methylomes must reflect adaptation processes (Stein et al., 2018). Nevertheless, the analyzed O. sativa genotypes were not clustered according to the subspecies division. It is possible that methylome differences at the sub‐species level could be strongly affected by the environmental conditions in which they have evolved, or by stable epialleles segregation, generated by spontaneous DNA methylation (Zhang et al., 2018). The segregating epimutations have been reported, in several cases, as factors that contribute to inheritable variation, independently of DNA sequence changes (Cortijo et al., 2014; Richards, 2006).

Considering DNA methylation level in different rice genotypes, we wondered whether tolerance to aluminotoxic conditions in contrasting rice varieties could be associated with differential methylation patterns. When analyzing contrasting genotypes in both, O. sativa and O. glumaepatula, we found several DMRs between tolerant and susceptible varieties, suggesting a possible role of DNA methylation in Al‐response regulation. For O. sativa we found mainly hypomethylated regions in the CHH context whereas hypermethylated ones were found in the CHH context for O. glumaepatula. Likewise, most of found DMRs between Al‐tolerant and susceptible varieties were unique for each rice species, suggesting that epigenetic regulation mechanisms related to Al‐tolerance are different between O. sativa and O. glumaepatula. These variations might be associated with different adaptation processes in contrasting environments, for example, O. glumaepatula being a wild species, endemic from South America, has been permanently subjected to acid soils guiding an adaptation process to Al toxicity. Likewise, epimutation accumulation in plants is also a source of heritable epigenetic and phenotypic diversity in plants such as the generation of epialleles related to stress tolerance (Graaf et al., 2015).

Our experimental strategy depicted 91 shared DMRs and 83 common DAGs between O. sativa and O. glumaepatula. Additionally, we found DAGs that have been reported as differentially expressed under Al‐stress conditions (Arbelaez et al., 2017; Arenhart et al., 2014) (Tables S1 and S2). Notably, 11 of these genes were differentially methylated between Al‐tolerant and susceptible genotypes for both species (Table 1). These convergent features represent promising genetic elements that might structure a core machinery in the epigenetic regulation of Al‐tolerance. Further work must elucidate their functional role in Al‐tolerance by evaluating the transcriptional regulation of identified DAGs under stress conditions, helping to confirm the causal relationship between methylation and transcriptional repression or activation. Our results suggest epigenetic mechanisms of aluminum tolerance based on a set of limited core genes rather than generalized regulation mechanisms exerting phenotypic effects through thousand genes.

CONCLUSIONS

According to our results the methylation patterns associated with genes and TEs for both rice species are conserved. Interestingly, we found a positive correlation between the methylation level and the density of Gypsy TEs. But a negative correlation was found between the methylation level and the density of Mite TEs and genes. Likewise, we found for Gypsy TEs a clear difference between methylation levels according to their distance to the closest gene, whereas no differences for Mite TEs were observed. However, there exist several genomic regions with species‐specific methylation patterns, reflecting the own evolutionary histories of O. sativa and O. glumaepatula. Complementarily, we found several regions potentially regulated through epigenetics that are related to Al‐tolerance in rice. We reported different DMRs and DAGs in O. sativa and O. glumaepatula independently. Few of these reported DAGs between Al‐tolerant and susceptible genotypes for both rice species were shared between them, suggesting different mechanism of Al‐response between cultivated and wild rice species, nevertheless, these convergent features, represent promising genetic elements that might structure a core machinery in the epigenetic regulation of Al‐tolerance. Our findings represent a first approach in the understanding of the differential methylation patterns between wild and cultivated genotypes and suggest the role of DNA methylation in the regulation of Al tolerance in rice.

MATERIALS AND METHODS

Plant material and genomic DNA extraction

This study presents the evaluation of the global methylation patterns for the O. sativa genotypes: Azucena (AZU), Nipponbare (NIP), IR64 and BGI9311 (BGI), and Oryza glumaepatula: Og131 and Og97 (Table S1). Seeds used in this study were given by the plant physiology laboratory from ICESI University, Cali, Colombia. After sterilization and break of dormancy, seeds were subjected to dark conditions at 30°C for 4 days. Later, seeds were grown in a culture room at 30°C and 12:12 dark/light conditions for 10 days. Seedlings were transferred to a hydroponic medium with a Kimura B solution (pH 7) and Arnon micronutrients. Roots from 3‐week‐old seedlings were collected and stored at −80°C. Total genomic DNA was extracted from frozen root tissue by CTAB 2X protocol with modifications (Maropola et al., 2015). Genomic DNA quality was evaluated on agarose gels and DNA quantity was measured using a Nanodrop spectrophotometer (Thermo Scientific).

Whole‐genome bisulfite sequencing

Bisulfite‐seq (BS‐seq) libraries were made from genomic DNA isolated from O. sativa and O. glumaepatula seedlings roots. DNA from three independent seedlings was pooled as one sample, and two samples were sequenced per genotype. DNA was first fragmented by sonication from 100 to 300 base pairs (bp) in size, followed by end‐blunting, dA addition at the 3′ end, and ligation of adapters. Next, adaptor‐ligated molecules of 200 to 300 bp were isolated by agarose gel electrophoresis and subjected to a treatment of sodium bisulfite conversion using the ZYMO EZ DNA Methylation‐Gold kit (ZYMO Research Corporation, Irvine, California, USA). Finally, the Polymerase chain reaction (PCR) enriched libraries were purified and subjected to high‐throughput sequencing with Illumina HiSeq X ten platform to achieve an approximately 30X sequencing depth for each sample. Bisulfite‐seq was performed in CD Genomics (CD Genomics Inc., Shirley, New York, USA).

Data filtering, read alignment, and general statistics

The FastQC tool was used to perform basic statistics on the quality of the raw reads. Then, sequencing adapters and low‐quality data of the sequencing data were removed by Trimmomatic (version .36). Pre‐processed reads were aligned to the Os‐Nipponbare‐Reference‐IRGSP‐1.0 genome downloaded from de Rice annotation project database (RAP‐DB) (https://rapdb.dna.affrc.go.jp/download/irgsp1.html) using Bismark v..16.3 (Krueger & Andrews, 2011) and the Bowtie2 v.2.2.8 tool (Langmead & Salzberg, 2012) with default parameters. Only the uniquely aligned reads were maintained and all the samples were de‐duplicated using the Bismark deduplication module. The methylation calling data obtained from Bismark was used for further analysis. The methylation level of each cytosine was defined as the proportion of reads displaying mCs among all the reads covering the same cytosine position.

Analysis of DNA methylation landscape between rice species

A comparative analysis among the rice genotypes methylome was made using the methylation level per cytosine position throughout the complete genome. The methylation level was also calculated for each sequence context (CG, CHG and CHH) inside gene and TE‐body as well as their 2 kb upstream and 2 kb downstream regions. Each region was divided into 20 bins of equal size and the average methylation level for each bin was calculated and plotted. In addition, the number of unmethylated (uM), gene body methylated only in CG context (gbM), and TE‐like methylated (teM) genes were determined using the method reported by (Kawakatsu et al., 2016).

Methylation patterns in relation to genes and TEs

The rice genome was divided in windows of 300 kb and for each one of them we compute (i) the number of genes, Gypsy TEs and Mite TEs; and (ii) the average methylation level per sequence context for each species. We compare the correlation among the characteristics computed per window using a PCC. For the two families of TEs, Gypsy and Mite, the distance between the TE and the closest gene was calculated. Based on this information, TEs were classified into two categories: close TEs as those located less than 2 kb to a gene, and distant TEs as those located more than 2 kb away from a gene. Finally, we defined whether TEs were located inside the centromere region or in the chromosome arms. The coordinates used to define the area of influence of the centromere in each chromosome was 2Mb on each side of the window with the highest frequency of CentO (AA) sequences (Lee et al., 2005).

Identification of differentially methylated regions (DMRs) between rice species

A comparison between methylation level of the cytosines were performed using the Pearson Correlations Coefficients (PCC). Thus, the methylome of each species is compared with the others. The hierarchical clustering between all the comparisons was made using 1‐PCC distance metric. Likewise, differentially methylated regions (DMRs) between rice species were identified using the tiling window approach with a window's size of 200 bp and step size of 100 in the software Methylkit (v.1.16.1.). Only the methylated cytosines covered by ≥10 reads and windows with at least five cytosines each were considered for the analysis. A linear regression model was used to determine statistical significance (q‐value) followed by a Sliding Linear Model (SLIM) correction. A window with methylation difference of 75%, 50%, and 15% for CG, CHG, and CHH respectively, and a q‐value ≤.01 as compared with the reference samples was considered as a DMR. Neighboring DMRs with a gap less than 100 bp were merged. Pairwise comparisons were made between all the O. sativa genotypes against the O. glumaepatula genotypes to determine DMRs (O. sativa genotype was considered as the reference group). Finally, overlapping DMRs for all the pairwise comparisons were selected. The location of DMRs in the genome was defined with at least 1 bp overlapping between the DMR and a functional feature (Genebody, upstream −2 kbps, downstream +2 kbps and TEs) (Sun et al., 2019; C. Wang et al., 2018). Genes overlapping with DMRs in the functional or promoter region were defined as DMR‐associated genes (DMGs).

Pairwise comparisons were also made between Al‐tolerant and susceptible genotypes for O. sativa (NIP‐BGI, AZU‐BGI) and O. glumaepatula (Og97–Og131) using the same procedure described above. The susceptible genotype was considered as the reference group. Overlapping DMRs for the two comparisons made inside O. sativa were selected for further analyses.

DMR‐associated genes identified among Al‐tolerant and susceptible genotypes for O. sativa and O. glumaepatula were compared with a set of 250 genes previously associated to Al stress‐response in rice (Arbelaez et al., 2017; Arenhart et al., 2014). According to the reported experimental data, the 250 genes showed significant changes in expression after Al exposure (upregulated genes Log2FC ≥ 1, downregulated genes Log2FC ≤ −1) (Arbelaez et al., 2017; Arenhart et al., 2014). This comparison allowed us to explore the association between expression and methylation patterns of Al responsive genes.

CONFLICT OF INTERESTS

The authors declare that they have no competing interests.

AUTHOR CONTRIBUTIONS

J.J.G.‐F., M.A.Q. and T.G.‐H conceived the original idea and designed the experiments, J.J.G.‐F. and J.C. performed the experiments, J.J.G.‐F., M.A.R. and F.T.‐T. designed and developed the bioinformatic methods, J.J.G.‐F., M.A.Q. and F.T.‐T. analyzed the data, J.J.G.‐F., M.A.Q., T.G.‐H, M.A.R. and F.T.‐T. write, review, and edit the manuscript. All authors have read and agreed to the published version of the manuscript.

Figure S1. Methylated cytosines (mCs) per window (300Kpbs windows) throughout the genome with respect to the number of genes, Mite TEs and Gypsy TEs.

Click here for additional data file.

Figure S2. Boxplot showing the methylation level variation for Gypsy and Mite TEs in the chromosome arms and centromeric regions.

Figure S3. Pearson correlation coefficients (PCC) heatmap calculated using the methylation level for each cytosine position along the genome among all and O. glumaepatula genotypes A.) Heatmap for the CG context, B) for the CHG context and C) for the CHH context.

Figure S4. Clustering analysis of rice genotypes according to the Pearson Correlations Coefficients (PCC) among rice samples using methylation levels average for all the sequence contexts. A) genes and B) TEs.

Figure S5. Clustering analysis of rice genotypes according to the Pearson Correlations Coefficients (PCC) among rice samples using the average methylation levels for all the sequence contexts. A) genebody, B) upstream region (‐2Kb) and C) downstream region (+2Kb). The analysis was performed using the 250 genes reported as differentially expressed under aluminum exposure in rice (Arbelaez et al., 2017; Arenhart et al., 2014).

Click here for additional data file.

Table S1. Rice methylome alignment statistics with respect to the reference genome.

Table S2. DMR‐associated genes found between Al‐tolerant and susceptible genotypes for

Table S3. DMR‐associated genes found between Al‐tolerant and susceptible genotypes for O. glumaepatula

Click here for additional data file.