Fine-tuning of gene expression through the Mettl3-Mettl14-Dnmt1 axis controls ESC differentiation

Abstract

The marking of DNA, histones, and RNA is central to gene expression regulation in development and disease. Recent evidence links N6-methyladenosine (m⁶A), installed on RNA by the METTL3-METTL14 methyltransferase complex, to histone modifications, but the link between m⁶A and DNA methylation remains scarcely explored. This study shows that METTL3-METTL14 recruits the DNA methyltransferase DNMT1 to chromatin for gene-body methylation. We identify a set of genes whose expression is fine-tuned by both gene-body 5mC, which promotes transcription, and m⁶A, which destabilizes transcripts. We demonstrate that METTL3-METTL14-dependent 5mC and m⁶A are both essential for the differentiation of embryonic stem cells into embryoid bodies and that the upregulation of key differentiation genes during early differentiation depends on the dynamic balance between increased 5mC and decreased m⁶A. Our findings add a surprising dimension to our understanding of how epigenetics and epitranscriptomics combine to regulate gene expression and impact development and likely other biological processes.

Keywords: DNA methylation; DNMT1; ESCs; METTL14; METTL3; differentiation; epigenetics; epitranscriptomics; gene expression; m(6)A.

Figure 1.. The DNA methyltransferase DNMT1 is recruited to chromatin by METTL3-METTL14

(A) NanoBRET energy transfer indicates that NanoLuciferase-METTL3 and -METTL14 are in close proximity to Halo-DNMT1 in living HEK293 cells (HaloTag as negative control, n = 4). (B) Immunostaining showing co-localization of transiently expressed DNMT1 with METTL3 or METTL14 in the nuclei of COS-7 cells. Scale bars: 5 μm, horizontal bar: median. (C) Endogenous co-immunoprecipitation (coIP) of DNMT1 with METTL3 and METTL14 in HeLa cells (n = 3). (D) In vitro coIP shows increased DNMT1-METTL14 interaction with rising METTL3 levels (n = 2). (E) Direct interaction between recombinant METTL14 and DNMT1, observed by in vitro pull-down followed by western blotting (n = 3). (F and G) Recruitment of DNMT1 to chromatin assessed by ChIP-qPCR in the GAL4–5XUAS system. In HEK2935XUAS cells, DNMT1 is recruited by the METTL14 RGG domain or by full-length METTL14, but not by the NLS domain nor by METTL14-ΔRGG2 (n = 4). (H) RNase treatment does not impair the interaction of METTL14-FLAG with DNMT1-Myc, as assessed co-immunoprecipitation in HeLa cells (n = 3). (I) Recruitment of DNMT1 to chromatin by METTL14 is independent of m⁶A activity (n = 3). ChIP-qPCR in the GAL4–5XUAS system performed as in Figure 1G but using a GAL4-tagged METTL14 mutant (R254/255A) unable to support METTL3-mediated m⁶A methylation. All data are means ± SEM (A) or means ± SD (F, G, and I). p values by two-way ANOVA (A), Pearson’s correlation analysis (B), and two-tailed unpaired t test (F, G, and I). See also Figure S1.

Figure 2.. METTL3 facilitates DNMT1-dependent gene-body DNA methylation

(A) Reduced 5mC levels (M value) in METTL3 KO HeLa cells, based on the 50,000 most variable CpGs (n = 3). (B) Volcano plot of 5mC changes in METTL3 KO HeLa cells, with significantly hypo- and hypermethylated sites as light and dark blue dots (corrected p value < 0.05 and Δβ value > 0.2). (C) Percentage of gene-body 5mC changes in METTL3 KO HeLa cells (measured as in Figure 2B). (D) Dot blotting in HeLa cells treated with METTL3 inhibitor STM-2457 (1, 5, and 25 μM) shows no difference in 5mC levels (n = 3; ns: not significant). Data (quantified for 25 μM) as means ± SEM. (E and F) Volcano plot showing no changes in 5mC following METTL3 (STM-2457; E) or FTO (FB23–2; F) inhibitor treatment (n = 3). (G and H) DNMT1 does not bind H3K36me3 by coIP in HeLa cells transiently expressing DNMT1-Myc (G) and in vitro pull-down of recombinant FLAG-tagged DNMT1 by unmodified and H3K36me3-modified nucleosomes (H). DNMT3B and recombinant glutathione S-transferase (GST)-PWWP domain of DNMT3B as positive controls (n = 3). (I) SETD2 knockdown in HeLa cells has very limited impact on METTL3-dependent gene-body DNA methylation (n = 3). (J) Knockdown of DNMT1, DNMT3A, or DNMT3B in HeLa cells reveals that METTL3-dependent gene-body 5mC mostly relies on DNMT1 (n = 3). (K) METTL3 KO impairs DNMT1 binding to chromatin, but not that of DNMT3A or DNMT3B (n = 3). (L) DNA methylation changes among METTL3-dependent 5mC genes in HeLa cells after knockdown of DNMT1, DNMT3A/DNMT3B, or SETD2, or perturbation of m⁶A (by METTL3i or FTOi). p values by two-tailed t test (A and D). See also Figure S2.

Figure 3.. 5mC and m⁶A contribute together to regulating the expression of common target genes

(A) In HeLa cells, gene-body 5mC strongly co-occurs with coding-sequence m⁶A. (B) Venn diagram of the overlap between 5mC-marked gene-bodies (mean β value > 0.25) and m⁶A-marked coding sequences (peak in m⁶A-seq). (C) Proportion of m⁶A-associated CpGs by 5mC level (0%–25%, 25%–50%, 50%–75%, 75%–100% β value), determined in the TSS and CDS by bootstrapping (STAR Methods). (D) Association between 5mC and m⁶A in TSS (top) and CDS (bottom). CpGs grouped by 5mC levels (0%–25%, 25%–50%, 50%–75%, 75%–100% β value) and by m⁶A association or lack thereof. (E) 5mC levels of the top 2,000 most variable intragenic CpGs among 5mC-m⁶A target genes show dependence on METTL3 and DNMT1. (F) Percentage of gene-body CpGs within 5mC-m⁶A targets with decreased or increased 5mC in METTL3 KO HeLa cells. (G) Genes exhibiting hypomethylation (Δβ < −0.2) upon METTL3 KO are globally downregulated by RNA-seq (n = 3). Two-tailed Wilcoxon test. (H and I) Genes showing reduced m⁶A (1.5-fold decrease) in METTL3 KO show increased gene expression, as determined by IP and input of m⁶A-seq, respectively (n = 2). Spearman correlation analysis and two-tailed t test. See also Figure S3.

Figure 4.. In ESCs, Mettl3-Mettl14 partner with Dnmt1 for DNA methylation deposition in gene-bodies

(A) Reduced 5mC levels in Mettl3 KO ESCs (n = 3). (B) 5mC changes (Δβ value) in Mettl3 KO ESCs, with hypo- and hypermethylated sites in light and dark blue (corrected p value < 0.05 and Δβ > 0.2). (C) Reduced 5mC levels in gene-bodies (corrected p value < 0.05 and Δβ > 0.2). (D) Gene-bodies whose 5mC marking is Mettl3-dependent (defined as 5mC loss ≥10% in Mettl3 KO) are barely affected in Dnmt3b KO ESCs. From published data (GEO: GSE72856). (E) High signal in Mettl14 ChIP-seq, but not in H3K36me3 ChIP-seq among sites of Mettl3-dependent 5mC deposition. From published data (GEO: GSE206730 and GSE31039). (F) Reduced 5mC across all genes in Dnmt1 KO ESCs (n = 3). (G) Targets of Mettl3-mediated 5mC deposition depend on Dnmt1 for gene-body methylation. (H) Subcellular protein fractionation indicates that Mettl3 knockout, but not catalytic inhibition, reduces Dnmt1 binding to chromatin in ESCs (n = 2). (I) Chart showing, for genes with METTL3-dependent gene-body 5mC (total column height), the proportion (in blue) of genes showing Dnmt1 and/or Dnmt3b dependency (defined as 5mC loss ≥10% in the corresponding DNA methyltransferase knockout) or genes associated with Mettl14 binding to H3K36me3 (Mettl14 and H3K36me3 ChIP-seq overlap). (J) Dnmt1 binding (by ChIP-seq), Dnmt1 largely overlaps with Mettl3 and Mettl14 binding in genes with “Mettl3-dependent” gene-body 5mC. From published data (GEO: GSM2059182, GSE202848, and GSE206735). (K) Among genes with “Mettl3-dependent” gene-body 5mC, there is a strong overlap of Mettl3-Mettl14 ChIP-seq targets with Dnmt1, but minimal overlap with Dnmt3b. p values by two-tailed unpaired t test (A) and one-sided hypergeometric test (J). See also Figure S4.

Figure 5.. The transcriptional effect of gene-body 5mC and the post-transcriptional effect of m⁶A combine to regulate gene expression

(A) Distribution of 5mC sites (TSS vs. gene-body), according to corresponding transcript m⁶A status. (B) Strong overlap between genes with m⁶A-marked transcripts and gene-body 5mC (mean β value > 0.25). (C) Distribution of m⁶A peaks (TSS vs. CDS) in ESCs (left) and proportions of 5mC-marked and -unmarked CpGs with m⁶A-marked transcripts (right). (D) SLAM-seq in Dnmt1 KO ESCs highlights changes in nascent transcription (n = 3), with significant up- and downregulation in red and blue, respectively (fold-change > 1.5 and corrected p value < 0.05). (E) Decreased gene-body methylation in Dnmt1 KO ESCs coincides with reduced nascent transcription and reduced steady-state transcript levels. (F) SLAM-seq in ESCs treated with 50 μM STM-2457 shows that while inhibition of Mettl3 catalytic activity increases steady-state levels of m⁶A-marked RNAs, it does not affect nascent transcript formation (n = 3). (G) Depletion of Mettl3 (from GEO: GSE86336) or Mettl14 (in-house data, n = 3) increases the stability of transcripts that are normally m⁶A-marked (by actinomycin D assay followed by RNA-seq). (H) Comparison of nascent and steady-state transcript levels (SLAM-seq) for 5mC-m⁶A targets indicates that Dnmt1 KO mostly affects transcription, whereas the effect of Mettl3 inhibition is post-transcriptional (top, see Figure S5H). Post-transcriptional regulation can be quantified by the difference between steady-state and nascent RNA levels (bottom). (I) Mettl3 inhibition and KO display similar post-transcriptional effects (by SLAM-seq and actinomycin D assay) on 5mC-m⁶A targets but different steady-state regulation (RNA-seq). (J) Gene expression (input m⁶A-seq), m⁶A (m⁶A-seq), and 5mC (Infinium array) were tracked during the transition from naive to formative pluripotency (at 0, 3, and 24 h, n = 3). (K) Increased 5mC (Δβ > 0.1) is associated with gene upregulation 24 h after induction of formative pluripotency (left). Precision nuclear run-on sequencing (PRO-seq) indicates a concomitant rise in active transcription (right, n = 3). p values by chi-squared test (A), hypergeometric test (C), two-way t test (E, F, I, and K), Kolmogorov-Smirnov test (G). See also Figure S5.

Figure 6.. 5mC and m⁶A are both required during ES-to-EB differentiation

(A) Schematic model of ES-to-EB differentiation upon LIF removal, with all ESC lines used. (B) Immunostaining of differentiation markers (Sox17, Gata4, Gata6, Foxa2) indicates that impaired EB formation in Mettl14 KO cells is fully rescued by expression of Mettl14^WT but only partially by Mettl14^ΔRGG expression. Scale bars, 100 μm. (C) Mettl14 KO reduces 5mC levels in EBs, and these are rescued by expression of Mettl14^WT, but not Mettl14^ΔRGG (n = 3). (D and E) Altered 5mC levels in Mettl3 KO EBs (D), with hypo- and hypermethylated CpGs in light and dark blue (corrected p value < 0.05 and Δβ > 0.2), and heatmap of the top 2,000 most variable intragenic CpGs (E). (F) PCA plot for all m⁶A peaks found by m⁶A-seq, showing that EBs and ESCs have distinct m⁶A profiles (n = 3). (G) Changes in m⁶A between WT EBs and ESCs, with significant changes highlighted in orange (left, fold-change > 1.5 and corrected p value < 0.05) and number of significant peaks/genes (right). (H) Heatmap of the top 2,000 most variable m⁶A regions during the ESC-to-EB transition. (I) Hierarchical clustering showing that EBs lacking Mettl3 or Mettl14 cluster with Mettl14^ΔRGG-rescued EBs, from RNA-seq (n = 3). (J and K) Comparing Mettl14 KO with Mettl14^WT and Mettl14^ΔRGG allows to assess the effects of reducing both 5mC and m⁶A vs. m⁶A alone in EBs. Gene set enrichment analysis (GSEA) of gene expression indicates that Mettl14 KO strongly impairs differentiation, with Mettl14^ΔRGG showing an intermediate phenotype. See also Figure S6.

Figure 7.. Dynamic, coordinated adjustments of 5mC and m⁶A fine-tune expression of key genes in differentiating ESCs

(A and B) Identification of 1,732 candidate genes (see Figure S7A and STAR Methods) and changes in expression upon Mettl3 inhibition (Exp: steady-state, Tx: nascent transcription, by SLAM-seq) or knockout (by RNA-seq). (C) Expression of candidate genes in ESCs and EBs in WT and Mettl3 KO conditions. (D) Heatmap showing that candidate genes upregulated in EBs (top 500 by EB/ES fold-change in WT cells, referred to as “up” candidates) show impaired expression in Mettl3 KO EBs. (E and F) Divergent regulation of differentiation-involved candidate genes following Mettl3 inhibition or Mettl3 KO in ESCs by RNA-seq (E) and validation by RT-qPCR (F, n = 6). (G) Candidate genes reach higher gene-body 5mC levels than other genes during EB differentiation. (H) Gene-body methylation of candidate genes depends on Mettl3 and Dnmt1. (I) Dnmt1 KO reduces nascent transcription of candidate genes, especially for genes that are upregulated in WT EBs (“up,” as defined in Figure 7D). (J) “Up” candidate genes display a global loss of m⁶A during differentiation, with significant peaks highlighted in orange (fold-change > 1.5 and corrected p value < 0.05). (K) “Up” candidate genes display increased transcript half-lives in Mettl14 KO ESCs. (L and M) Comparing Mettl14 KO, Mettl14^WT, and Mettl14^ΔRGG conditions reveals the effects of reducing 5mC alone, m⁶A alone, or both marks in EBs. “Up” candidate genes display opposite effects for 5mC and m⁶A, with the influence of 5mC prevailing during differentiation. Data as mean ± SD (F). p values by two-way t test (B, F, G, and K), one-sample t test (I and L). See also Figure S7.