Quantitative analysis questions the role of MeCP2 as a global regulator of alternative splicing

Summary

MeCP2 is an abundant protein in mature nerve cells, where it binds to DNA sequences containing methylated cytosine. Mutations in the MECP2 gene cause the severe neurological disorder Rett syndrome (RTT), provoking intensive study of the underlying molecular mechanisms. Multiple functions have been proposed, one of which involves a regulatory role in splicing. Here we leverage the recent availability of high-quality transcriptomic data sets to probe quantitatively the potential influence of MeCP2 on alternative splicing. Using a variety of machine learning approaches that can capture both linear and non-linear associations, we show that widely different levels of MeCP2 have a minimal effect on alternative splicing in three different systems. Alternative splicing was also apparently indifferent to developmental changes in DNA methylation levels. Our results suggest that regulation of splicing is not a major function of MeCP2. They also highlight the importance of multi-variate quantitative analyses in the formulation of biological hypotheses.

Keywords:

DNA methylation – Gene expression – Gene regulation – Neurons – Developmental neuroscience – RNA sequencing – Alternative splicing – RNA splicing

Introduction

MeCP2 (methyl CpG binding protein 2) is an important mediator of epigenetic regulation in mammalian cells, and particularly in neurons [1–3]. Mutations that lead to depletion or overexpression of the MECP2 gene are associated with severe neurological diseases, most notably Rett syndrome, a devastating autism-spectrum disease affecting approximately 1 in 10,000 female births [4]. MeCP2 is a chromatin protein that binds to modified cytosine residues, primarily mCG [5], the most abundant modification in mammalian genomes, but also mCA, particularly mCAC [6], which is a feature of neuronal cells [7, 8]. This association with chromatin enables MeCP2 to interpret epigenetic signals to modulate gene expression events. By far the best documented role of MeCP2 is repression of transcription: several large-scale studies have confirmed this role, and proposed mechanistic models that can account for the effect of MeCP2 on transcription [9–11]. Other roles for MeCP2 have also been hypothesised, however, including chromatin compaction, micro-RNA processing, and regulation of alternative splicing (reviewed in reference [12]). These hypotheses highlight alternative roles of MeCP2 in orchestrating gene expression at transcriptional and post-transcriptional levels, but they have varying degrees of support from data.

In this paper, we focus on the possibility that MeCP2 is a regulator of alternative splicing. This was first proposed after a physical interaction with the splicing factor YB1 was observed [13] in 2005. Subsequent studies provided additional support; most recently via a modest but significant association between changes in DNA methylation and intron retention in the transition between promyelocytes and granulocytes [14]. A role of MeCP2 in splicing is plausible from a mechanistic perspective [9–11] if binding of the protein serves as a “brake” for polymerase, leading to a reduced transcription rate, which in turn can alter splicing preferences [15]. These studies relied on correlative, univariate analyses of data from experiments which were not designed to specifically test the association of MeCP2 with splicing. This renders the conclusions vulnerable to confounding effects, which could explain both changes in splicing and DNA methylation/MeCP2 occupancy. Although a measure of significance of the association (assuming no confounding factors) was deduced, the effect size was not quantified. As a result, questions regarding the extent to which variation in splicing can be attributed to changes in DNA methylation or MeCP2 levels remain open.

Here we revisit the evidence for this phenomenon, leveraging recent, high-quality data sets which track transcriptomic changes in experiments specifically designed to reflect changes in MeCP2 activity. These new data sets offer an unprecedented opportunity to tackle this question statistically, providing a quantitative estimate of MeCP2’s effect while carefully controlling for confounding effects.

Results

It has been suggested that DNA methylation regulates alternative splicing through different mechanisms [16]. We evaluated the relationship from two perspectives: either dependent on MeCP2 or independent of MeCP2. Considering the former first, we sought transcriptomic data sets from cells or tissues with different controlled levels of MeCP2 while keeping DNA methylation constant. Several comprehensive data sets have recently been collected in situations that closely mimic this idealised scenario (Table 1). In particular, our previous work [9] captured transcriptomic data sets in Lund Human Mesencephalic (LUHMES)-derived human dopaminergic neurons that expressed widely different levels of MeCP2, while total levels of DNA methylation remained constant. Additionally, Boxer et al. [17] comprehensively captured transcriptomes and methylomes from wildtype (WT) and MeCP2 (KO) mouse brain tissues. To investigate the latter, we examined the Dnmt1/3a/3b triple-knockout mouse embryonic stem cells (DNMT-TKO mESCs) that lack DNA methyltransferase activity [18, 19]. Yearim et al. performed splicing analyses on a relatively low-coverage data set and suggested that DNA methylation may influence splicing of alternate exons [18]. Therefore, we sought to revisit this conclusion using a more recent, higher coverage DNMT-TKO dataset [19] to investigate the role of DNA methylation in influencing alternative splicing independent of MeCP2. Finally, it is also possible that splicing regulation is due to a combination of changes in MeCP2 levels and DNA methylation. In early postnatal brain, the level of genomic mCA, a major determinant of MeCP2 binding, concurrently increases with the level of MeCP2 [5, 20, 21]. Stroud et al. [20] profiled the transcriptomes and methylomes of differentiating mouse neurons in the first few weeks after birth. Therefore, developing neurons can be considered an ideal scenario where one can investigate the combined effects of MeCP2 and DNA methylation. For completeness, we also re-analysed the non-neuronal data sets from promyelocytes and granulocytes [14] which first led to suggest the potential link between MeCP2, DNA methylation and alternative splicing. MeCP2 RNA expression levels show predominant expression in neurons and brain tissue while mESCs and granulocyte express significantly lower levels (S1 Fig). All of these data sets provide a good level of replication and high sequencing depth, with the exception of the early developmental time points of [20], where experimental constraints limited the achievable depth.

Adequate sequencing coverage is important for robust alternative splicing analysis, and some of the data sets we employ, particularly the early developmental times in the Stroud et al. data [20], have relatively low coverage. To address this, we use a recently proposed splicing quantification strategy, BRIE (Bayesian Regression for Isoform Estimation [22]), which uses Bayesian sequence-derived, informative priors to provide robust estimates of splicing ratios even at very low coverage levels such as encountered in single cell RNA-seq data.

Summary statistics for all the data sets used, including sequencing depth and replication level, are provided in Table 1.

Differential splicing analysis questions the function of MeCP2 and DNA methylation as a global regulator of alternative splicing

To explore the effects of MeCP2 on splicing and gene expression, we initially performed differential splicing and differential gene expression analyses for a number of data sets which investigated changes in MeCP2 levels (Table 1). Varying the abundance of MeCP2 in cultured human neurons over a wide range (absence (KO), 1x (WT), 4x and 11x wildtype (OE) levels) is accompanied by many changes in gene expression relative to wildtype (KO = 178; 4xOE = 2548; 11xOE = 2348 at a significance threshold of p-adjusted value < 0.05), but relatively few differential splicing events (KO = 101; 4xOE = 120; 11xOE = 73 at a significance threshold of Bayes factor ≥ 3) (S2A and S2B Fig). We observe lower number of differential splicing events in KO mouse brain (57) compared to changes in gene expression (2943) (S2C and S2D Fig). This suggests that changing the levels of MeCP2 (in vitro or in vivo) results in sizable changes in transcription but negligible changes in splicing.

In WT mouse neurons, where levels of DNA methylation and MeCP2 increase during development, we found that the number of differential splicing events (489, 408) is lower compared to the number of differential gene expression events (1122, 538) when 8 weeks old and 1 week old Parvalbumin (Pv) and Vasoactive intestinal peptide (Vip) neurons are compared (S2E and S2F Fig). Comparisons between 8 weeks old and 3 weeks old Pv and Vip neurons showed more differential splicing events (60, 527) than differentially regulated genes (4, 118) (S2E and S2F Fig). Based on this analysis, we could not establish whether MeCP2 and DNA methylation predominantly perturbs transcription or alternative splicing in developing neurons. Because of the increase in CA methylation [20, 21] and the dynamic nature of alternative splicing [24] in brain development, it is necessary to investigate the role of DNA methylation independent of developmental effects and MeCP2. Therefore, we performed similar analysis for DNMT-TKO mESCs in comparison to WT mESCs. In the absence of DNA methylation, many genes are differentially regulated (7342), but there are relatively few differential splicing events (134) (S2G and S2H Fig). Based on this analysis, DNA methylation predominantly perturbs transcription, with a much smaller effect on splicing. In a separate data set for DNMT-TKO mESCs [18], the number of differentially regulated genes (102) was only marginally higher than number of differential splicing events (71) (S2I and S2J Fig). Differential gene expression analysis was not possible on the RNA-seq data set of Promyelocytes and Granulocytes [14, 23] due to the absence of replicates. The lower number of differential splicing events compared to differentially expressed genes is consistent at different significance thresholds (S2 Fig) in all comparisons except in developing neurons. To investigate whether differential splicing events are directly associated with context specific DNA methylation or MeCP2 binding, we performed correlation analysis on significanctly altered events (S2B, S2D, S2F, S2H and S2J Fig). We found that there is no apparent linear relationship between changes in splicing and changes in DNA methylation or MeCP2 binding (S1 Table) except a modest one for Pv neurons (highlighted in S1 Table). By performing classical differential gene expression and differential splicing analysis, we find that MeCP2 primarily functions as a regulator of gene expression and has minimal effect on splicing. Whether the primary role of DNA methylation and neuronal development is in influencing gene expression or alternative splicing, is uncertain. To test this in a more quantitative manner, we used machine learning techniques in our subsequent multivariate analysis.

Sequence features are highly predictive of splicing ratios

Machine learning approaches have previously shown that genomic sequence features are highly predictive of exon inclusion/exclusion ratios (splicing ratios) and the splicing code assembled using these features helped characterise alternative splicing on a genome-wide scale [25, 26]. The splicing ratios were calculated using several independent methods, and our analysis indicated a strong correlation between results obtained with different methods (S3A and S3E Fig). Methods that directly model splicing ratios such as Miso [27] and BRIE [22] showed particularly strong agreement (S3A Fig). To quantify the explanatory power of sequence features in predicting splicing ratios in the current data sets, we used two independent multivariate regression models: Linear Regression and Random Forest Regression, the latter being a flexible non-parametric approach. We first validated that the relationship between the splicing code and splicing ratios [25, 26] could be confirmed for all the data sets described in Table 1 We find that both models explain 50%-65% of the variation (Fig 1A, 1B, 1D, 1E and 1F left panels). This effect is quantified by the metric fraction of variation explained (FVE) which is analogous to R². When using splicing ratios from MISO, we still retain 35%-50% of the explanatory power of sequence-derived features (Fig 1A, 1B, 1D, 1E and 1F right panels). We note that predictive models achieved lower accuracy in Vip neurons and Pv neurons during early mouse brain development (Fig 1C), possibly due to the highly dynamic nature of splicing in early development [24]. We confirm that splicing ratios can be accurately regressed against sequence-derived features and this analysis provided positive controls for further analysis.

MeCP2 changes account for a negligible fraction of splicing variation

We then probed whether changes in inclusion/exclusion ratios between different cell types could be explained by changing levels of MeCP2. We refer to differential exon inclusion/exclusion between different cell types as splicing variation (not to be confused with splicing ratios, which are the ratio of inclusive/exclusive isoform abundance within the same sample). Since global changes in transcription correlate with MeCP2 levels [9] and sequence-derived features can predict splicing ratios (Fig 1) [25, 26], our hypothesis is that any role of MeCP2 in regulating splicing should be reflected in improved predictions of differential exon inclusion/exclusion. Since global DNA methylation remains constant in neurons expressing multiple levels of MeCP2 (S4A Fig), we assume that DNA methylation at exons and introns also remain constant. To allow for the possibility that MeCP2 regulation of splicing might be dependent on the local sequence context, we used regression models with both MeCP2 levels and sequence features as input variables.

When we quantified the explanatory power of contrasting MeCP2 levels (in vitro) by regressing splicing variation against sequence features and MeCP2 binding, we could only explain 0-3% of the splicing variation regardless of the quantification method or model (Fig 2A left panel). In some instances, the model fit is extremely poor, as evidenced by negative levels of splicing variation explained (Fig 2A right panel). To recapitulate a similar situation in vivo, we used data sets from [17] which profiled transcriptomes and methylomes of WT and MeCP2 KO mouse brain tissues at very high sequencing depth. We found that DNA methylation levels, quantified from WGBS-seq, are statistically indistinguishable between WT and MeCP2 KO mice brains (S4B and S4C Fig): correlations between methylomes of different WT replicates are identical to correlations between methylomes in WT and KO samples. This further supports our assumption that DNA methylation remains constant in neurons expressing different levels of MeCP2. In mouse brain data sets, when we quantified the explanatory power of MeCP2 by regressing splicing variation against sequence features, the model fit was once again poor, with negligible or even negative levels of splicing variation explained (Fig 2B). Importantly, the 57 differentially spliced genes in this data set did not show any significant level of differential methylation (S1 Table), ruling out a potentially confounding effect of a compensatory local change in methylation. Therefore, multiple levels of MeCP2 (in vitro and in vivo) combined with sequence features cannot predict splicing variation, suggesting that MeCP2 does not regulate splicing, even when sequence context is taken into account.

DNA methylation changes account for a negligible fraction of splicing variation

MeCP2 binding in neurons is strongly associated with DNA methylation in a CG and CA context [6]. We asked whether changes in context specific DNA methylation in developing brain, the primary determinant of MeCP2 binding, affects splicing variation. To do so, we calculated the changes in DNA methylation at alternative exons and surrounding regions for both dinucleotides CG and CA (Fig A in S1 Appendix) in developing mouse neurons [20]. In this instance, splicing variation is the difference in splicing ratio across Pv and Vip neurons at different ages. We quantified the explanatory power of methylation features by regressing splicing variation against features derived from DNA methylation, as well as sequence features. The fraction of variance explained by the models is 0-10% (Fig 3A left panel). The analysis fails to capture any relationship (linear/non-linear) when splicing variation is regressed against DNA methylation alone (S5A Fig). Notice, however, that the inclusion of sequence features accounts for context-specific effects of DNA methylation, in a scenario where MeCP2 is essential for normal function of nerve cells. Therefore, context specific DNA methylation and sequence features in developing neurons are only at best modest predictors of splicing variation.

We then specifically set out to test whether changes in DNA methylation and their sequence context can predict splicing variation independently of MeCP2. To investigate the relationship between DNA methylation and alternative splicing independent of MeCP2 and developmental effects, we chose DNMT-TKO mESCs which lack DNA methylation [18, 19]. As these cells also express low levels of MeCP2 (S1 Fig), this analysis is important to rule out any confounding effects on our conclusions possibly caused by the high correlation between MeCP2 binding and DNA methylation in neuronal cells. Here, splicing variation is the difference in splicing ratio between DNMT-TKO mESCs and WT mESCs. We quantified the explanatory power of methylation features by regressing splicing variation against DNA methylation and associated DNA sequence features. Both regression models are able to account for a minimal fraction (0-10%) of the observed variance in splicing ratios (Fig 3C left panel). In a separate DNMT-TKO mESCs data set, the accuracy was lower (0-5%) (Fig 3B left panel). Therefore, DNA methylation when considered in a scenario independent of MeCP2 and developmental effects cannot predict splicing variation.

The promyelocytes-granulocytes data set is the only one where methylation levels combined with sequence features do appear to explain a small fraction of variance in splicing variation (5-10%) (Fig 3D). However, even in this case, the fraction of variance explained remains modest when compared to the explanatory power of sequence features in predicting splicing ratios (65% at best and 25% at worst, Fig 1 left panels). It is important to stress that both regression models perform worse when DNA methyation features alone are used to predict splicing variation (S3B Fig), rendering the initial observation moot [14]. In addition to the absence of replicates, as noted above, a further reason for exercising caution is that this study did not explicitly control for changes in MeCP2 abundance.

Discussion

MeCP2 is a major interpreter of the cell’s epigenetic state whose effects on transcription are extensively categorised and recognised as consistent in multiple systems [2, 6, 9–11]. In addition to its effects on gene expression, multiple other potential functions for MeCP2 have been proposed. However, its high abundance in neurons and pervasive presence in the chromatin environment create difficulties in attributing functions to MeCP2 based purely on observational studies. One often cited function of DNA methylation and its cognate interpreter, MeCP2, is regulation of alternative splicing. Here we have used multivariate analysis of recent data sets to examine this hypothesis in systems where either DNA methylation levels or MeCP2 levels are varied either alone or in combination. The results reveal that MeCP2’s role as a global regulator of RNA splicing is minimal and therefore of questionable biological significance.

While statistically our conclusion is robust, it must be qualified in several ways. First of all, we focused on detecting a global role in splicing for MeCP2. The fact that MeCP2 variation explains a negligible fraction of global splicing variation does not imply that MeCP2 plays no role in regulating splicing for specific transcripts. Our analysis would not detect a role of this kind, which would require more targeted experiments to characterise and validate any affected genes. It should however be pointed out that our analysis of splicing variation of WT and MeCP2 KO mouse brains highlighted that changes in splicing were not associated with any significant changes in DNA methylation at the same loci, suggesting that even a gene-specific role of MeCP2/DNA methylation in regulating splicing is likely to be complex and/ or indirect. A second important reservation is that the small splicing effects that we did detect are not necessarily unimportant. Diseases associated with MECP2 mutations are long-term pathologies, indicating that cells can survive for years in the presence of MeCP2 aberrations. Therefore, it is likely that the molecular effects of MeCP2 will be subtle, a point well-made previously [2, 6, 9–11]. Thirdly, it should be noted that our best-controlled data sets and most of our analyses apply to neuronal cells where MeCP2 is both abundant and strongly implicated in neuronal cell physiology. While it is appropriate that questions relating to MeCP2 function should focus on neuronal data sets, this might overlook effects in other cell types. Indeed, we did detect a modest statistical association between combined DNA methylation and DNA sequence features in the promyelocyte-granulocyte transition, although there was no relationship between DNA methylation alone and alternative splicing in this system. MeCP2 was initially implicated in regulation of alternative splicing in this myeloid lineage, but genetic depletion of MeCP2 specifically in non-brain tissues has no obvious phenotypic consequences, indicating that the protein is dispensable outside the central nervous system [28].

On a more general note, our analysis sounds a note of caution regarding the risks of over-interpretation of experiments that rely on high-throughput data. While such experiments are key in expanding our understanding of molecular physiology, univariate analyses can mask confounding factors. Hence, we argue that biological hypotheses arising from observational studies should always be further validated by carefully designed experiments and supported by multivariate statistical analyses.

Materials and methods

Sequencing data sets

Table 1 details the data sets used in the analysis. These include RNA sequencing (RNA-seq) and Bisulfite sequencing (BS-seq) libraries quantifying DNA methylation levels in Dnmt1/3a/3b triple-knockout (DNMT-TKO) and wildtype (WT) mouse embryonic stem cells (mESCs). Also, data sets concerning MeCP2 include RNA-seq and Chromatin immunoprecipitation followed by sequencing (ChIP-seq) libraries from post-mitotic (day 9 of differentiation) human neurons expressing multiple levels of MeCP2 and BS-seq libraries quantifying DNA methylation patterns from wildtype MeCP2-expressing neurons; RNA-Seq and BS-Seq libraries from mouse neuronal subtypes paravalbumin (Pv-) and vasoactive-intestinal-peptide (Vip-) expressing interneurons in mouse brain at different post-natal ages; RNA-seq and BS-seq libraries from cells of the mouse myeloid lineage, namely Granulocytes and Promyelocytes. Information about genome assemblies, annotation, software versions and command-line arguments are detailed in S1 Appendix.

Quantification of splicing ratios

RNA-Seq data sets are downloaded and sequencing reads are mapped to the reference genome using HISAT2 [29] aligner. Splicing ratios are quantified using Miso [27], BRIE [22] and Cufflinks [30]. We found Cufflinks to perform poorly in simulated data with known Ψ values (S3D Fig) when the sequence reads were single-end (S1 Appendix). The Pearson correlation between Cufflinks and Miso quantification in our data sets reflects this (representative scatter plot S3E Fig). Because two of the data sets that we have considered are sequenced with single-end reads and the resultant suboptimal modelling (S3F Fig), we decided to omit Cufflinks quantifications and regression models from further analyses.

Differential gene expression

Gene counts are extracted using featureCounts [31] from alignment files. Subsequently, differential gene expression is estimated using DESeq2 [32].

Differential splicing

Differential splicing is estimated using Miso [27].

Quantification of DNA methylation patterns

Processed BS-Seq data sets described in Table 1 are downloaded and DNA methylation ratio (mC basecalls to the count of all reads C as m ^ C = m C C) at individual nucleotide for CG and CA dinucleotides is accumulated as BigWig [33] tracks. DNA methylation is calculated in two metrics: 1) Region-specific mean methylation defined as m ^ C = ∑ m C N C where N_C is the number of C’s within the region 2) Region-specific methylation density defined as M ^ C = N C · m ^ C L where L is the length of the region. These two metrics are calculated for all the genes such that every gene has separate quantifications for CA and CG dinucleotides. For each gene, several specific regions are considered 1) ASE (alternative spliced exon) 2) combination of all exons containing the ASE 3) combination of all exons without the ASE 4) combination of all introns 5) combination of all introns without the ASE. A visual representation of the specific regions is described in S1 Appendix.

Quantification of ChIP-seq signal

Processed ChIP-seq data sets described in Table 1 are downloaded and l o g 2 C h I P I n p u t signal is calculated using bigWigAverageOverBed [33] for regions around the alternative exon. A visual representation of the specific regions is described in S1 Appendix.

Regressing splicing ratios with sequence features

First, we use an inverse probit transformation of our splicing ratios Ψ i.e, y = Φ⁻¹(Ψ) to map splicing ratios from the interval [0, 1] to the whole of the real numbers. Subsequently, we model the transformed splicing ratio y as a function of a set of j covariates (sequence features) X ∈ R j. We split the data, fit the model and compute the coefficient of determination (FVE) which is analogous to R², 5 consecutive times (with different splits each time) to avoid overfitting. The FVE calculated is the amount of variation in the test set that can be explained by the covariates (sequence features). F V E ( y , y ^ ) = 1 - ∑ i = 0 n ( y i - y ^ i ) 2 ∑ i = 0 n ( y i - y ¯ ) 2 where y ¯ = 1 n ∑ i = 0 n y i. This is applied for both regression approaches.

Linear regression

We model the y (transformed splicing ratio) as a linear function of X (sequence features). To regularize the weight and noise parameters in the regression problem, we follow the regularization approach [34, 35] which is implemented in the scikit-learn [36] package. Information about arguments and parameters are detailed in S1 Appendix.

Random forest regression

We generate random decision trees that are fitted on sets of training data from X (sequence features) and y (transformed splicing ratio) to form a meta estimator known as random forests [37], implemented in scikit-learn [36] package. Using averaging from a number of trees, we predict the transformed splicing ratio y on the test data.

Regressing differential splicing ratios with differential DNA methylation and sequence features

As discussed in the previous section, we use the same approach for the modelling. In this instance, we add more covariates in the form of DNA methylation features. For the mean methylation metric, we use an inverse probit transformation and take the difference.

We calculate the differential splicing ratio i.e, Δy = ΔΦ⁻¹(Ψ)

Supporting information

S1 Table [pdf]
Pearson () correlation between changes in splicing ratios and changes in DNA methylation for all data sets.

S1 Fig [tif]
MeCP2 expression across data sets considered for analysis (refer ) as counts per million quantified from RNA-seq data.

S2 Fig [a]
Number of differential gene expression and splicing events across data sets.

S3 Fig [a]
Modelling transformed splicing ratios = Φ(Ψ), estimated using different quantification methods, against sequence features .

S4 Fig [a]
DNA methylation levels in neurons expressing multiple levels of MeCP2.

S5 Fig [a]

S1 Appendix [pdf]
Detailed information about software, software versions and command line arguments used for analysis.