Non-CG methylation is dynamically remodeled during rice male gametogenesis

To study DNA methylation dynamics during male gametogenesis in rice, we manually isolated male meiocytes, unicellular microspores and sperms of the Zhonghua11 (ZH11) variety as previously described [17, 20,21,22]. About 400 meiocytes, 300 unicellular microspores and 100 sperm cells were collected for bisulfite sequencing (BS-seq) analysis using a protocol developed for small numbers of cells (Additional file 1: Fig. S1a, b) [23]. Data with two biological replicates were obtained (Additional file 1: Fig. S1c, Additional file 2: Table S1). Violin plots of the BS-seq data revealed that the overall CG methylation (mCG) levels (especially in TE and TE-related genes, TEG) gradually increased during male gametogenesis (Fig. 1a), whereas CHG methylation (mCHG) levels were first increased in unicellular microspore (UM) but subsequently decreased in sperm (S) to the lowest levels (Fig. 1a). Density plots confirmed the mCHG variations between microspore (UM) and meiocyte (Me) and between sperm (S) and microspores (UM) (Fig. 1b). A similar trend of CHH methylation (mCHH) variation was also observed (Fig. 1b). Scanning of differentially methylated regions (DMRs, defined within 100-bp windows, see Methods) between microspore and meiocyte (UM-Me) detected more hyper than hypo CG (876 hyper versus 646 hypo), CHG (10,720 hyper versus 3,764 hypo), and CHH (31,501 hyper versus 20,983 hypo) DMRs, indicating a clear gain of non-CG methylation in microspore (Fig. 1c, d). In sperm relative to microspore, there were more hypo- than hyper-DMRs at non-CG, especially CHG context (30,455 hypo- compared to 1,727 hyper-DMRs), confirming a clear decrease of mCHG in sperm. The mCHG levels in sperm were the lowest when compared with those in egg and zygote or somatic tissues of the same rice variety (Additional file 1: Fig. S2a, see below), which could be also observed in other rice varieties (Additional file 1: Fig. S2b) [4, 16, 17, 24]. Analysis of sex cell methylomes obtained from the Nipponbare (NIP) variety indicated that sperm mCHG and mCHH levels were lower than pollen vegetative cell levels, but comparable to central cells of the female gametophyte (Additional file 1: Fig. S2b) [16, 24]. About 2/3 (5,221/7,775) of the increased mCHG (hyper DMRs) in microspore (relative to meiocyte) were erased in sperm, the other 1/3 (2,554/7,775) microspore hyper DMRs remained hyper-methylated in sperm (Additional file 1: Fig. S3a, clusters A and B). These two clusters of the microspore hyper DMRs were located mainly in genic and intergenic regions (Additional file 1: Fig. S3b), and showed enrichment of euchromatin histone marks (H3K36me3, H3K9ac) (Additional file 1: Fig. S3d). By contrast, many sperm hypo DMRs (12,588) that showed no change between microspore and meiocyte (Additional file 1: Fig. S3a), which corresponded mainly to TE and TEG (Additional file 1: Fig. S3b), and were enriched for the heterochromatin mark H3K9me2 (Additional file 1: Fig. S3d). The data indicate that during male gametogenesis mCHG and mCHH at both euchromatin and heterochromatin loci are dynamically remodeled to the lowest levels in sperm.

Fig. 1

DNA cytosine methylation in rice meiocytes, unicellular microspores, sperms. a Violin plots showing overall cytosine methylation levels (mCG, mCHG, and mCHH) in transposable elements (TE), transposable gene (TEG) and protein coding gene (Gene) of rice Zhonghua 11 (ZH11) meiocyte (Me), unicellular microspore (UM) and sperm (S). Values of the methylomes are averages from the two replicates. The average methylation levels (white dots) and median values (black bars) are indicated. b Density plot showing the frequency distribution of fractional methylation difference between the indicated samples. c Numbers of differentially methylated regions (DMRs) of between the indicated comparisons, distributed in TE (> 500 bp), TEG, gene, and intergenic regions. DMRs located in TE (red), gene (light green), intergenic region (pink), and TEG (yellow) are shown. d Genome browser screenshots of mCG, mCHG, and mCHH in meiocytes (Me), unicellular microspore (UM), sperm (S). Differentially methylated regions are grey colored

From cluster A of the microspore hyper DMRs, 315 genes were identified. These genes are enriched for translation, ribonucleoprotein complex, and cellular protein metabolic functions, implying that mCHG in protein translation and RNA-binding pathway genes was particularly dynamic during meiocyte to sperm development (Additional file 1: Fig. S3e, Additional file 4: Table S3). Several genes showed lower expression in microspore than meiocyte and/or sperm (Additional file 1: Fig. S4a), suggesting that hypermethylation might play a role in their repression in microspore.

CMT3a and CMT3b function during male gametogenesis

CMT3a is a major CHG methyltransferase gene expressed at high levels throughout the sporophytic development in rice [25, 26]. However, its expression became low or undetectable in sperm in several rice varieties (Fig. 2a). By contrast, CMT3b expression was low or undetectable in vegetative tissues/organs but showed expression in reproductive cells including meiocyte, microspore, and sperm (Fig. 2a). To study the function of CMT3 genes during male gametogenesis, we produced cmt3a and cmt3b knockout (KO) plants in the ZH11 background using the CRISPR technique and two independent lines for each gene were obtained (Additional file 1: Fig. S5a). The cmt3a mutants produced mainly defective pollens and were infertile (Fig. 2b). Cytology sections revealed that the cmt3a pollen development was arrested likely at the bicellular microspore stage (Additional file 1: Fig. S5b).

Fig. 2

Effects of cmt3a and cmt3b mutations on DNA methylation in meiocyte, microspore and sperm. a Transcript levels in FPKM of rice CMT3a and CMT3b in seedling (Se), roots (Ro), meiocyte (Me), unicellular microspore (UM), sperm (S), egg (E), zygote (Z), endosperm nuclei (En, 1.5 days after fertilization) and globular embryo (GE, 3 days after fertilization) from RNA-seq data. The sperm (Kit-S) in Kitaake background was reported by Anderson et al., (2013). b The pollen grains of wild type and cmt3a and cmt3b mutants were I2-KI stained. Bars = 50 μm. c Violin plots comparing overall cytosine methylation levels of wild type and cmt3a and cmt3b mutant meiocyte (Me), unicellular microspore (UM) and sperm (S). The average methylation levels (white dots) and median values (black bars) in transposable elements (TE) are shown. Values of the methylomes are averages from the two replicates. d Number of differential methylated regions (DMR) in cmt3a and cmt3b relative to wild type. Relative portions in TE (> 500 bp), TEG, gene, and Intergenic regions are indicated by different colors. e Venn diagrams showing overlapping of hypo-CHG DMRs in cmt3a and cmt3b meiocyte (left) and sperm (right) relative to wild type cells. f Box plots of DNA methylation levels of hypo-CHG DMRs in meiocyte (Me) versus microspore (UM) (upper) and sperm (S) relative to microspore (UM) (lower) in wild type, cmt3a (3a) and cmt3b (3b) cells. The significance was calculated with multiple comparison tests. Different letters on top of the bars indicate a significant difference (p < 0.05). g Genome Browser screen captures showing high CHG methylation sites in microspore relative to meiocyte and sperm decreased in cmt3b mutants (highlighted by grey)

By contrast, the cmt3b mutants showed no clear plant morphological phenotype. To check the effects of cmt3 mutations on DNA methylation during male gametogenesis, we performed BS-seq of meiocyte, unicellular microspore, and sperm isolated from two independent CRISPR/Cas9-free lines of cmt3a and/or cmt3b at T3-4 generation (as well as a tissue culture-regenerated wild type line) (Additional file 2: Table S1). Violin plots of the data revealed that mCHG was almost absent from cmt3a meiocyte (Fig. 2c, Additional file 1: Fig. S6a), consistent with the CMT3a loss-of-function effects in somatic tissues [25]. To a lesser extent, mCG was also reduced in cmt3a meiocyte. However, cmt3a sperm mCHG (as well as mCHH) levels became higher than the mutant meiocyte (Fig. 2c), suggesting that additional activities partially restored mCHG and/or ectopically mediated mCHH in the mutant sperm. DRM2 and CMT2 being highly expressed in rice sperm (Additional file 1: Fig. S7a), the residual mCHG level in cmt3a sperm could be maintained by RdDM or CMT2. This hypothesis is supported by the observation that the drm2/cmt2 mutations also reduced the mCHG levels of those loci in leaves (Additional file 1: Fig. S7b).

By contrast, the cmt3b mutation led to a clear loss of overall mCHG in microspore and sperm but the mutation effect was less clear in meiocyte (Fig. 2c, Additional file 1: Fig. S6a). The cmt3b mutation also resulted in some increases of mCHH in sperm. Density plots confirmed the observations (Additional file 1: Fig. S6b). The increases of mCHH in cmt3a/b sperm might be of indirect effects to compensate mCHG loss in the mutants. Nearly all of the cmt3b hypo-CHG DMRs in meiocyte and sperm overlapped with those of cmt3a (Fig. 2d, e), indicating that CMT3b functioned to maintain mCHG on a fraction of the CMT3a targets. The cmt3b mutation resulted in a large number (33,412) of hypo-CHG DMRs in microspore (Fig. 2d), and diminished the mCHG difference between microspore and sperm observed in wild type (Additional file 1: Fig. S6c). In fact, the methylation levels of hyper-CHG DMRs in wild type microspore versus meiocyte were decreased to the meiocyte levels in cmt3b microspore, and the methylation levels of the hypo-CHG DMRs in wild type sperm versus microspore were decreased to the sperm levels in cmt3b microspore (Fig. 2f, g). For the three clusters of DMRs shown in Additional file 1: Fig. S3a, the cmt3b mutation largely reduced the mCHG levels in microspore (Additional file 1: Fig. S3c). In addition, the 315 genes (from cluster A) showed higher mCHG in exons than introns in microspore (Additional file 1: Fig. S4b), while the cmt3b mutation reduced the mCHG levels from both exons and introns, suggesting that CMT3b may preferentially target gene exons in microspore (Additional file 1: Fig. S4a, b). The analysis indicates that CMT3b is required for the increase of mCHG in microspore.

Histone demethylases JMJ706 and JMJ707 reduce CHG methylation

As CMT3a is not or lowly expressed in sperm, mCHG diluted by mitosis may not be maintained in sperm. Alternatively, active DNA demethylation may be involved, as mutations of DNA demethylases locally remodeled DNA methylation in sperm [17]. Since mCHG is linked to H3K9me2 through a positive feedback loop [27, 28], we investigated whether H3K9me2 demethylases were also involved in the decease of mCHG in sperm. JMJ706 was shown to function as a H3K9 demethylase in rice [29]. JMJ707 is closely related to JMJ706 [29], of which JMJ707 showed high expression in sperm (Additional file 1: Fig. S8a). To test whether the genes were involved in mCHG during male gametogenesis, we produced jmj706 and jmj707 double knockout (KO) plants and obtained two independent lines (Additional file 1: Fig. S8b). The KO lines (j67) showed a reduced pollen viability and seed setting rate (Additional file 1: Fig. S8c, d). We analyzed the DNA methylome of male meiocyte and sperm of the mutant lines (Cas9-free at T3-4 generation) and found that the mutations had no drastic effect on the overall methylation in the cells (Fig. 3a, Additional file 1: Fig. S9a). However, the methylation levels of the hypo-DMRs in wild type meiocyte versus microspore were increased in j67 meiocyte (Fig. 3b, c, d). Similarly, the methylations levels of the hypo- DMRs in wild type sperm versus microspore were augmented in j67 sperm (Fig. 3b, c, d). The jmj706/7 mutations clearly elevated mCHG levels of cluster B DMRs in meiocyte and cluster C in sperm (Additional file 1: Fig. S3c). There was no overlap between cmt3b and jmj706/7-affected DMRs (Additional file 1: Fig. S9b). The analysis indicated that JMJ706/707 play a role to reduce mCHG at a set of genomic loci (mainly TE or TEG) in male meiocyte and sperm.

Fig. 3

Effects of jmj706/707 mutations on DNA methylation in male sex cells. a Comparison of overall TE methylation levels in wild type and jmj706/707 mutant meiocyte (Me) and sperm (S). Values of the methylomes are averages from the two replicates. The average methylation levels (white dots) and median values (black bars) are shown. b Density plots of CHG methylation differences between jmj706/707 mutant and wild-type meiocyte (upper) and sperm (lower) (black lines). The red traces are density plots confined to the hypo DMRs between meiocyte and microspore (Me-UM) (upper) or between sperm and microspore (S-UM) (lower). c Box plots showing DNA methylation levels of 50-bp hypo-CHG methylation regions in wild type meiocyte (Me) (upper) and sperm (S) (lower) relative to microspore (UM) and in jmj706/707 meiocyte (j67-Me) and sperm (j67-S). The significance was calculated with multiple comparison tests. Different letters on top of the bars indicate a significant difference (p < 0.05). d Genome browser screen shots showing low CHG methylation sites in meiocyte (upper) or sperm (lower) relative to microspore but elevated in j67 mutants. Grey illustrates differentially methylated regions

Effect of the cmt3a/b and jmj706/7 mutations on sexual lineage-specific methylation

It is shown that Arabidopsis male sex cells show sexual-lineage-specific methylation (SLM) or sexual-lineage-hypermethylation (SLH) [18]. Using the published methods [18], we identified 555 SLH loci in rice meiocyte, microspore and sperm relative to somatic cells (seedling) (Additional file 1: Fig. S10a). The cmt3b mutations reduced the SLH levels in all three male sex cell types, especially in microspore (Additional file 1: Fig. S10b, c). By contrast, the jmj706/7 had no clear effect on SLH in the sex cells (Additional file 1: Fig. S10b). Further analysis could divide the 555 SLH into 340 canonical SLH and 215 SLM loci (Additional file 1: Fig. S11a) [18]. The canonical SLH loci corresponded mainly to TEs, while the SLM loci were located mainly in genes (body and promoter regions) (Additional file 1: Fig. S11b), suggesting that SLM mainly targets genic regions during male gametophyte and sperm development. In total, 132 genes were targeted by SLM, which are enriched for translational function (Additional file 1: Fig. S11c, Additional file 5: Table S4). The 132 genes appeared to be repressed in sperm compared to meiocyte or microspore (Additional file 1: Fig. S11d, e). The cmt3b mutation reduced SLM and increased expression of some of the genes in sperm (Additional file 1: Fig. S11d, e), suggesting that CMT3b may be involved in SLM and repression of some of the genes in sperm.

Function of CMT3a and CMT3b in egg and zygote DNA methylation

Unlike in sperm, CMT3a is highly expressed in egg and zygote. CMT3b transcripts are detected in Egg and zygote (Fig. 2a). To study CMT3 function in egg and zygote, we compared wild type, cmt3a and/or cmt3b egg and zygote methylomes by BS-seq analysis (Additional file 2: Table S1). In wild type, egg and zygote mCHG levels were higher than sperm (Fig. 4a, Additional file 1: Fig. S12a). Because cmt3a was infertile, we only analyzed cmt3a egg methylome. As in male meiocyte and somatic tissues [25], the cmt3a mutation eliminated almost all mCHG in egg (Fig. 2c and Fig. 4a, Additional file 1: Fig. S12a). The cmt3b mutation had a limited effect on overall mCHG in egg, but caused a clear decrease of mCHG in zygote (Fig. 4a, Additional file 1: Fig. S12a). The cmt3b mutation produced more hypo-CHG DMRs (23,460) in zygote than egg (13,249) or sperm (17,576) (Fig. 4b). About 24% (5,623/23,460) of the hypo-DMRs in cmt3b zygote overlapped with the hyper-DMRs in wild type zygote versus sperm (Fig. 4c). In fact, the cmt3b mutation reduced the methylation differences of the DMRs between wild type zygote and sperm (Z-S) or egg (Z-E) (Fig. 4d, e). Together, the data indicate that CMT3b participates in reestablishing mCHG methylation at a subset of the Z-S and Z-E DMRs in the zygote.

Fig. 4

Effect of cmt3a and cmt3b mutations on DNA methylation in zygote and/or egg cells. a Comparison of overall TE methylation levels in sperm (S), egg (E) and zygote (Z) of wild type and cmt3a, cmt3b mutants. Values of the methylomes are averages from the two replicates. The average methylation levels (white dots) and median values (black bars) are shown. b Number of differential methylated regions (DMR) in the indicated mutant cells relative to wild type. Different colors indicate the distribution of DMR in TE (red), gene (light green), intergenic region (pink), and TEG (yellow). c Venn diagrams showing overlapping of hyper-CHG DMRs in zygote relative to sperm (Z-S) (upper) or egg (Z-E) (lower) and hypo-CHG DMRs in cmt3b zygote relative to wild type zygote (3bZ-Z). d Box plots showing DNA methylation levels of Z-S (upper) or Z-E (lower) hyper-CHG DMR of the indicated cell type. The significance was calculated with multiple comparison tests. Different letters on top of the bars indicate a significant difference (p < 0.05). e Screenshots of CHG methylation levels of 5 representative genes in the indicated cell

Function of rice CMT3b in gene expression in reproductive cells

To study the cmt3 mutation effects on gene expression in sex cells, we first performed RNA-seq of wild type and cmt3b male meiocyte and sperm. Three replicates (two replicates for WT sperm) were performed (Additional file 1: Fig. S13a, Additional file 3: Table S2). PCA analysis indicated a high reproducibility of the replicates (Additional file 1: Fig. S13b). The WT meiocyte transcriptome showed high correlation with previously published rice meiocyte transcriptomes (r > 0.85) (Additional file 1: Fig. S13c). In cmt3b meiocyte, 2,259 and 1,639 genes were respectively up- and downregulated (Fig. 5a). Similar numbers of differentially expressed genes (DEGs) were detected in the mutant sperm (Fig. 5a). Upregulated genes in cmt3b meiocyte were enriched for gene transcription function, while upregulated genes in cmt3b sperm were mainly enriched for translational function (Additional file 1: Fig. S14a, b). A small number of upregulated DEGs were found to associate with hypo-DMRs in the mutant cells (Fig. 5b, c; Additional file 6: Table S5). The analysis suggests that CMT3b plays a role in shutdown of transcriptional activity in meiocyte for preparation of meiosis and may contribute to the low translational activity in sperm [30].

Fig. 5

The cmt3b mutation affected non-TE gene expression in meiocyte and sperm.a Number of differentially expressed non-TE genes (pink) and TE-related genes (TEG) (red) in cmt3b mutant meiocyte and sperm relative to wild type (padj < 0.01, FC > 2). The highly expressed genes (TPM > 10) in sperm were filtrated for comparison. b Number of upregulated genes overlapping with hypo-CHG methylated genes (DMG) detected in cmt3b meiocyte and sperm relative to the wild type cells. c Genome browser screenshots of the methylation and expression levels of representative genes of the 148 (left) and 122 (right) genes shown in (b)

In parallel, we analyzed the egg and zygote transcriptomes of wild type and cmt3a/b plants. Because cmt3a was infertile, we only analyzed cmt3a egg transcriptome. PCA analysis indicated high levels of reproducibility of the replicates. The cmt3b egg and zygote transcriptomes were close to, but distinct from, the wild type (Additional file 1: Fig. S13b). By contrast, the cmt3a egg transcriptome was largely distal from the wild type (Additional file 1: Fig. S13b), consistent with the drastic effect of cmt3a mutation on mCHG in egg. There were in total 4,868 upregulated genes (> two-fold, Q < 0.01), of which 1,648 were hypomethylated at CHG sites in cmt3a egg (Additional file 1: Fig. S12b, c). Interestingly, among the up-regulated genes, 1,461 were TEGs, of which 982 (67.2%) were hypomethylated at CHG sites in cmt3a egg (Additional file 1: Fig. S12b, c). The analysis indicates that CMT3a-mediated mCHG is required mainly for TEGs repression in egg, consistent with previous results showing that cmt3a mutation resulted in burst of TE expression [26]. The cmt3b mutation resulted in totally 3,022 upregulated genes (> two-fold, Q < 0.01) in egg, of which few were TEGs and only 120 (including 26 TEG) were identified as hypo-CHG methylated genes in the mutant egg (Additional file 1: Fig. S12b, c), indicating that unlike cmt3a, the cmt3b mutation de-repressed mainly non-TE-related genes, which was likely independent of a clear loss of mCHG. However, about 40% of the upregulated genes in cmt3b overlapped with those detected in cmt3a eggs (Additional file 1: Fig. S12d), suggesting that both CMT3 genes are required for gene (mainly non-TEGs) repression in egg. The cmt3b mutation resulted in upregulation of 2,179 and downregulation of 1,870 genes in zygote (Additional file 1: Fig. S12b). Similar to that observed in cmt3b egg, relatively few upregulated genes were TEGs or showed hypo mCHG in cmt3b zygote (Additional file 1: Fig. S12c).

CMT3b represses zygotic genes in egg cells

In wild type zygote we identified 1804 down- and 2628 upregulated (> two-fold, Q < 0.01) genes relative to egg (Fig. 6a). Among the upregulated genes, 416 overlapped with previously identified genes expressed in rice zygotes (Fig. 6b) [31]. The 416 genes were enriched for chromatin replication and cell division functions (Fig. 6c), consistent with zygotic genome activation to promote zygote cell division in plants [32]. Within the 416 zygotic genes, 83 were upregulated in cmt3b egg. By contrast, although the cmt3a mutation caused a larger number of upregulated genes in egg, only 39 were among the 416 zygotic genes (Fig. 6d). Among the 83 zygotic genes upregulated in cmt3b egg were those encode histones, chromatin proteins (HMGs, SMC2, TOP2), DNA methyltransferases (MET1 and CMT3a), transcription factors (E2F, HAP3, NAC), and cell division-related proteins (Fig. 6e, f, g, Additional file 7: Table S6). The analysis suggested that CMT3b has a function to repress zygote gene expression program in egg.

Fig. 6

The cmt3b mutation de-represses zygote-expressed genes in egg. a Number of differentially expressed genes in wild type zygote relative to egg (FC > 2, p < 0.01). b Overlaps of high zygotic expression gene relative to egg in NIP, DJ and ZH11 varieties (Anderson et al., 2017; Zhou et al., 2021). c GO enrichment of zygote-expressed genes detected in three cultivars. d Transcript heatmaps of the 416 zygote-expressed genes in cmt3b egg (3b-E) and zygote (3b-Z) compared with wild type egg (E) and zygote (Z). e Transcript heatmaps of several representative zygote genes upregulated in cmt3b egg. f GO enrichment of the 83 upregulated genes in cmt3b egg shown in (d). g Integrative Genomics Viewer screenshots of six examples of the genes described in (d)