Methylation status was significantly reduced after AzaD treatment
To assess the impact of AzaD on rice DNA methylation, we treated Kitaake (Oryza sativa ssp. geng/japonica) rice seedlings with AzaD and performed whole-genome bisulfite sequencing (WGBS). The results of WGBS revealed a global reduction of CG, CHG, and CHH methylation (Figure 1A). Among the three contexts, CG methylation experienced the most dramatic decline (57.6%), from 68.4% in Kitaake to 29.0% in AzaD treated plants. This reduction was of a lesser magnitude than that of osmet1-2 with respect to its wild-type (WT) Nipponbare (Oryza sativa ssp. geng/japonica) plants, in which CG methylation was reduced by 77.4%, from 44.7% in the WT to 10.1% in osmet1-2 mutant plants (Hu et al, 2014) (Figure 1A). Compared with CG methylation, AzaD treatment had a lesser effect on CHG and CHH methylation, but these were also significantly reduced (CHG: from 36.2% to 19.3%; CHH: from 3.9% to 2.8%). In addition, the decrease of CHH methylation in AzaD treated plants was lower than that in osmet1-2 (CHH from 5.21% to 2.98%), and conversely, the decrease of CHG methylation in AzaD treated plants was significantly greater than that in osmet1-2 mutants (25.8% to 22.1%) (Figure 1A).
To better understand the effects of AzaD treatment, we measured changes in DNA methylation levels between AzaD treated and control plants (Figure 1B, S1 and S2). The results suggested that after AzaD treatment, the loss of CG and CHG methylation was relatively more prevalent throughout the whole chromosome and showed no distinct preference for heterochromatin or euchromatin (Figure 1B, S1 and S2). On the contrary, CHH methylation levels were slightly elevated in heterochromatin regions and declined in euchromatin regions. Compared to AzaD treatment, osmet1-2 mutants showed less CHG methylation decline, and CHH methylation was reduced throughout the whole chromosome (Figure 1B, S1 and S2). The reduction of methylation levels in heterochromatin regions was milder than in euchromatin regions in osmet1-2, but methylation levels became evenly reduced after AzaD treatment (Figure 1B, S1 and S2). The methylation levels of genes and TEs were also examined, respectively (Figure 1C). CG methylation levels declined more severely in gene bodies and TEs of osmet1-2 mutants than AzaD treated plants. However, in CHG methylation an opposite pattern was observed, in which the extent of methylation reduction was greater in the AzaD treatment than in osmet1-2, both in genes and TEs, with osmet1-2 methylation levels remaining slightly reduced in gene bodies (Figure 1C). CHH methylation levels of both genes and TEs experienced little change in AzaD treated plants, but was significantly reduced in osmet1-2 mutants (Figure 1C). The difference in CHG and CHH methylation level changes between the AzaD treatment and osmet1-2 mutant plants suggested that AzaD treatment and mutation of OsMET1-2 may affect distinct pathways involved in methylation establishment and maintenance. Real-time quantitative PCR (RT-qPCR) results showed that the expression levels of the primary methyltransferases were significantly down-regulated, which is in correspondence with the whole genome demethylation (Figure 1D). CMT3 was up-regulated in osmet1-2 mutants, which explains the more severe loss of CHG methylation in AzaD treated plants in comparison with osmet1-2 mutants (Hu et al, 2014).
DMR identification and comparison with osmet1-2 mutants
With the aim of studying methylation changes more intuitively, DMRs in AzaD treated and osmet1-2 mutant plants were identified with respect to the check (CK, Kitaake without AzaD treatment) and WT (wild-type Nipponbare) plants, respectively, according to specific standards for each methylation context (Figure 2A). DMRs was determined when the differences in methylation level reach 40% for CG, 20% for CHG, and 10% for CHH (Higo et al, 2020, Shi et al, 2021), respectively. 135,264 CG hypo DMRs were identified in AzaD treated plants. There were almost no CG hyper DMRs identified after AzaD treatment. Even more CHG hypo DMRs (150,256) were identified after AzaD treatment, due to the different standards used between CG and CHG DMR identification. Moreover, comparing to CG and CHG, fewer CHH DMRs were identified: 11,338 CHH hypo DMRs and 3,137 CHH hyper DMRs. (Figure 2A).
The number of CG hypo DMRs in AzaD treated plants, for which 135,264 CG hypo DMRs were identified, was lower than that in osmet1-2 mutants, indicating a more moderate impact of AzaD treatment on CG methylation than osmet1-2 mutation. There was a more extensive CHG reduction and more CHG hyper DMRs in the AzaD treated plants than in osmet1-2. Interestingly, a higher proportion of CHH hyper DMRs were found in AzaD treated plants (AzaD treated rice: 21.7%, osmet1-2 mutants: 4.1%) (Figure 2A).
The mutation of OsMET1-2 mainly affected CG methylation, but according to previous studies, CHG and CHH methylation were reduced as well if CG methylation is first depleted in the same region (Ji et al, 2018). A frequency distribution of CHG and CHH methylation levels in CG hypo DMRs was created to determine if AzaD treatment would cause CHG and CHH methylation reductions in CG hypo DMRs (Figure 2B). The results showed that in the CG hypo DMRs, CHG methylation was down-regulated in the AzaD treated plants, consistent with osmet1-2 mutants. After AzaD treatment, the number of CG hypo DMRs with CHG levels greater than 60% was reduced more than in osmet1-2 mutants, indicating a greater reduction of CHG methylation resulting from AzaD treatment. Surprisingly, the CHH level in CG hypo DMRs was not completely down-regulated. The number of CG hypo DMRs with CHH methylation levels of 5-20% was increased, suggesting that in specific CG hypo DMRs, CHH methylation was activated, which was not observed in osmet1-2 mutants (Figure 2B).
The distribution of DMRs in AzaD treated and osmet1-2 plants are quite similar in many aspects except for some key distinctions (Figure S3). Each kind of DMR shared a similar distribution ratio in the 1 kb regions upstream and downstream of genes, with the highest ratio being that of CHH hypo DMRs (CG hypo: 11% in AzaD, 13% in osmet1-2; CHG hyper: 15% in AzaD, 20 % in osmet1-2; CHG hypo: 12% in AzaD, 15% in osmet1-2; CHH hyper: 13% in AzaD, 16% in osmet1-2; CHH hypo: 36% in AzaD, 33% in osmet1-2). Additionally, CG and CHG DMRs distributed in TEs also shared similarities (CG hypo: 18% in AzaD, 17% in osmet1-2; CHG hyper: 9% in AzaD, 7% in osmet1-2; CHG hypo: 23% in AzaD, 25% in osmet1-2). A similar pattern was also observed in CHH hypo DMRs within genes, between AzaD treated and osmet1-2 plants(CHH hypo: 17% in AzaD, 15% in osmet1-2). Moreover, variations in distribution ratios were also found between the AzaD treatment and osmet1-2. More CG and CHG hypo DMRs in osmet1-2 were found to locate within genes (CG hypo: 20% in AzaD, 34% in osmet1-2; CHG hypo: 7% in AzaD, 21% in osmet1-2), while the number of CHH hyper DMRs located in TEs was larger in AzaD treated plants than that in osmet1-2 (CHH hyper: 15% in AzaD, 7% in osmet1-2) (Figure S3). This indicates that AzaD treatment may affect different methylation processes compared to mutation of OsMET1-2.
To further compare the similarities and differences between AzaD treated and osmet1-2 plants, we identified syntenic regions between the Nipponbare and Kitaake rice cultivar genomes and compared the overlap ratio between DMRs by judging whether the DMR corresponding regions in the two backgrounds overlapped (Figure 2C). As a result, CG hypo DMR showed the highest overlap ratio (87.91%) between AzaD treated plants and osmet1-2 mutants (Figure 2C), suggesting high consistency in the area affected by AzaD treatment and osmet1-2 mutation on CG context. In addition, CHG and CHH hypo DMRs showed moderate overlap levels, with overlap ratios of 39.22% and 57.87%, respectively (Figure 2C). However, CHG and CHH hyper DMRs had overlap ratios of only 2.56% and 6.88%, respectively (Figure S4). The high overlap ratio of hypo DMRs indicated a similar demethylation pattern.
In all, though AzaD treatment primarily caused demethylation in the CG context, CHG and CHH methylation levels were significantly decreased as well. Comparative analysis with osmet1-2 mutants also revealed similarities in the CG methylome between the AzaD treated plants and osmet1-2 mutants. These results suggest that AzaD treatment in rice can mimic osmet1-2 mutation to resolve the difficulty of obtaining CG demethylated materials.
TEs were significantly activated according to distance from genes
TEs are an important part of the genome and have the ability to move from one location in the genome to another (Lanciano and Cristofari, 2020). DNA methylation is a primary method used to suppress transposon expression. Although methylation levels were reduced both in the AzaD treated and osmet1-2 plants, overall gene expression levels showed no significant elevation (Figure S5). However, TEs in AzaD treated and osmet1-2 plants both showed significant activation (Figure 3A), which is in accordance with the observation that TE expression is mainly affected by methylation levels (Deniz et al, 2019, Ewing et al, 2020, Wambui Mbichi et al, 2020). Further experiments were carried out to reveal how TEs respond to AzaD treatment.
TEs are often separated into two major classes: class I TEs, also called retrotransposons, utilize an RNA intermediate that is reverse transcribed before genomic reinsertion; class II TEs, or DNA transposons, move via excision from one location in the genome followed by insertion into another (Capy et al, 1997, Takata et al, 2007). An inherent difference in methylation levels was found in CK plants between class I retrotransposons and class II DNA transposons. CG and CHG methylation levels in class I retrotransposons are higher than that in class II DNA transposons (Figure 3B) and both were down-regulated in AzaD treated plants. Contrastingly, CHH methylation levels were higher in class II DNA transposons than in class I retrotransposons. Interestingly, after AzaD treatment, CHH methylation levels were increased in class I retrotransposons, while CHH methylation of class I TEs in osmet1-2 plant were significantly down-regulated (Figure S6). In addition, class II DNA transposons were demethylated in the CHH context, to the same methylation levels as the CG and CHG contexts (Figure 3B). This was accompanied with preferential distribution of class I TEs in heterochromatin, explaining the CHH hyper methylation in heterochromatin (Figure S7). Demethylation may cause TE activation because of the negative association between methylation and TE expression. We characterized the effects of AzaD-induced methylome changes on TE expression. The results showed that 5,118 of 19,788 TEs (22.53%) showed significantly different expression levels after AzaD treatment, of which 4,459 TEs were up-regulated (log2FoldChange(FPKM+1) > 1) and 659 TEs were down-regulated (log2FoldChange(FPKM+1) < -1).
We further compared the response of the two classes of TE to AzaD treatment. Some ‘dead’ TEs were removed as they weren’t expressed both in CK and AzaD treated rice, while the remaining TEs consisted of 5,291 live class I retrotransposons and 2,207 class II DNA transposons. After AzaD treatment, 63.03% of live class I TEs was activated (log2FoldChange(FPKM+1)>2), which was higher than that of class II TEs activated (50.92%) (Figure 3C left). This indicates that class I TE expression is more sensitive to AzaD treatment than that of class II TEs. The different methylation changes between class I and II TEs brings into question the essential reason for this phenomenon and whether more undetermined differences between class I retrotransposons and class II DNA transposons exist.
TE expression patterns were found to be relevant to their distance to the closest gene, with class I TEs being distributed more in heterochromatin regions (Rebollo et al, 2011, Eichten et al, 2012), so we hypothesized that the difference in activating ratios of the two TE classes was related to their distance from genes. To confirm this hypothesis, TEs were classified into two groups: close (distance to closest gene less than 2 kb) and far (distance to closest gene more than 2 kb), and higher activating ratios were found in the far TEs (Figure 3C right). More precisely, TEs were classified into further groups according to their superfamilies and sorted by their average distance to the closest gene (Figure S8, 3D). In each group, activated TEs had a farther distance from the closest gene (Figure 3D), the discovery of which further confirmed our previous hypothesis that the activating ratio of the two TE classes was due to their distance from genes.
After classification into different families, TEs from different families showed divergent methylation patterns. Short interspersed nuclear elements (SINEs), belonging to class I retrotransposons, showed the highest CHH methylation levels among all class I retrotransposons (Zhang et al, 2015), which was also illustrated in our data (Figure S9). We speculated the high CHH methylation levels in SINEs may also be associated with their distance to the nearest genes. To explain the phenomenon of SINEs having abnormally higher CHH methylation levels, the methylation level of each TE family was calculated to reveal correlations between TE methylation levels and distance to the nearest gene. TE families farther from genes showed higher CG and CHG methylation levels (Figure S9A, B). An opposite pattern was observed for CHH methylation levels, with more distant TEs exhibiting relatively lower CHH methylation levels (Figure S9C). SINEs, as a family from class I retrotransposons, were closer to genes unlike other families of class I TEs (Figure S9C). This abnormal distance to genes led to their higher CHH methylation levels than other families of class I retrotransposons. The correlation between methylation changes and distance to the nearest gene of each TE was also found after AzaD treatment. Overall, CG methylation levels decreased more severely in more distant TEs (Figure S10), which might be the reason for activation of these TEs. Conversely, more distant TEs exhibited CHH hypermethylation after AzaD treatment (Figure 3E). A similar pattern in CG methylation also happened in osmet1-2mutants (Figure S10), while the CHH hyper methylation was not found in osmet1-2 mutants (Figure S11).
24-nt siRNAs involved in the RdDM pathway were altered as a rapid response to CG methylation loss
RdDM pathway, which mainly involves 24-nt siRNAs, is a vital pathway in establishing de novo methylation, especially CHH methylation (Matzke and Mosher, 2014, Cuerda-Gil and Slotkin, 2016, Wendte and Pikaard, 2017). To explore the mechanism of CHH methylation changes, we performed small RNA sequencing (sRNA-seq) on rice leaf tissue of CK and AzaD treated rice plants. The abundance of 21-22-nt siRNAs were found to be slightly increased while that of 24-nt siRNAs was decreased (Figure 4A). The overall changes in 24-nt siRNA abundance was in correspondence with CHH demethylation, with a correlation coefficient of 0.389 (Figure 4B, S12). Next, we calculated regional abundance of 24-nt siRNAs in DMRs of each sequence context (CG, CHG, and CHH) in CK and AzaD treated rice plants. The highest RPKM of siRNAs was found in CHH hypo DMRs before treatment and the abundance was significantly down-regulated after AzaD treatment. A significant increase was also found in the abundance of 24-nt siRNA in CHH hyper DMRs, while no obvious differences in 24-nt siRNA abundance were found in other DMRs (Figure 4C). These results demonstrate that changes in CHH were regulated by 24-nt siRNAs. Considering the correlation between 24-nt siRNA abundance and CHH methylation, we constructed a boxplot of 24-nt siRNA RPKM for far and close TEs of each family (Figure 4D), and the distribution was highly consistent with that of CHH methylation levels described above (Figure 3E). These finding further confirmed the close correlation between CHH methylation and 24-nt siRNA abundance in TE regions.
Next, we compared CG methylation changes in CHH hyper/hypo DMRs and random regions, and the results showed that CG methylation levels of CHH hyper DMRs were more significantly decreased than CHH hypo DMRs and random regions (Figure 4E). This illustrated that CHH hyper methylation was mainly found distributed across regions with severe CG methylation loss. This result also explains that the rice plants likely induced de novo CHH methylation with 24-nt siRNAs through the RdDM pathway to compensate for the severe loss of CG methylation and to stabilize the whole genome after AzaD treatment.
AzaD treatment causes severe phenotype changes in Kitaake rice plants
Rice seedling and root growth was severely retarded after AzaD treatment, and the stem and leaves became etiolated (Figure 5A). To characterize the effects of AzaD-induced methylome changes on gene expression, we performed RNA-sequencing (RNA-seq) on leaf tissues of CK and AzaD treated rice plants. Differentially expressed genes (DEGs) were identified in AzaD treated plants (Figure 5B). In all, 1,876 genes showed significantly up-regulated expression levels and 1,730 genes were down-regulated. We further compared the overlap rates of DEGs in AzaD treated plants and osmet1-2 mutants. In total, 1,876 genes were up-regulated in the AzaD treatment and 422 genes were up-regulated in osmet1-2 mutants; only a few genes (162) showed opposite regulation patterns and 1,089 genes showed no obvious change in osmet1-2 mutants (Figure 5B). In all, 1,730 genes were down-regulated after AzaD treatment, of which 322 genes were also down-regulated in osmet1-2 mutants; 115 genes were up-regulated in osmet1-2 mutants and the remaining 1,195 genes showed no change (Figure 5B). GO analysis showed that genes up-regulated in both backgrounds were involved in ‘secondary metabolic biosynthetic process’, ‘secondary metabolite process’, ‘organic acid biosynthetic process’, ‘carboxylic acid biosynthetic process’, among others (Figure S13A). Genes down-regulated in both backgrounds were involved in ‘regulation of hormone levels’, ‘hormone metabolic process’, and ‘meta ion homeostasis’, among others (Figure S13B). Genes only up-regulated in AzaD treated plants were involved in ‘response to water deprivation’, ‘response to water’, and ‘response to heat’, among others, suggesting that AzaD treatment activated stress response genes (Figure S14A). Genes only up-regulated in osmet1-2 mutants were involved in ‘cell cycle process’, ‘chromosome’, ‘cytoskeleton’, and ‘DNA replication’, among others, which are all essential processes and components for survival (Figure S15A). However, genes involved in biological processes such as ‘response to drug’, ‘response to water’, ‘response to lipid’, and ‘response to cold’ were down-regulated in osmet1-2 mutants (Figure S15B). Genes involved in ‘thylakoid’ and ‘envelope’ were down-regulated in AzaD treated plants (Figure S14B). Genes related to abiotic stress were mainly up-regulated after AzaD treatment but were down-regulated in osmet1-2 mutants. These different gene regulation patterns suggest there are different impacts of AzaD treatment and mutation of OsMET1-2 on the whole genome.
According to RNA-seq analysis, yellow-green leaf 1 (YGL1) and YGL3 were found to be down-regulated after treatment, as was the gibberellin-enhanced gene OsGAE1 (Figure 5C). Interestingly, the repressed expressions of these genes were not caused by methylation changes, since methylation was down-regulated rather than up-regulated, which should have caused up-regulated gene expression (Figure S16). These down-regulated genes shared a common characteristic in that a transposon was demethylated in the body flanking the gene, indicating that the transposons might disturb the expression of their flanking gene and cause corresponding phenotypes.