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

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.

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.

Table 1

Sequencing data sets considered for multi-variate quantitative analysis.

Accession	Source	Total RNA-seqCoverage	Reps	Type
GSE125660 [9]	MeCP2 WT Neurons	218.2 × 106 pairs	4	BS-seq, RNA-seq, ChIP-seq
	MeCP2 KO Neurons	214.2 × 106 pairs	4	RNA-seq, ChIP-seq
	MeCP2 OE4x Neurons	217.5 × 106 pairs	4
	MeCP2 OE11x Neurons	220.4 × 106 pairs	4
GSE128186 [17]	WT brain tissue	713.9 × 106 reads	10	BS-seq, RNA-seq
	MeCP2 KO brain tissue	796.5 × 106 reads	10
GSE64910 [18]	WT mESCs	49.1 × 106 reads	2	BS-seq, RNA-seq
	DNMT-TKO mESCs	50.3 × 106 reads	2
GSE67867 [19]	WT mESCs	357.3 × 106 reads	3	RNA-seq
	DNMT-TKO mESCs	364.8 × 106 reads	3
GSE103214 [20]	Pv Neurons (1 week)	37.7 × 106 reads	2	BS-seq, RNA-seq
	Pv Neurons (3 weeks)	27.6 × 106 reads	1
	Pv Neurons (8 weeks)	97.5 × 106 reads	2
	Vip Neurons (1 week)	17.9 × 106 reads	1
	Vip Neurons (3 weeks)	57.9 × 106 reads	2
	Vip Neurons (8 weeks)	45.3 × 106 reads	2
GSE48307 [23]GSE85517 [14]	Promyelocytes	63.5 × 106 reads	1	BS-seq, RNA-seq, ChIP-seq
	Granulocytes	59.7 × 106 reads	1

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 R2. 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.

Fig 1

Fraction of variance explained by regressing sequence features against splicing ratios.

(A) Cultured human neurons expressing multiple levels of MeCP2 [9](B) Wildtype (WT) and MeCP2 Knockout (KO) brain tissue [17](C) Pv and Vip neurons at 1 week, 3 weeks and 8 weeks [20] (D) DNMT-TKO and WT mESCs [18](E) DNMT-TKO and WT mESCs [19] (F) Promyelocytes and Granulocytes [14, 23]. See Table 1 for detailed information about data sets.

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.

Fig 2

Explanatory power of varying MeCP2 levels (in vitro and in vivo).

(A) Fraction of variance explained by regressing splicing variation against sequence features combined with MeCP2 binding in cultured human neurons [9]. (B) Fraction of variance explained by regressing splicing variation against sequence features combined with changes in DNA methylation in WT and MeCP2 KO mouse brain [17].

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.

Fig 3

Fraction of variance explained by regressing splicing variation against sequence features along with DNA methylation features.

(A) Pv and Vip neurons at 1 week, 3 weeks and 8 weeks [20](B) DNMT-TKO and WT mESCs [19](C) DNMT-TKO and WT mESCs [18](D) Promyelocytes and Granulocytes [14, 23].

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 ) 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 where N is the number of C’s within the region 2) Region-specific methylation density defined as 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 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 = Φ−1(Ψ) 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) . We split the data, fit the model and compute the coefficient of determination (FVE) which is analogous to R2, 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). where . 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. where t = treatment/time − point and c = control/time − point. For the methylation density metric, we calculate the difference directly.

We calculate the differential splicing ratio i.e, Δy = ΔΦ−1(Ψ) and model it to a set of covariates which now consists of sequence features and DNA methylation feautures.

Pearson (r) correlation between changes in splicing ratios and changes in DNA methylation for all data sets.

(PDF)

Click here for additional data file.

MeCP2 expression across data sets considered for analysis (refer Table 1) as counts per million quantified from RNA-seq data.

(TIF)

Click here for additional data file.

Number of differential gene expression and splicing events across data sets.

(A) Number of differentially regulated genes and (B) Number of differential splicing events in cultured neurons expressing multiple levels of MeCP2 [9] (C) Number of differentially regulated genes and (D) Number of differential splicing events in MeCP2 KO mouse brains compared to MeCP2 WT [17] (E) Number of differentially regulated genes and (F) Number of differential splicing events in developing mouse neurons [20] (G) Number of differentially regulated genes and (H) Number of differential splicing events in DNMT-TKO mESCs [19] (I) Number of differentially regulated genes and (J) Number of differential splicing events in DNMT-TKO mESCs [18].

(TIF)

Click here for additional data file.

Modelling transformed splicing ratios y = Φ−1(Ψ), estimated using different quantification methods, against sequence features X.

(A) Representative example of correlation between quantification methods BRIE and Miso. (B) Representative scatter plot showing prediction of splicing ratios using linear regression model learned from BRIE’s quantification and sequence features. (C) Representative scatter plot showing prediction of splicing ratios using linear regression model learned from Miso’s quantification and sequence features. (D) Comparison of splicing ratio estimation of known Ψ from simulated data at different coverage. (E) Representative example of correlation between quantification methods Cufflinks and Miso. (F) Representative scatter plot showing prediction of splicing ratios using linear regression model learned from Cufflinks’ quantification and sequence features.

(TIF)

Click here for additional data file.

DNA methylation levels in neurons expressing multiple levels of MeCP2.

(A) Referencing results from our previous work [9]. High Performance Liquid Chromatography (HPLC) quantification of methylated cytosines in neurons expressing wild-type, 4 times and 11 times MeCP2. (B) CG and (C) CA methylation across introns quantified from bisulfite sequencing of WT and MeCP2 KO mouse brains (GSE128172 [17]).

(TIF)

Click here for additional data file.

Fraction of variance explained by regressing DNA methylation features against differential splicing ratios in (A) developing mouse neurons [20](B) Promyelocytes and Granulocytes [14, 23].

(TIF)

Click here for additional data file.

Detailed information about software, software versions and command line arguments used for analysis.

(PDF)

Click here for additional data file.

13 Jul 2020

Dear Dr Sanguinetti,

Thank you very much for submitting your Research Article entitled 'Quantitative analysis questions the role of MeCP2 in alternative splicing' to PLOS Genetics. Your manuscript was fully evaluated at the editorial level and by independent peer reviewers. The reviewers appreciated the attention to an important problem, but raised some concerns about the current manuscript. Based on the reviews, we will not be able to accept this version of the manuscript, but we would be willing to review again a much-revised version. We cannot, of course, promise publication at that time.

Should you decide to revise the manuscript for further consideration here, your revisions should address the specific points made by each reviewer. We will also require a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript.

If you decide to revise the manuscript for further consideration at PLOS Genetics, please aim to resubmit within the next 60 days, unless it will take extra time to address the concerns of the reviewers, in which case we would appreciate an expected resubmission date by email to plosgenetics@plos.org.

If present, accompanying reviewer attachments are included with this email; please notify the journal office if any appear to be missing. They will also be available for download from the link below. You can use this link to log into the system when you are ready to submit a revised version, having first consulted our Submission Checklist.

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see our guidelines.

Please be aware that our data availability policy requires that all numerical data underlying graphs or summary statistics are included with the submission, and you will need to provide this upon resubmission if not already present. In addition, we do not permit the inclusion of phrases such as "data not shown" or "unpublished results" in manuscripts. All points should be backed up by data provided with the submission.

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool.  PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

PLOS has incorporated Similarity Check, powered by iThenticate, into its journal-wide submission system in order to screen submitted content for originality before publication. Each PLOS journal undertakes screening on a proportion of submitted articles. You will be contacted if needed following the screening process.

To resubmit, use the link below and 'Revise Submission' in the 'Submissions Needing Revision' folder.

[LINK]

We are sorry that we cannot be more positive about your manuscript at this stage. Please do not hesitate to contact us if you have any concerns or questions.

Yours sincerely,

Dirk Schübeler

Associate Editor

PLOS Genetics

John Greally

Section Editor: Epigenetics

PLOS Genetics

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: In this study, Chhatbar and colleagues investigate the role of MeCP2 in regulating alternative splicing. MeCP2 is a methyl-DNA binding protein, and over the last decades multiple functions have ben assigned to this protein, including regulation of alternative splicing. The implication of MeCP2 in splicing regulation was mainly based on correlative analysis. In this study, the authors utilise new datasets to identify changes in alternative splicing introduced by MeCP2 and DNA methylation. This is an important contribution towards understanding MeCP2 biology and its contribution to neuronal function. However, there are a few points that are not completely clear.

First, the authors investigate the contribution of MeCP2 to splicing through analysis of transcriptomic data from cells with different levels of MECP2 but constant total methylation data. They find that MeCP2 primarily functions as a regulator of gene expression and has minimal effect on splicing. In these comparisons, the authors mention that global DNA methylation is constant, however there is no analysis provided that supports that conclusion. Differential methylation could be restricted to a small number of regions and the authors need to exclude that the alternative splicing events (even if low in number) are not occurring at these differentially methylated sites. The authors utilise DNMT-TKO ES cells to address how much differential methylation would influence splicing, however there could be profound differences between ES cells and neurons in DNA-methylation-mediated splicing regulation, and this needs to be addressed.

Next they apply machine learning to identify the contribution of DNA methylation and MECP2 to alternative splicing. They first show that they retrieve a good prediction based on DNA sequence features, which serves as positive control. Finally they show that neither MeCP2 or DNA methylation allow prediction of splicing events, globally. This assumes that, if methylation or MeCP2 would contribute to splicing, this effect should be visible for all transcripts or exon/intron boundaries throughout the genome. This also assumes that DNA methylation and binding of MeCP2 is constant across the entire genome and over all exon/intron boundaries - which is not the case. Therefore, a context-dependent role can not be excluded where MeCP2 activity or DNA methylation would influence splicing only if certain criteria are met. E.g. only a at few splicing sites where high MeCP2 binding or high DNA methylation levels are present to influence splicing. I am not sure if the model applied by the authors fully takes these variations in consideration.

Reviewer #2: MeCP2 is an abundant meCpA and meCpG-binding protein with essential roles in neurons. In humans, its mutations cause Rett syndrome. Determining the molecular basis of the disease has been challenging. A role in transcription repression now seems established, but this does not rule out other functions. In particular, as DNA methylation is linked to splicing control, the possibility that MeCP2 regulates splicing has been raised in several publications. Most recently, Wong et al showed in a 2017 Nat Comm paper (PMID 28480880) that MeCP2 inhibits intron retention.

The current manuscript re-examines the role of MeCP2 in splicing using an advanced bioinformatic approach. More specifically, the authors use methods that allow them to tease apart the contribution of DNA methylation per se, from the contribution of MeCP2 itself. Their conclusion is that, once the confounders are removed, there is little evidence to suggest that MeCP2 actually regulates splicing extensively.

The paper is well written and argued. The figures are clear and support the conclusion. I liked the discussion, which is lucid and insightful.

I cannot judge the adequacy of the machine learning approaches used, and I have no technical comments as I am not a bioinformatician.

My main question is biological: I understand using the TKO cells as they are devoid of DNA methylation, but aren’t they also quite low for MeCP2 expression? How can we rule out that there is no effect on splicing because MeCP2 is absent or low? If the authors could show the FPKMs for MeCP2 in ES cells, in comparison to the cells used in the Nat Comm paper (promyelocytes/granulocytes), that would help interpretation.

Reviewer #3: In the current manuscript Chhatbar et al test the effect of a DNA methylation-binding protein, MeCP2, on alternative splicing. The manuscript is based on re-analyses of published data in diverse model systems: ESCs, developing neurons, and hematopoietic differentiation. The authors conclude that MeCP2 has a minimal effect on splicing, and advocate for caution when drawing conclusions that arise from observational high throughput studies. Overall, the manuscript is well-written, easy to follow and the statistical approaches seem appropriate. Nevertheless, there are caveats in this study that the authors need to address before this manuscript is made available for publication.

In the introduction, the authors claim: “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.”

However, many of the datasets used in the study are not “high quality”. For example, some of the ESC data (both wt and KO) from the Ast lab are sequenced at a depth of 40M reads (SE, 50bp) which would be insufficient or barely sufficient for proper AS data analyses. Similarly, the brain neuron RNA-seq data is sequenced at 30M or less (50bp SE). Vip neuron data are sequenced at 20M, which is fairly low and insufficient for any meaningful AS analysis. The authors should provide sequencing depth for each sample (ideally in millions of mapped reads) in Table 1.

It should also be explained better why the selected samples were chosen and how their suitability for AS analyses was determined. As sequencing depth plays a major role in AS detection, a few sentences on the role and effects of depth in the current analyses would be most welcome.

The y axes in figures S1B and S1D need to be adjusted. I don’t think it is necessary to use the same scale for differential gene expression and differential splicing. I understand that the authors want to show that the number of differential splicing events is minimal when compared to the number of differentially expressed genes, however it would be much better to represent the splicing data on a separate scale (ie. 1-100 or something like that).

It would be beneficial if the authors could produce a figure with MeCP2 expression levels across different tissues / samples used in the current manuscript.

Finally, as this study only discusses the global impact of MeCP2 on AS and does not try to confirm or refute potential locus specific AS events impacted by MeCP2 binding, the title needs to be adjusted. Something like: “Quantitative analysis questions the role of MeCP2 as a genome-wide (global?) regulator of alternative splicing” would be a better fit.

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Genetics data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: None

Reviewer #2: Yes

Reviewer #3: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: No

11 Aug 2020

Submitted filename: Response_to_reviewers.pdf

Click here for additional data file.

28 Aug 2020

Dear Dr Sanguinetti,

We are pleased to inform you that your manuscript entitled "Quantitative analysis questions the role of MeCP2 as a global regulator of alternative splicing" has been editorially accepted for publication in PLOS Genetics. Congratulations!

Before your submission can be formally accepted and sent to production you will need to complete our formatting changes, which you will receive in a follow up email. Please be aware that it may take several days for you to receive this email; during this time no action is required by you. Please note: the accept date on your published article will reflect the date of this provisional accept, but your manuscript will not be scheduled for publication until the required changes have been made.

Once your paper is formally accepted, an uncorrected proof of your manuscript will be published online ahead of the final version, unless you’ve already opted out via the online submission form. If, for any reason, you do not want an earlier version of your manuscript published online or are unsure if you have already indicated as such, please let the journal staff know immediately at plosgenetics@plos.org.

In the meantime, please log into Editorial Manager at https://www.editorialmanager.com/pgenetics/, click the "Update My Information" link at the top of the page, and update your user information to ensure an efficient production and billing process. Note that PLOS requires an ORCID iD for all corresponding authors. Therefore, please ensure that you have an ORCID iD and that it is validated in Editorial Manager. To do this, go to ‘Update my Information’ (in the upper left-hand corner of the main menu), and click on the Fetch/Validate link next to the ORCID field.  This will take you to the ORCID site and allow you to create a new iD or authenticate a pre-existing iD in Editorial Manager.

If you have a press-related query, or would like to know about one way to make your underlying data available (as you will be aware, this is required for publication), please see the end of this email. If your institution or institutions have a press office, please notify them about your upcoming article at this point, to enable them to help maximise its impact. Inform journal staff as soon as possible if you are preparing a press release for your article and need a publication date.

Thank you again for supporting open-access publishing; we are looking forward to publishing your work in PLOS Genetics!

Yours sincerely,

Dirk Schübeler

Associate Editor

PLOS Genetics

John Greally

Section Editor: Epigenetics

PLOS Genetics

www.plosgenetics.org

Twitter: @PLOSGenetics

----------------------------------------------------

Comments from the reviewers (if applicable):

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The authors have addressed all my concerns. The inclusion of new datasets now better support their initial conclusions. I furthermore appreciate the change of the title.

Reviewer #2: My concerns on the first version were rather minimal and have been addressed in the revision

Reviewer #3: The authors have successfully addressed all of the reviewers' concerns. I have no further queries.

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Genetics data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: No

----------------------------------------------------

Data Deposition

If you have submitted a Research Article or Front Matter that has associated data that are not suitable for deposition in a subject-specific public repository (such as GenBank or ArrayExpress), one way to make that data available is to deposit it in the Dryad Digital Repository. As you may recall, we ask all authors to agree to make data available; this is one way to achieve that. A full list of recommended repositories can be found on our website.

The following link will take you to the Dryad record for your article, so you won't have to re‐enter its bibliographic information, and can upload your files directly:

http://datadryad.org/submit?journalID=pgenetics&manu=PGENETICS-D-20-00858R1

More information about depositing data in Dryad is available at http://www.datadryad.org/depositing. If you experience any difficulties in submitting your data, please contact help@datadryad.org for support.

Additionally, please be aware that our data availability policy requires that all numerical data underlying display items are included with the submission, and you will need to provide this before we can formally accept your manuscript, if not already present.

----------------------------------------------------

Press Queries

If you or your institution will be preparing press materials for this manuscript, or if you need to know your paper's publication date for media purposes, please inform the journal staff as soon as possible so that your submission can be scheduled accordingly. Your manuscript will remain under a strict press embargo until the publication date and time. This means an early version of your manuscript will not be published ahead of your final version. PLOS Genetics may also choose to issue a press release for your article. If there's anything the journal should know or you'd like more information, please get in touch via plosgenetics@plos.org.

2 Oct 2020

PGENETICS-D-20-00858R1

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

Dear Dr Sanguinetti,

We are pleased to inform you that your manuscript entitled "Quantitative analysis questions the role of MeCP2 as a global regulator of alternative splicing" has been formally accepted for publication in PLOS Genetics! Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out or your manuscript is a front-matter piece, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Genetics and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Matt Lyles

PLOS Genetics

On behalf of:

The PLOS Genetics Team

Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom

plosgenetics@plos.org | +44 (0) 1223-442823

plosgenetics.org | Twitter: @PLOSGenetics