Nuclear m6A reader YTHDC1 regulates the scaffold function of LINE1 RNA in mouse ESCs and early embryos

Abstract

N⁶-methyladenosine (m⁶A) on chromosome-associated regulatory RNAs (carRNAs), including repeat RNAs, plays important roles in tuning the chromatin state and transcription, but the intrinsic mechanism remains unclear. Here, we report that YTHDC1 plays indispensable roles in the self-renewal and differentiation potency of mouse embryonic stem cells (ESCs), which highly depends on the m⁶A-binding ability. Ythdc1 is required for sufficient rRNA synthesis and repression of the 2-cell (2C) transcriptional program in ESCs, which recapitulates the transcriptome regulation by the LINE1 scaffold. Detailed analyses revealed that YTHDC1 recognizes m⁶A on LINE1 RNAs in the nucleus and regulates the formation of the LINE1-NCL partnership and the chromatin recruitment of KAP1. Moreover, the establishment of H3K9me3 on 2C-related retrotransposons is interrupted in Ythdc1-depleted ESCs and inner cell mass (ICM) cells, which consequently increases the transcriptional activities. Our study reveals a role of m⁶A in regulating the RNA scaffold, providing a new model for the RNA-chromatin cross-talk.

Keywords: 2-cell; H3K9me3; LINE1-scaffold complex; YTHDC1; retrotransposons.

Figure 1

YTHDC1 is essential for mouse ESCs. (A) Strategy for functional studies of Ythdc1 in ESCs. All cell lines were treated with 4-OHT for 3 days before harvest to ensure the depletion of endogenous Ythdc1. (B) Schematic of mouse wild-type (WT) YTHDC1, truncated YTHDC1 after the recombination (KO YTHDC1) and mutant YTHDC1 (W378A YTHDC1). aa, amino acid. (C) Growth curve showing that Ythdc1 cKO and W378A ESCs exhibited a poor proliferation rate. Cell numbers on the last day were used to assess the significance. (D and E) Colony formation abilities of Ythdc1 cKO and W378A ESCs were impaired revealed by AP staining. (F) RT-qPCR analysis showing the relative RNA level of key pluripotent markers in Ythdc1 f/f and cKO ESCs. (G) RT-qPCR analysis showing that EBs derived from Ythdc1 cKO ESCs exhibited the abnormal expression level of differentiation markers 7 days after in vitro differentiation. (H) Ythdc1 cKO and W378A ESCs exhibited a weak ability to generate chimeric mice. (I) Principal component analysis (PCA) showing the transcriptome differences between each ESC line. (J) GO analysis of genes dysregulated in both Ythdc1 cKO and W378A ESCs defined in Fig. S2E. Fold enrichment of each term is labeled in the plot. Data are presented as means with SDs (n = 3 in (C, F and G) and n = 4 in (D). Significance was calculated with unpaired two-tailed Student’s t test (**P < 0.01, ***P < 0.001, ****P < 0.0001). See also Figs. S1 and S2

Figure 2

Behaviors of YTHDC1-targeted RNAs in ESCs. (A) YTHDC1 RIP-seq signal is enriched around the center of m⁶A peaks defined in the control ESCs. (B) YTHDC1 RIP peaks preferred to localize at introns and retrotransposons in genome in the control ESCs. PROMPT, promoter upstream transcript. (C) Usage of exons (calculated as PSI, see METHODS) adjacent to YTHDC1 RIP peaks was decreased upon Ythdc1 depletion. (D) Expression level of YTHDC1-targeted coding genes is relatively lower in Ythdc1 deficient ESCs. A total of 941 coding genes possessing RIP peaks (fold enrichment > 5) in their exons or introns were defined as YTHDC1 RIP targets. (E) Venn diagram showing the overlap of YTHDC1 RIP targets and dysregulated genes defined in Fig. S2E. (F) GO analysis of YTHDC1-targeted coding genes defined by RIP-seq in the control ESCs. Fold enrichment of each term is labeled in the plot. (G and H) YTHDC1 RIP-qPCR and m⁶A MeRIP-qPCR analysis show the preference of YTHDC1 binding and m⁶A deposition on nuclear LINE1 transcripts in the control ESCs, respectively. Relative enrichment was calculated as the percent of input relative to the negative control antibody IgG. Two independent reactions for each IP were performed. Data are presented as means with SDs (n = 3 technical replicates). Primers used in these assays target the consensus sequences of LINE1 ORF regions (Su et al., 2012). We have checked the m⁶A signal on 18S rRNA to ensure that 18S could be used as a negative control using the published nuclear RNA m⁶A MeRIP-seq data (Liu et al., 2020). (I) Profiles generated by deepTools showing the distribution of m⁶A IP signal on bodies of METTL3-sensitive m⁶A-marked (left) and METTL3-insensitive m⁶A-marked (right) LINE1 caRNAs in WT and Mettl3 KO ESCs. M3, METTL3 or Mettl3. (J) Density plots and boxplots showing the occurrence frequency of RRACH (left) and ABAG (right) motifs in METTL3-sensitive /insensitive m⁶A peaks. M3, METTL3. (K) YTHDC1 RIP peaks were enriched on both METTL3-sensitive and METTL3-insensitive m⁶A-marked LINE1 transcripts in the nucleus. M3, METTL3. L1, LINE1. (L and M) METTL3-sensitive m⁶A-marked LINE1 RNAs increased their abundance on chromatin upon Ythdc1 depletion, whereas METTL3-insensitive m⁶A-marked LINE1 caRNAs remained unchanged in Ythdc1-depleted ESCs. Scatter plot showing the relative level of m⁶A-marked LINE1 caRNAs in the control ESCs (x-axis) and fold change of the level upon Ythdc1 depletion (y-axis) in L. A total of 1,989 METTL3-sensitive m⁶A-marked LINE1 transcripts, 1,664 METTL3-insensitive m⁶A-marked LINE1 transcripts (including 319 transcripts marked with both kinds of m⁶A) and 8,053 LINE1 transcripts without m⁶A sites were analyzed. M3, METTL3. caRNA m⁶A-seq data of WT and Mettl3^−/−-1 ESCs (GSE133600) published by Liu et al. (2020) were used for the analyses in (I, J, K, L and M). Input data of caRNA m⁶A-seq in control and Ythdc1 cKO ESCs were also downloaded from GSE133600 to quantify the relative abundance of LINE1 transcripts on chromatin in (L and M). Means of replicates were used to generate the summarized data in (C, D, L and M). Significance (**P < 0.01, ***P < 0.001, ****P < 0.0001) was calculated with two-tailed Student’s t test (paired in (C, D and M) and unpaired in (G, H and J)). See also Fig. S3

Figure 3

YTHDC1 is involved in the function of LINE1-scaffold complex in ESCs. (A) DEGs defined in LINE KD ESCs were similarly dysregulated in Ythdc1 cKO and W378A ESCs. 498 LINE1 KD-upregulated genes and 248 LINE1 KD-downregulated genes with |log₂(fold change)| > 0 and P value < 0.01 were defined by edgeR function in R using the published RNA-seq data (Percharde et al., ; GSE100939). Down, LINE1 KD-downregulated genes. Up, LINE1 KD-upregulated genes. (B) GSEA showing the global upregulation of LINE1 RNA-targeted genes (left, n = 2,397) and LINE1 sequence-enriched genes (right, n = 1,480) in Ythdc1 cKO ESCs. These genes were defined by Lu et al. (2020). NES, normalized enrichment score. (C) Many 2C-related genes were consistently upregulated in LINE1 KD ESCs and Ythdc1 cKO ESCs. A total of 76 genes with log₂(fold change) > 1 (2C::tdTomato⁺ versus 2C::tdTomato⁻) were considered 2C-related genes based on the published RNA-seq data (Macfarlan et al., 2012). RC ASO, reverse complement of LINE1 ASO (served as a negative control). L1 ASO, LINE1 RNA KD by ASO. (D) 2C-related genes (Macfarlan et al., ; defined as in (C) were significantly upregulated in Ythdc1 cKO ESCs but modestly downregulated in Mettl3 KO ESCs. (E) RT-qPCR analysis showing that representative 2C-related retrotransposons were upregulated in Ythdc1 cKO and W378A ESCs. (F) Scatter plot showing the relative abundance of retrotransposon subfamilies in Ythdc1 f/f ESCs (x-axis) and the expression fold change upon Ythdc1 depletion in ESCs (y-axis). 61 subfamilies with log₂(fold change) > 0.3 and P value < 0.05 were defined as Ythdc1 cKO-upregulated retrotransposons, and 304 subfamilies with |log₂(fold change)| < 0.1 were defined as Ythdc1 cKO-unchanged retrotransposons using edgeR function in R. Representative 2C-related retrotransposons consistently upregulated in LINE1 KD ESCs and Ythdc1 cKO ESC are labeled in the plot. (G) 2C-related retrotransposons were upregulated upon Ythdc1 depletion but not affected by Mettl3 depletion in ESCs. A total of 46 2C-related retrotransposons were identified as retrotransposons with log₂(fold change) > 1 (2C::tdTomato⁺ versus 2C::tdTomato⁻) based on the published study (Macfarlan et al., 2012), and other 565 repeat subfamilies serve as the control. (H and I) Proportion of 2C-like cells was significantly increased after Ythdc1 depletion in the ESC population. Fluorescence signal plot showing GFP signal distribution (x-axis) in a representative Ythdc1 cKO line with the Zscan4c promoter-GFP transgene in (I), and 2C-like cell population was determined according to the GFP signal in a WT control line without the transgene. (J) Co-IP in mouse ESCs showing that endogenous YTHDC1 interacted with both NCL and KAP1 proteins. DC1, YTHDC1. (K) NCL RIP-qPCR analysis showing that the association of NCL with nuclear LINE1 RNAs was attenuated in Ythdc1 cKO ESCs. Relative enrichment was calculated as the percent of input relative to the negative control antibody IgG. Two independent reactions for each IP were performed. Data are presented as means with SDs in (E, H and K) (n = 3). Means of replicates were used to generate the summarized data in (A, D, F and G). Significance (*P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001) was calculated with two-tailed Student’s t test (paired in (A, D and H) and unpaired in (E, G and K)). See also Fig. S4

Figure 4

YTHDC1 facilitates the genome recruitment of KAP1 in ESCs. (A) Scatter plot showing the KAP1 ChIP signal in Ythdc1 f/f ESCs (x-axis) and fold change of the signal upon Ythdc1 depletion in ESCs (y-axis) on 39807 KAP1 peaks identified in f/f ESCs. A total of 30217 Ythdc1 cKO-decreased KAP1 peaks were defined as KAP1 peaks with log₂(fold change) < −0.2. (B) KAP1 ChIP signal was significantly decreased on Ythdc1 cKO-decreased KAP1 peaks (defined in (A)) in Ythdc1 cKO and W378A ESCs. (C) Ythdc1 cKO-decreased KAP1 peaks (defined in (A)) were enriched at LINE1 RNA-targeted loci (Lu et al., 2020) and retrotransposons in genome. (D and E) RNA level was upregulated and KAP1 ChIP signal was decreased at LINE1 RNA-targeted loci (24202 regions defined by Lu et al., 2020) in Ythdc1 cKO and W378A ESCs. An equal number of random loci serve as the control. (F) Violin plot showing that Ythdc1 cKO-decreased KAP1 peaks (defined in (A)) were more enriched on Ythdc1 cKO-upregulated retrotransposons compared to Ythdc1 cKO-unchanged retrotransposons in ESCs. Ythdc1 cKO-upregulated/unchanged retrotransposons were defined as in Fig. 3F. (G) KAP1 ChIP signal was significantly decreased on Ythdc1 cKO-upregulated retrotransposons (defined as in Fig. 3F) in Ythdc1 cKO and W378A ESCs. (H) Scatter plot showing fold change of the retrotransposon level upon Ythdc1 depletion (x-axis) or Kap1 KD (y-axis) in ESCs. Representative 2C-related retrotransposons consistently upregulated in Ythdc1 cKO ESCs and Kap1 KD ESCs are labeled in the plot. Pearson’s correlation coefficient was calculated in R. shCtrl, control shRNA expressed. shKap1, Kap1 shRNA expressed. (I) Ythdc1 cKO-upregulated retrotransposons were significantly more derepressed compared to Ythdc1 cKO-unchanged retrotransposons upon Kap1 KD in ESCs. Ythdc1 cKO-upregulated/unchanged retrotransposons were defined as in Fig. 3F. (J) KAP1 ChIP-qPCR analysis showing the unaffected KAP1 binding at 5′ end of Dux locus, 3′ exons of Zfp11/Zfp13 genes and MERVL elements, and the significantly decreased KAP1 binding at certain 2C-related retrotransposons upon Ythdc1 depletion in ESCs. TSS of Rpl3 gene which is not targeted by KAP1 serves as a negative control. Relative enrichment was calculated as the percent of input normalized to enrichment at the negative control intergenic chromosome 11 (int-chr11). Two independent reactions for each IP were performed. Data are presented as means with SDs (n = 3 technical replicates). Means of replicates were used to generate the summarized data in (A, B, D, E and G–I). Significance (*P < 0.05, **P < 0.01, ****P < 0.0001) was calculated with two-tailed Student’s t test (paired in (B, D, E and G) and unpaired in (F, I and J)). See also Fig. S5

Figure 5

YTHDC1 promotes H3K9me3 installation to regulate the repression of 2C retrotransposons in mouse ESCs and early embryos. (A) Scatter plot showing fold change of KAP1 ChIP signal (x-axis) and H3K9me3 ChIP signal (y-axis) upon Ythdc1 depletion in ESCs on KAP1 peaks identified in Ythdc1 f/f ESCs. Pearson’s correlation was calculated in R. (B) H3K9me3 ChIP signal was more severely decreased on Ythdc1 cKO-decreased KAP1 peaks (defined in Fig. 4A) upon Ythdc1 depletion in ESCs. (C) H3K9me3 ChIP signal was significantly decreased on LINE1 RNA-targeted loci (24,202 regions defined by Lu et al., 2020) in Ythdc1 cKO and W378A ESCs. An equal number of random loci serve as the control. (D) H3K9me3 ChIP signal was significantly decreased on Ythdc1 cKO-upregulated retrotransposons (defined as in Fig. 3F) in Ythdc1 cKO and W378A ESCs. (E) Scatter plot showing fold change of MMETn RNA level upon Ythdc1 depletion (x-axis) or Kap1 KD (y-axis) in ESCs. A total of 225 expressed MMETn elements are included in this plot. (F) Scatter plot showing fold change of KAP1 ChIP signal (x-axis) and H3K9me3 ChIP signal (y-axis) upon Ythdc1 depletion in ESCs on genome loci of MMETn elements. A total of 225 expressed MMETn elements are included in this plot. (G) H3K9me3 ChIP-qPCR analysis showing the significantly decreased H3K9me3 level at certain 2C-related retrotransposons including MMETn, ETnERV3 and MuRRS elements and the unaffected H3K9me3 level at MERVL elements in Ythdc1 cKO ESCs. TSS of Rpl3 gene which lacks H3K9me3 marks serves as a negative control. Relative enrichment was calculated as the percent of input normalized to enrichment at the negative control intergenic chromosome 11 (int-chr11). Two independent reactions for each IP were performed. (H) Proportion of hatched blastocysts was decreased at E4.5 upon Ythdc1 depletion in mouse embryos. Blastocysts with embryonic cells moving out of the zona pellucida were considered as having hatched. (I) Scatter plot showing the relative abundance of retrotransposon subfamilies in Ctrl ICM (x-axis) and fold change of the level upon Ythdc1 depletion in ICM (y-axis). 33 subfamilies with log₂(fold change) > 0.3 and P value < 0.05 were defined as retrotransposons upregulated in Ythdc1 KO ICM using edgeR function in R. Representative 2C-related retrotransposons consistently upregulated in Ythdc1 cKO ESCs and Ythdc1 KO ICM are labeled in the plot. (J) Venn diagram showing the overlap of upregulated retrotransposons in Ythdc1 cKO ESCs (defined as in Fig. 3F) and in Ythdc1 KO ICM (defined as in (J)). Group C1, group C2 and group C3 retrotransposons are upregulated only in Ythdc1 cKO ESCs, in both Ythdc1 cKO ESCs and Ythdc1 KO ICM, and only in Ythdc1 KO ICM, respectively. Hypergeometric test was performed to calculate the significance of the overlap with phyper function in R. (K) RT-qPCR analysis showing that representative 2C-related retrotransposons were upregulated in E4.5 Ythdc1 KO blastocysts. (L) Ythdc1 KO ICM-decreased H3K9me3 peaks (defined in Fig. S6L) were enriched at LINE1 RNA-targeted loci (Lu et al., 2020) and retrotransposons. (M) Violin plot showing that Ythdc1 KO ICM-decreased H3K9me3 peaks (defined in Fig. S6L) were more enriched at group C2 retrotransposons compared to group C1 and C2 retrotransposons. Group C1-C3 retrotransposons were defined in (J). (N) H3K9me3 ChIP signal was more severely decreased on group C2 retrotransposons compared to group C1 and C2 retrotransposons in Ythdc1 KO ICM. Group C1-C3 retrotransposons were defined in (J). (O) Model showing the indispensable role of the m⁶A reader YTHDC1 in regulating the function of nuclear LINE1 scaffold in pluripotent cells. The LINE1 RNA scaffold is located at genomic loci of 2C-related retrotransposons and recognized by YTHDC1 through m⁶A modifications. YTHDC1 further promotes the recruitment of NCL-KAP1 and facilitates the deposition of H3K9me3, ensuring that these sites remain at a silent state and inhibiting the activation of 2C program. Data are presented as means with SDs in (G and K) (n = 3 technical replicates). Means of replicates were used to generate the summarized data in (A–F, I and N). Significance (*P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001) was calculated with two-tailed Student’s t test (paired in (C and D) and unpaired in (B, G, K, M and N)). See also Fig. S6