Intragenic DNA methylation modulates alternative splicing by recruiting MeCP2 to promote exon recognition

Abstract

Although the function of DNA methylation in gene promoter regions is well established in transcriptional repression, the function of the evolutionarily conserved widespread distribution of DNA methylation in gene body regions remains incompletely understood. Here, we show that DNA methylation is enriched in included alternatively spliced exons (ASEs), and that inhibition of DNA methylation results in aberrant splicing of ASEs. The methyl-CpG-binding protein MeCP2 is enriched in included ASEs, particularly those that are also highly methylated, and inhibition of DNA methylation disrupts specific targeting of MeCP2 to exons. Interestingly, ablation of MeCP2 results in increased histone acetylation and aberrant ASE-skipping events. We further show that inhibition of histone deacetylase (HDAC) activity leads to exon skipping that shows a highly significant degree of overlap with that caused by MeCP2 knockdown. Together, our data indicate that intragenic DNA methylation operates in exon definition to modulate alternative RNA splicing and can enhance exon recognition via recruitment of the multifunctional protein MeCP2, which thereby maintains local histone hypoacetylation through the subsequent recruitment of HDACs.

Figure 1

DNA methylation is enriched in included ASEs. (A) The cartoon on the left indicates the cassette exons: an ASE flanked by two constitutively spliced exons. “N1 transcript” represents the spliced RNA product that includes the ASE (red) and “N2 transcript” represents the spliced RNA product that excludes the ASE (blue). To detect the alternative splicing events, RNA was isolated from IMR90 or HCT116 cells and profiled using RNA-Seq. Numbers of included and excluded annotated ASEs of expressed genes, which were identified using MISO from the RNA-Seq data, are listed in the table. (B) DNA methylation level is significantly higher over included ASEs than excluded ASEs. CpG methylation levels were calculated from published data of IMR90 cells, within cassette exons from the two ASE groups described in A and the methylation levels were displayed in ± 300-bp windows relative to the splice acceptor (acc) and donor (don) sites across the cassette exons. The higher level of CpG methylation in included than excluded ASEs observed from this analysis is consistent with the higher level of CpG methylation in included exons demonstrated in an earlier report. A schematic representation of the ASE and flanking exons is displayed, oriented to the direction of transcription. The P-values for the relative enrichment of CpG methylation in included compared to excluded exons were calculated using one-sided t-test for the normalized DNA methylation level in the exons.

Figure 2

DNA methylation functions in pre-mRNA splicing. (A) DNA methylation is involved in splicing. To detect the function of DNA methylation in splicing, RNAs isolated from 5azadC-treated and control (CTR) IMR90 cells were profiled using RNA-Seq. Similarly, RNAs isolated from HCT116-WT or HCT116-DKO cells were profiled using RNA-Seq. The significantly changed exons were identified using MISO from the RNA-Seq data. The table lists the numbers of the significantly downregulated and upregulated exons induced by DNA methylation inhibition. (B) Hypermethylated exons have a higher tendency to be downregulated than hypomethylated exons following inhibition of DNA methylation by 5azadC treatment in IMR90 cells. The degree of exon downregulation was calculated and compared between hypermethylated and hypomethylated included exons. High meDNA indicates ASE with > 70% CpG methylation and low meDNA indicates ASE with < 30% CpG methylation. P-value was calculated by one-sided KS test. (C) Inhibition of DNA methylation leads to skipping of exon 10 of the HAUS8 gene in the spliced product. UCSC browser display of indicated data (NCBI36/hg18 assembly) over coordinates chr19:17,021,565-17,029,029 for HAUS8 focusing on the last three RefSeq annotated exons. DNA methylation data from HCT116-WT and IMR90 cells were obtained from public sources and displayed accordingly^,. Exon junctions from paired-end RNA-Seq data from 5azadC-treated IMR90 cells and HCT116-DKO cells and their respective control cells are displayed, where red indicates aberrant skipping events and black indicates normal “full-length” transcript. The number of tags representing each exon junction is shown to the left of the junctions. Exons are labeled according to RefSeq. meDNA: DNA methylation. (D) Comparison of DNA methylation in unperturbed (HCT116-WT) and demethylated (HCT116-DKO) cells in the HAUS8 exon 10 region. DNA methylation was determined by the bisulfite-sequencing of the HAUS8 exon 10 region. Each colored circle represents individual CpG cytosine methylation status within the analyzed region (black: methylated; white: unmethylated). Percentage of methylation was determined based on the number of methylated CpG cytosines divided by the total number of CpG cytosines within the analyzed sequences from all the PCR clones. P-value shown was calculated using one-sided Student's t-test.

Figure 3

MeCP2 exhibits DNA methylation-dependent binding to exons. (A) MeCP2 is enriched in promoters and exonic regions. MeCP2-binding profile in IMR90 cells was determined using ChIP-Seq and its distribution is displayed by normalizing to the genomic size of each indicated category. (B) High MeCP2 density correlates with high levels of DNA methylation in IMR90 cells. The MeCP2 densities at exons sorted according to DNA methylation levels as indicated are displayed. P-value was calculated by one-sided KS test. (C) Co-localization of MeCP2 and methylated DNA within exons occurs significantly more often than expected by random chance. The set of top n MeCP2-enriched exons shares k exons with the set of top n DNA-methylated exons. For each value n (x-axis), the overlap percentage R (y-axis) of the two sets is shown by a red curve. The overlap percentage expected by random chance is shown as a blue curve. (D) UCSC Genome Browser tracks showing MeCP2 ChIP-Seq data in HCT116 and IMR90 cells. The region shown contains exons 9-11 of the HAUS8 gene. (E) MeCP2 binding at the aberrantly skipped exon 10 of HAUS8 is significantly reduced upon DNA demethylation. MeCP2 binding was measured relative to input chromatin in the indicated cells by ChIP-qPCR. Error bars represent SE. P-values shown were calculated using one-sided Student's t-test.

Figure 4

MeCP2 is significantly enriched in included ASEs. (A) MeCP2 density is higher in included compared with excluded ASEs in IMR90 cells. The MeCP2 densities at the included and excluded ASEs and their flanking constitutively spliced exons in IMR90 cells were measured by ChIP-Seq and displayed in ± 300-bp surrounding the displayed exon. Schematic of the ASE and the flanking exons is displayed below each graph. (B) MeCP2 density is higher in included compared with excluded ASEs in HCT116 cells. Data are displayed as in A. All P-values displayed in A and B were calculated by one-sided KS test.

Figure 5

MeCP2 is critically involved in ASE inclusion in the spliced RNA products. (A) Significantly reduced MeCP2 mRNA levels in IMR90 cells upon MeCP2 knockdown. MeCP2 expression in IMR90 cells infected with viruses containing the negative control shRNA-luciferase (shLuc) construct or shRNA-MeCP2 (shMeCP2) construct was measured in triplicate using real-time qPCR and calculated as a percentage relative to β-actin mRNA levels using the delta-C_t method. Error bars represent SE. P-value shown was calculated using one-sided Student's t-test. (B) Reduced MeCP2 protein levels in IMR90 cells infected with viruses containing a shRNA-MeCP2 construct. MeCP2 protein expression in IMR90-shLuc and IMR90-shMeCP2 cells was measured by the western blot assays using the same antibody for ChIP-Seq experiments. β-Actin serves an internal control. (C) Knockdown of MeCP2 results in aberrant alternative splicing events. Schematic representation of aberrant exon downregulation and upregulation upon MeCP2 knockdown is displayed above the table listing the total numbers of these events in IMR 90 and HCT116 cells. For both cell lines, shLuc serves as a control shRNA construct, MeCP2-KD represents MeCP2 knockdown mediated by a specific shRNA. Higher number of aberrant ASE downregulation than upregulation events occur upon MeCP2 knockdown for both IMR90 and HCT116 (P < 10⁻¹⁰⁰, one-sided binomial test). (D) The original MeCP2 density is higher in the ASEs that become aberrantly downregulated than those that become upregulated upon MeCP2 knockdown. MeCP2 densities of the ASEs were measured in shLuc-transduced IMR90 cells by ChIP-Seq. (E) Nucleosome density is not significantly different between exons that become aberrantly upregulated and those that become downregulated upon MeCP2 knockdown. The average nucleosome densities for each ASE group were calculated based on data from IMR90 cells from the GEO database (GSE21823). (F) The original MeCP2 density in HCT116-WT cells is higher in ASEs that become aberrantly downregulated than those that become upregulated in HCT116-DKO cells. Average MeCP2 densities were calculated for each indicated group from ChIP-Seq data of HCT116-WT cells. P-values displayed in D-F were calculated by one-sided Wilcoxon test. (G) Ser2P-Pol II is significantly enriched in included ASEs. Ser2P-Pol II density is higher in included than excluded ASEs in HCT116 cells. The Ser2P-Pol II densities at included and excluded ASEs and their flanking constitutively spliced exons in HCT116 cells were measured by ChIP-Seq and displayed in ± 300 bp surrounding the displayed exon. Schematic of the ASE and flanking exons is displayed below each graph. All P-values displayed were calculated by one-sided Student's t-test.

Figure 6

Inhibition of HDACs results in aberrant ASE exclusion events that significantly overlap with those caused by MeCP2 knockdown. (A) Genome-wide alteration of splicing events is induced by TSA treatment. The total numbers of significantly downregulated (Down-reg.; aberrantly excluded) and upregulated (Up-reg.; aberrantly included) ASEs in TSA-treated HCT116-WT cells and TSA-treated IMR90 cells in comparison to their respective untreated control cells are displayed. There is significantly higher numbers of downregulated ASEs than upregulated ASEs in both cell types (P < 10⁻¹⁰⁰, one-sided binomial test). (B) Venn diagram showing the overlap between the aberrantly excluded ASEs upon MeCP2 knockdown and aberrantly excluded ASEs induced by TSA treatment in IMR90 cells. (C) Venn diagram showing the overlap between the aberrantly excluded ASEs upon MeCP2 knockdown and aberrantly excluded ASEs induced by TSA treatment in HCT116 cells. (D) Downregulated rather than upregulated ASEs induced by either MeCP2-knockdown or TSA treatment frequently occur together. RNA-Seq data of altered ASE expression from control and MeCP2-KD or TSA-treated HCT116 and IMR90 cells were compared. P-values were calculated using χ²-test.

Figure 7

MeCP2 recruits HDAC activity, promoting local histone hypoacetylation and exon inclusion. (A) Schematic representation showing the partial structure of the FAM204A gene. The red squares (a, b, and c) below the gene indicate the locations of qPCR probes used for ChIP assays. (B) Comparison of expression level of FAM204A isoforms that excludes (red) or includes (blue) exon 2. Expression of the FAM204A isoforms was determined by RNA-Seq data analyzed by MISO (confirmed by qPCR, data not shown). The normalized expression levels are displayed as a percentage, where the total expression of the two isoforms sum to 100%. WT, HCT116-WT cells; DKO, HCT116-DKO cells; TSA, TSA-treated HCT116-WT cells; MeCP2-KD, MeCP2-knockdown HCT116-WT cells. (C) Knockdown of MeCP2 decreased MeCP2 binding to the exon 1 (probe a) and ASE (probe b) regions of the FAM204A gene. MeCP2 occupancy at different DNA regions was measured using the ChIP assays with chromatin harvested from control (shLuc) or MeCP2-knockdown (MeCP2-KD) HCT116 cells, followed by quantification using the qPCR assays. (D) Knockdown of MeCP2 resulted in a significant increase of histone acetylation over the ASE region (probe b). Acetylation level was measured using ChIP-qPCR assays with chromatin harvested from the cells as described in C. (E) DNA demethylation by DNMT-deficiency resulted in decreased binding of MeCP2 to the ASE region (probe b). MeCP2 binding was measured using ChIP-qPCR assays with chromatin harvested from HCT116-WT and HCT116-DKO cells. (F) DNA demethylation by DNMT-deficiency resulted in increased histone acetylation over the ASE region (probe b). MeCP2 binding was measured using ChIP-qPCR assays with chromatin harvested from HCT116-WT and HCT116-DKO cells. (G) Inhibition of HDAC activity did not significantly change the binding of MeCP2 to the ASE region. MeCP2 binding was measured using ChIP-qPCR assays with chromatin from untreated (CTR) and TSA-treated HCT116-WT cells. Error bars in C-G represent SE. P-values shown in C-G were calculated using one-sided Student's t-test. N.S., not significant. (H) Cartoon showing the three independent treatments that downregulate the inclusion of ASEs: (1) inhibition of DNA methylation, (2) knockdown of MeCP2 and (3) inhibition of HDAC activity using TSA treatment. (I) Model depicting the relationship of DNA methylation, MeCP2 binding, HDAC recruitment and Pol II elongation in exon recognition.