MBNL1 and RBFOX1 co-regulate alternative splicing events transcriptome-wide through a conserved buffering mechanism

Alternative splicing (AS) is controlled by cis-regulatory elements recognized by networks of trans-acting factors. Here we investigate modes and mechanisms of AS co-regulation by MBNL1 and RBFOX1, two RNA binding proteins (RBPs) critical for developmental AS transitions. We generated two cell models that express each RBP under separate inducible promoters. Transcriptome-wide categorization of the impacts of RBFOX1 expression on MBNL1 splicing revealed a common co-regulatory mode through which RBFOX1 buffers MBNL1 dose-dependent splicing regulation by reducing the total range of exon inclusion or exclusion. Minigene mutational analysis and in vitro binding experiments suggest that this buffering mechanism occurs through a shared cis-regulatory element previously unidentified as critical for MBNL1-dependent activity. Overall, our studies define a conserved co-regulatory mechanism through which RBFOX1 and MBNL1 can fine-tune and provide redundancy for AS outcomes. These studies indicate overlapping use of RNA motifs with potential implications for when activity of RBPs is disrupted.

HEK-293 cell line transiently transfected with a HA-RBFOX1 expression vector. Expression of HA-MBNL1 93 alone via doxycycline treatment following transfection of the minigene reporter led to significantly increased 94 INSR exon 11 inclusion as quantified by percent spliced in (Ψ) (i.e. percent exon inclusion) (Fig 1b). A similar, 95 albeit reduced level, of inclusion was observed following expression of HA-RBFOX1 in the absence of MBNL1 96 induction (Fig 1b). While co-expression of both splicing factors further enhanced exon inclusion above that of 97 each RBP individually, the overall ΔΨ was not additive as maximal exon inclusion (i.e. Ψ = 1.0) was not 98 achieved (Fig 1b). We observed a similar pattern of splicing activity for the MBNL1-dependent Nfix exon 8 99 minigene reporter 36 where expression of either RBP alone induced exon exclusion while co-expression led to a 00 small, but not additive combinatorial effect (Fig S1). The lack of an additive impact on exon inclusion under co-01 expression of RBFOX1 and MBNL1 suggests that these two RBPs likely do not function independently but 02 rather in a cooperative manner to modulate exon inclusion. 03 To further discern the effects of RBFOX1 and MBNL1 co-expression on AS of the INSR minigene, we 04 measured changes of exon 11 inclusion across a range of MBNL1 protein expression in the presence or 05 absence of HA-RBFOX1. Following transfection of a HA-RBFOX1 expression vector, MBNL1 expression was 06 induced by titration of doxycycline to produce a range of MBNL1 concentration and splicing activity (Ψ) plotted 07 with minimal impacts at high MBNL1 concentrations. This conclusion is supported by the fact that almost half 23 the amount of MBNL1 is required to affect an equivalent amount of splicing regulation when HA-RBFOX1 is co-24 expressed (log (EC50) = 0.22 ± 0.03 (+ RBFOX1) vs. 0.13 ± 0.03 (-RBFOX1)). Overall, MBNL1 and RBFOX1 25 appear to operate through a shared mechanism of AS regulation within the context of the INSR exon 11 26 minigene splicing reporter. 27 Figure 1: RBFOX1 buffers MBNL1 dose-dependent splicing regulation of INSR exon 11 minigene through binding to a shared cis-regulatory element. (a) Distribution of MBNL1 (YGCY) and RBFOX1 (UGCAUG) binding motifs in wild-type, MUT1, and MUT2 INSR exon 11 minigenes. (b) Bar-plot representation of cell-based splicing assay with wild-type INSR exon 11 minigene in HA-MBNL1 tetracycline-inducible HEK-293 cells 18 transfected with either empty vector or HA-RBFOX1 individually or in combination. HA-MBNL1 expression was induced via doxycycline (10 ng/mL). Data represented as mean percent exon inclusion (Ψ) ± SEM, n=3. (* p < 0.05, **** p < 0.0001, one-way ANOVA). (c-e) MBNL1 dose-response curves for wild-type (WT), MUT1, and MUT2 INSR minigenes in presence or absence of HA-RBFOX1 expression. Ψ was measured at each concentration and plotted as a function of log (dox) and fit to a four-parameter dose-curve. Data represented as mean ± SEM, n = 3. Representative native gels and quantitative parameters derived from dose-response curves are supplemented ( Fig S3, Table S1). (f) Representative competitive EMSA showing binding of both purified MBNL1 and RBFOX-RRM to a model RNA containing the INSR intron 11 putative UGCAUG RNA motif and surrounding sequence. (g) Quantification of percent RNA bound by either RBFOX1 or MBNL1 upon challenge with increasing concentrations of MBNL1 (n = 3). tested the function of the predicted UGCAUG RBFOX binding site within the downstream intron 11. This cis-31 regulatory motif was mutated to UCGAUG to generate the INSR MUT1 minigene (Fig 1a). Inversion of this GC 32 dinucleotide has been previously shown to abrogate RBFOX1 binding and prevent AS regulation [37][38][39] . As 33 predicted, co-expression of HA-RBFOX1 had no effect on splicing of the INSR MUT1 minigene (Fig 1d, red  34 curve). Surprisingly, this mutation significantly blunted MBNL1 dose-dependent splicing regulation compared to 35 WT INSR (Fig 1d, compare grey and blue curves). While this result is consistent with RBFOX1 regulating exon 36 inclusion through this specific cis-regulatory element, the substantial impact on MBNL1 splicing regulation was 37 unexpected given that (1) many MBNL1 YGCY binding motifs are present within intron 11 and have previously 38 been shown facilitate MBNL1-dependent exon 11 inclusion [32][33][34][35] and (2) (Fig 1a). This mutation disrupts the RBFOX binding site as in INSR MUT1 but generates a new, YGCY 45 MBNL binding site in the same position as in the WT reporter (CGCU). Alteration of this motif nearly fully restored 46 MBNL1's ability to regulate splicing and, as predicted, co-expression of HA-RBFOX1 had no impact on AS (Fig  47   1e). Characterization of the MBNL1 dose-response for all three INSR minigene variants in the presence or 48 absence HA-RBFOX1 expression are consistent with a model in which MBNL1 and RBFOX1 bind an overlapping 49 site to modulate exon 11 inclusion. Furthermore, these experiments support a buffering co-regulatory mechanism 50 whereby at low concentrations of MBNL1, RBFOX1 binds to the target site to promote exon inclusion. As the 51 cellular concentrations of MBNL1 protein increases, MBNL1 may outcompete RBFOX1 to bind the same, 52 overlapping site such that the maximum Ψ is achieved. 53 To evaluate this proposed model, we conducted a competitive electrophoretic mobility shift assay (EMSA) 54 using a section of INSR intron 11 containing the putative RBFOX UGCAUG cis-regulatory element; no YGCY 55 Figure 2: MBNL1 dose-dependent splicing is broadly conserved in novel ddHEK and ddMEF double dosing cell lines in which MBNL1 and RBFOX1 expression levels are independently and precisely controlled. (a) Model depicting inducible MBNL1 and RBFOX1 constructs in HEK-293 cells (ddHEKs) and Mbnl1 -/-/Mbnl2 -/-MEFs (ddMEFs). HA-MBNL1 (ddHEK) or GFP-MBNL1 (ddMEF) expression is controlled via a tetracycline-inducible system. Doxycycline (dox) treatment induces expression of the full-length MBNL1 mRNA by de-repressing inhibition of transcription by binding to the tet repressor (tetR). Independently, expression of mOrange-tagged RBFOX1 is controlled via an ecdysone-inducible system. (b) Representative immunoblot showing gradient of MBNL1 expression produced via titration of doxycycline in the presence or absence of RBFOX1 expression induced via ponasterone A treatment. (c) Quantification of MBNL1 immunoblot normalized to GAPDH loading control in the presence and absence of RBFOX1 expression plotted as a function of log (dox). Data represented as mean ± SEM, n = 4. (d-f) MBNL1 dose-response curves for orthologous cassette exons NUMA1/Numa1, EXOC1/Exoc1, and INSR/Insr in ddHEK and ddMEF cell lines. Data represented as mean Ψ ± SEM (n = 3 -5) plotted against log [MBNL1] levels as quantified by immunoblot (Fig 2c). Quantitative parameters derived from MBNL1 dose-response are available in Fig  S4b, d. Distribution of YGCY motifs within 250-300bp upstream and downstream of each orthologous cassette exon is also displayed. In the case of INSR/Insr, the distribution of the RBFOX1 UGCAUG binding motif is also annotated. independent of RBFOX1 co-expression. While previous work has identified and characterized MBNL1 dose-99 dependent splicing regulation 18,19 , little work has been performed to date to assess if concentration-dependent 00 splicing outcomes are conserved across species. We selected several endogenously expressed cassette exons 01 previously shown to be impacted by MBNL1 expression 41-43 , including INSR/Insr exon 11, and quantified exon 02 inclusion levels across the MBNL1 concentration gradient in both systems. Ψ values were then plotted against 03 log [MBNL1] for both MBNL1-regulated inclusion and exclusion events (Fig 2d-f and Fig S4). Consistent with 04 previous studies, estimated EC50 values for cassette exon events assayed within the ddHEKs indicated that a 05 range of MBNL1 concentrations are required to effectively regulate exon inclusion with variable magnitudes of 06 cooperativity as measured by slope of the dose-response 18 (Fig S4a-b). These same general trends were also 07 observed within the ddMEFs (Fig S4c-d). In fact, several orthologous cassette exons expressed in both systems 08 demonstrated similar regulation as quantified by EC50 and slope upon MBNL1 titration, possibly as a 09 consequence of similar YGCY motif distribution (Fig 2d-e). This assessment of a limited number of endogenous 10 splicing events indicated that (1) MBNL1 dose-dependent splicing regulation is preserved across our human and 11 murine cell models, and (2) a wide range of dose-response curves described by a range of quantitative 12 parameters can be observed across both systems. 13 Interestingly, we observed substantial differences in the MBNL1 dose-response for INSR/Insr exon 11 14 between our murine and human dosing models. MBNL1 strongly promoted exon inclusion in ddHEKs (ΔΨ = 15 0.58) (Fig 2f), consistent with our minigene results (Fig 1c). However, the opposite effect was observed in the 16 ddMEFs where MBNL1 expression mildly stimulated exon 11 exclusion (ΔΨ = -0.10) (Fig 2f). Sequence analysis 17 of the mouse Insr exon 11 downstream intron uncovered a low number of putative MBNL1-binding motifs within 18 300 nucleotides of exon 11 (Fig 2f), consistent with the observed mild effect of MBNL1 expression and 19 established patterns of MBNL positional dependent AS regulation 17 . Interestingly, while the overall number of 20 YGCY motifs was reduced within the murine upstream intron 11 compared to the corresponding human intron, 21 the UGCAUG RBFOX1 binding motif previously identified as critical for MBNL1-dependent exon inclusion activity 22 (Fig 1f) was preserved. Despite the presence of this cis-regulatory element, the contrasting effect of MBNL1 23 expression within the ddMEF system underscores the complexity of AS captured by these dosing models. 24 Furthermore, these studies highlight the capacity of these companion cell lines to also define differences in acting factors. Selecting to exclusively focus on skipped exon (SE) splicing events, a list of common events between these 39 three pairwise comparisons was compiled to capture a subset of events in which MBNL1-dependent splicing is 40 affected by RBFOX1 expression. We then subjected this compiled subset to further filtering to identify events 41 where Ψ was significantly altered by MBNL1 expression (I IΔΨII -II  0.1, FDR < 0.05). This analysis yielded a 42 total of 866 and 287 events in ddHEKs and ddMEFs, respectively, which we deemed MBNL1-regulated. 43 . (e-f) AS events in category E do not show depletion of GCAUG motifs as when compared to all other co-regulatory groups (e). YGCAUG motifs are enriched in co-regulatory mode F compared to group E in ddHEKs, but not ddMEFs (f). Data represented as frequency (total number of motifs/total length of intronic sequence for the 500bp intronic region upstream and downstream of cassette exon) ± SEM, * p < 0.05, n.s. = not significant, unpaired t-test. set of mock MBNL1 dose-response curves for which RBFOX1 positively, negatively, or does not significantly 48 alter Ψ at low, high, or both endpoints (Fig 3b). These mock curves are binned into 9 inclusion/exclusion groups 49 (A -I) whereby equivalent impacts on the MBNL1 dose-response for inclusion or exclusion events are placed in 50 matching groups (Fig 3b). To categorize the 866 ddHEK and 287 ddMEF SE events into these 9 groups in an 51 unbiased manner, we devised a binary methodology in which splicing changes that exceeded a significance 52 threshold (ΔΨIII-I / ΔΨIV-II  0.1 or ΔΨIII-I / ΔΨIV-II  -0.1) were assigned "1" or "-1", respectively, and non-significant 53 changes 0.1 > ΔΨIII-I / ΔΨIV-II > -0.1) were assigned "0" (Fig 3b, brackets). Each event was assigned to a group 54 based upon the change induced by RBFOX expression at low and high MBNL1 expression. For example, if 55 RBFOX expression resulted in a ΔΨ > -0.1 (scored as -1) at low MBNL expression and a ΔΨ > 0.1 (scored as 56 1) at high MBNL expression, the event was binned as [-1,1] or group A (Fig 3b). 57 After all SE events were appropriately binned (Fig 3c) we assessed the general distribution across each 58 co-regulatory group. AS events were in general disproportionately binned across the various postulated groups 59 with a large portion of events in a co-regulatory group (all excluding E). Interestingly in both cell lines, significantly 60 more events in which RBFOX1 exclusively shifted Ψ in the same direction as MBNL1 (groups B, C, F) were 61 identified as compared to those in which RBFOX1 displayed negative co-regulation (groups D, G, H), with twice 62 as many positive co-regulated events observed in the ddHEKs (68.1% vs. 31.9%, respectively). This observation 63 is consistent with established patterns of developmental AS regulation whereby both RBPs have been shown to 64 promote fetal to adult AS transitions 13,29 . When looking at each group individually we found the highest number 65 of events (35% ddHEKs; 25% ddMEFs) bin into group E representing those AS events not impacted by RBFOX1 66 expression (Fig 3c, d). Of co-regulatory groups, the highest number of events binned into group F (164 in ddHEKs 67 and 77 in ddMEFs), followed by co-regulatory groups H and I (Fig 3c, d). Interestingly, group F most closely 68 resembles the dose-response associated with the buffering effects of RBFOX1 observed with the INSR minigene 69 (Fig 1). Importantly, this suggests that a comparatively large portion of co-regulated events may utilize a similar 70 shared motif mechanism to facilitate the observed patterns of RBFOX1 and MBNL1 AS co-regulation. 71 was not as pronounced for group F, especially in the ddHEK system. Group F generally retained high numbers 76 of events (133 events for ddHEKs and 42 for ddMEFs) at the highest threshold (ΔΨ = 0.30, Fig S5c) indicating 77 that most AS events that bin into co-regulatory group F are a consequence of a substantial impact on the MBNL1 78 dose-response by RBFOX1 expression at low MBNL1 concentrations. Conversely, this observation is consistent 79 with a more subtle effect of RBFOX1 co-expression on other types of co-regulatory modes. 80 Finally, we assessed conservation based on translated sequence homology for all individual grouped 81 events in both systems. We identified 71 and 74 splicing events of the 866 ddHEK and 287 ddMEF events, 82 respectively, in which the alternatively regulated exonic translated sequence had  70% identity compared to the 83 matching one-to-one orthologous gene. Of those orthologous events, 16 were binned into matching groups, 9 of 84 which were specific to group F (Table S3). While the total number of orthologous events identified was small, 85 likely due to the distinct species and tissue origin of both cell systems, the large proportion of events with 86 matching patterns of AS co-regulation that binned into group F further emphasizes the prominence of this co-87 regulatory mechanism. Additionally, these 9 events are good candidates for future targeted mechanistic studies 88 to evaluate the relative importance of specific of cis-regulatory motifs that may contribute to the shared patterns 89 of AS outcomes across systems. 90 Overall, by using postulated mock curves to bin AS events into 9 distinct groups we showed that there is 91 an unequal distribution of events between groups with the majority binned into groups demonstrating co-92 regulation. Furthermore, in both cell systems the most prevalent amongst those co-regulated (group F) having a 93 similar curve architecture to the INSR minigene experiment suggesting group F highlights an important co-94 regulatory mechanism in human and mouse. 95 Presence of putative RBFOX1 RNA binding motifs does not accurately predict distribution of AS events 97 into specific co-regulatory modes 98 distribution of events within AS co-regulation groups (Fig 3c), we examined the 500bp intronic sequence region 00 flanking each skipped exon event for individual motifs. Given that events were initially filtered based on a 01 significant change in Ψ upon MBNL1 expression, we confirmed that >99.5% of events in both systems contained 02 putative YGCY regulatory motifs (data not shown). We first chose to focus on events that bin into group E, which 03 display no co-regulatory capacity. We predicted that these MBNL1-dependent events are unable to be co-04 regulated by RBFOX1 due to the absence of an RBFOX binding GCAUG motif. Interestingly, analysis revealed 05 that only 40.5% and 25% of events in ddHEKs and ddMEFs, respectively, lacked a GCAUG motif in the flanking 06 intronic regions (data not shown) suggesting a majority of events are in group E not due to the lack of GCAUG 07 motifs. Alternatively, it is possible that the lack of co-regulation is due to an overall relative depletion of GCAUG 08 motifs within the upstream and downstream intronic sequences. To test this alternative hypothesis, we compared 09 the frequency of GCAUG motifs (number of GCAUG motifs/total sequence length) in group E to all other co-10 regulatory modes. We found that events in group E did not exhibit a significant depletion of GCAUG motifs in 11 either cell system (Fig 3e). This suggests that other factors beyond the presence or absence of GCAUG motifs 12 contribute to lack of response to RBFOX1. 13 In a similar manner to group E, we next analyzed the distribution of motifs in group F, which has the 14 highest number of binned events of all other co-regulation groups. Furthermore, the dose-response curve 15 architecture of group F most closely resembles the dose response associated with the buffering effect of our 16 INSR minigene experiment (Fig 1) in which we showed MBNL1 and RBFOX1 use a shared UGCAUG motif to 17 regulate splicing. To assess the potential for other events in group F to also utilize this shared motif mechanism, 18 we scanned the flanking intronic space for the presence of YGCAUG motifs. The leading pyrimidine nucleotide 19 (Y = C or U) allows for an MBNL binding motif (YGCA) to overlap with a RBFOX binding element. Surprisingly, 20 we found that only 30.5% in ddHEKs and 44% in ddMEFs of events in group F contain a YGCAUG motif (data 21 not shown) suggesting a portion of events use a different co-regulatory mechanism. Again, we assessed the 22 alternative hypothesis of whether the frequency of YGCAUG motifs is enriched in this group compared to non 23 co-regulated events (group E). Interestingly, we found that in our ddHEK cell system YGCAUG motifs are in fact 24 significantly enriched in group F compared to E; this enrichment was not observed in the ddMEFs (Fig 3f). This 25 inconsistency is likely due to the significantly lower overall number of SE events (866 ddHEK vs 287 ddMEF, Fig  26   3c) between our cell models. 27 Overall, we show here that the absence or presence and distribution of cis-regulatory motifs alone is not 28 sufficient to accurately predict splicing outcomes of individual splicing events for MBNL1 and RBFOX1. Splicing 29 events that are or are not co-regulated by MBNL1 and RBFOX1 exhibit an association beyond that of cis-30 regulatory factors, and both co-regulatory and non-co-regulatory groups are likely dictated by other trans-acting 31 factors or unknown RNA secondary/tertiary structure. 32 33 Full dose-response curves illustrate buffering effects of induced RBFOX1 expression on MBNL1 dose-34 dependent splicing regulation 35 To further explore and validate the authenticity of the buffering effects of RBFOX1 expression and cis-36 motif distribution on MBNL1 dose-dependent splicing regulation, endogenous dose-response curves were 37 generated for several splicing events in group F across the complete, inducible MBNL1 concentration range 38 within both dosing models. Based on our experiments with the INSR minigene reporter (Fig 1) and ΔΨ thresholds 39 used to categorize co-regulated events, we hypothesized that MBNL1-dependent SE exons within the F module 40 are buffered by RBFOX1 due to an overall reduced ΔΨ (IΔΨ-RBFOX1I > IΔΨ+RBFOX1I) covered over the MBNL1 41 concentration gradient through shared use of a shared YGCAUG RNA regulatory motif. 42 ITGA6/Itga6, EXOC1/Exoc1, and NUMA1/Numa1 are examples of three orthologous SE events classified 43 within group F in both ddHEK and ddMEF cell models (Fig 4a-4c). These three events have been previously 44 reported to be regulated by both RBPs, although in an independent context 24,41-45 . Consistent with our RNAseq 45 data, RBFOX1 expression buffered the endogenous MBNL1 dose-response by positively shifting exon inclusion 46 at low MBNL1 concentrations with minimal impacts at high MBNL1 levels, leading to an overall reduction in ΔΨ 47 (Table S4). While log (EC50) values were not significantly altered, the slope of the dose-response curves in all 48 cases was increased (ITGA6/Itga6, EXOC1/Exoc1) or only minimally decreased (NUMA1/Numa1), likely as a 49 consequence of the reduced ΔΨ+RBFOX1 covered over a minimal [MBNL1] range (Table S4). Interestingly, this 50 result indicates an increase in cooperative MBNL1 regulation, potentially as a result of dual use of shared cis-51 regulatory motifs with RBFOX1. In fact, sequence analysis within the flanking regions of these events revealed 52 predicted putative YGCY and GCAUG motifs positionally organized to drive the observed patterns of exon 53 inclusion or exclusion for each RBP (Fig 4a-c). Additionally, as predicted for this co-regulatory mode, the 54 available GCAUG motifs possess a leading pyrimidine base, specifically a UGCAUG, embedded within several 55 putative YGCY motifs. This configuration is predicted, based on our earlier data, to be utilized by both RBPs to 56 facilitate AS co-regulation (Fig 4a-4c). 57 58 While the distribution of motifs and splicing outcomes was generally preserved for these events within 59 both the ddHEK and ddMEFs (Fig 4a-4c), we observed differential patterns of splicing regulation for SYNE1 in 60 the ddHEK and ddMEF cell lines. Syne1 displayed a change in the MBNL1 dose-response curve consistent with 61 Figure 4: Full splicing dose-response curves illustrate variable impacts of RBFOX1 on MBNL1 dosedependent splicing regulation. (a-d) MBNL1 dose-response curves in the presence or absence of RBFOX1 expression for ITGA6/Itga6, EXOC1/Exoc1, NUMA1/Numa1, and SYNE1/Syne1, respectively. MBNL1 expression was generated via titration of doxycycline. Additionally, cells were treated with a vehicle control (-RBFOX1) or a single concentration of ponasterone A to induce RBFOX1 expression. Percent exon inclusion (Ψ) was measured for each endogenously expressed mRNA target, plotted against log [MBNL1] levels (Fig 2c), and then fit to a four-parameter dose curve. Data represented as mean ± SEM, n = 3 -5. The assigned co-regulatory group (A-I) as derived from RNAseq (ΔΨ threshold = 0.10, Fig 3f) is listed for each event and the associated mock curvature displayed below. Quantitative parameters derived from dose-response curves are listed in Table S4. Distribution of putative YGCY, YGCAUG, and RGCAUG binding motifs within and flanking the regulated exon (up to 500 nucleotides) are displayed.
flanking intronic regions revealed that SYNE1 contained two additional UGCAUG motifs in the upstream intron 63 in the ddHEK that would be predicted to contribute to the observed buffering effect; these motifs were absent in nucleotide (R = A or G) predicted to prevent MBNL1 binding while not affecting RBFOX1. 66 Several other orthologous SE events assayed also did not uniformly sort into the same co-regulatory 67 mode (NUMB/Numb (E/F), PLOD2/Plod2 (C/F), ADD3/Add3 (F/I), and INF2/Inf2 (D/F) (Fig S6). While some 68 events contained the expected patterns of cis-motif distribution to facilitate the observed patterns of exon 69 inclusion across the MBNL1 concentration gradient, no uniform or distinct differences in the distribution of cis-70 regulatory elements adjacent to each cassette exon explained the variability in observed AS co-regulation across 71 systems ( Fig S6). Collectively, the endogenous dose-response curves for multiple orthologous SE events 72 underscores the complex variability of AS co-regulation by RBFOX1 and MBNL1 across distinct systems. These 73 differences continue to highlight our lack of understanding of the interplay between both putative cis-motifs 74 landscapes and trans-acting RBP concentrations and their implications on AS co-regulation. 75 76 Discussion: 77 RBFOX1 buffers MBNL1 dose-dependent splicing in two distinct cellular systems 78 In this study we investigated potential modes and mechanisms of AS co-regulation between two 79 splicing factors, MBNL1 and RBFOX1. Using a minigene reporter system and both human and mouse dual-80 inducible cell lines we evaluated the effects of RBFOX1 expression on MBNL1 dose-dependent splicing 81 regulation. The results demonstrated that MBNL1 and RBFOX1 co-regulate an extensive collection of AS 82 events transcriptome-wide in an interdependent manner. To effectively analyze patterns of AS co-regulation 83 we developed a grouping system to categorize AS events into 9 inclusion/exclusion groups to compare 84 patterns of AS events across cell lines or conditions (Fig 5a). Analysis from our two dual-inducible cell lines 85 revealed that the most predominant co-regulatory group (F) is one in which RBFOX1 buffers the MBNL1 dose-86 response to fine-tune exon inclusion (Fig 3c-d). Based on mutagenesis of the INSR exon 11 minigene and in 87 vitro binding experiments, we propose a novel mechanism by which the RBFOX UGCAUG cis-regulatory RNA 88 redundancy for AS outcomes, which provides flexibility in maintaining AS patterns beyond that of each RBP 94 alone. 95 Complex interplay of cis-regulatory landscape use by MBNL and RBFOX to facilitate AS outcomes 96 While MBNL1 and RBFOX1 have previously been shown to bind to the same RNAs using a common 97 motif 31 , here we present evidence that MBNL1 and RBFOX1 compete for use of the same UGCAUG cis-98 regulatory element containing a sub-optimal MBNL binding element (YGCA) (Fig 1). While in the context of the 99 INSR minigene, mutation of this site had profound effects on splicing regulation by both RBPs, predicting 00 patterns of global AS co-regulation based on presence of specific cis-regulatory elements proved difficult. In 01 fact, minimal RBFOX1 motifs (GCAUG) were not significantly depleted within non co-regulated events (group 02 E) (Fig 3e). Conversely, events that were buffered (group F) by RBFOX1 showed general enrichment for 03 Figure 5: Use of dual-dosing cell models and unbiased categorization of AS co-regulation indicates widespread buffering of MBNL1 dose-dependent splicing by RBFOX1 mediated through use of shared cisregulatory motifs. (a) RNAseq was performed on dual MBNL1 and RBFOX1 inducible cell lines at minimal and high expression levels of one or both splicing factors. Following identification of AS events regulated by MBNL1 alone, the effects of RBFOX1 expression at the endpoints of the MBNL1 dose-response curve were categorized into 9 predicted co-regulatory modes based on impacts of RBFOX 1 on endpoints of the MBNL1 dose-response curve. (b) Analysis of the impacts of RBFOX1 using the methodology described above in two cell systems revealed the prevalence of coregulated events in which RBFOX1 buffered MBNL1 concentration dependent-splicing. Experiments with the INSR minigene reporter indicate that this co-regulatory mechanism can be facilitated by shared use the canonical UGCAUG RBFOX binding element. At low concentrations of MBNL, RBFOX binds this cis-regulatory motif to drive exon inclusion. However, at increased levels of MBNL proteins, RBFOX can be displaced and the low-specificity UGCA element is bound by MBNL to maximize exon inclusion. Furthermore, the general percentage of events in group E that lacked a GCAUG or events in group F that 05 contained a YGCAUG was much lower or higher, respectively, than expected. Collectively, these results 06 suggest that global patterns of AS co-regulation by MBNL1 and RBFOX1 are likely driven by a highly 07 populated and diverse cis-regulatory landscape rather than few, high affinity regulatory motifs. Specifically, the 08 expansion of the functional regulatory landscape through use of suboptimal or intermediate affinity motifs may 09 serve to enhance binding to higher-affinity, primary motifs or may only occur at high cellular concentrations to 10 narrow AS activity. Furthermore, additional trans-acting protein factors and/or unknown but biologically relevant 11 RNA secondary or tertiary structure may be important to understanding global patterns of AS. These results 12 speak to the complex use of finite cis-regulatory space by a multitude of trans-acting factors, especially as it 13 relates to how the regulatory landscape can be altered as concentrations of such factors shift. 14 15 Implications of MBNL1 and RBFOX1 cooperatively regulating AS outcomes in development and splicing switches that are completed upon MBNL expression, a phenomenon that is mechanistically consistent 23 with our proposed buffering model, in which initial expression of RBFOX1 alters exon inclusion that is 24 maximized as MBNL1 expression increases. Furthermore, many of the events (e.g. PLOD2, ITGA6, SYNE2) 25 previously shown to facilitate developmental maturation via cooperative regulation were also identified and 26 categorized into buffering group F in our systems. Together, these results provide a potential mechanism 27 where differential expression of these two RBPs cooperate during development to interdependently regulate 28 AS transitions. 29 study highlights the potential for MBNL and RBFOX protein families to functionally compensate for the other to 31 preserve AS patterns in instances where the cellular concentrations or functionality of the other RBP is altered 32 or compromised 31 . This ability may be especially impactful in the context of disease, specifically in the context 33 of Myotonic Dystrophy type 1 (DM1) and type 2 (DM2) in which toxic expansion CUG(G) repeat RNA 34 sequesters and reduces functional MBNL protein levels [47][48][49] . The presence of RBFOX in MBNL-depleted 35 tissues may compensate for the impacts of functional MBNL loss, potentially mitigating the degree of mis-36 splicing observed in these tissues. In fact, several events mis-spliced in DM that are associated with disease 37 symptoms are also classified as part of the co-regulatory group F, including INSR and SCN5A 50-52 . In fact, 38 increased RBFOX2 expression highly correlates with reduced inferred levels of functional, non-sequestered 39 MBNL proteins based on RNAseq data from a large cohort of DM1 patient muscle 53 . Previous work has also 40 shown DM1/DM2 induced pluripotent stem cell derived cardiomyocytes have increased RBFOX1 protein 41 levels 54 . These results suggest that RBFOX expression may be upregulated in DM as a compensatory 42 mechanism to mitigate disease-specific mis-splicing. As such, changes in expression of the family of RBFOX 43 splicing factors and its impacts on MBNL-regulated splicing events may in part explain phenotypic variation 44 observed in DM patients 55 . 45 Furthermore, this functional compensation may have important consequences on DM splicing 46 biomarker selection for therapeutic trials. Given the close association between mis-splicing and DM 47 phenotypes, one of the goals of the DM field is the development of effective biomarkers that estimate changes 48 in functional levels of MBNL1 upon therapeutic intervention 18,56 . Naturally, many of these studies focus on 49 MBNL1-dependent AS events as candidate biomarkers. Candidate AS biomarkers that are buffered by 50 RBFOX1 expression such that the changes in splicing are no longer significantly altered by progressive 51 changes in free MBNL1 levels may no longer serve as adequate therapeutic biomarkers. Therefore, it is 52 important to consider RBFOX expression when developing DM biomarkers in model systems. cells were plated in twenty-four well plates at a density of 1.5 x 10 5 cells/well. Cells were transfected 24 hours 07 later at approximately 80% confluence. Plasmids were transfected into each well with 1.5 μl of TransIT-293 08 (Mirus Bio) as per the manufacturer's protocol. 250 ng of empty vector (pcDNA3.1+, mock) or the pcDNA3.1+ 09 HA-RBFOX1 expression vector were co-transfected with 250ng of a single minigene reporter (total = 500 10 ng/well). Fresh doxycycline (Sigma Aldrich) was prepared at 1 mg/mL, diluted, and added to the cells at the 11 appropriate concentrations four hours post transfection (10ng/mL). 24 hours post-transfection cells were 12 harvested for experimental analysis. 13 ddHEKs were routinely cultured as a monolayer in DMEM (Corning) supplemented with 10 % FBS 14 (Corning), 10 μg/mL blasticidin (Gibco), 150 μg/mL hygromycin (Gibco), and 2 μg/mL puromycin (Gibco) at 37 15 °C under 5 % CO2. ddMEFs were also cultured as a monolayer in DMEM (Corning) supplemented with 10% 16 FBS (Corning), 1% penicillin-streptomycin (Gibco), and 2 μg/mL puromycin (Gibco). For experimental setup, 17 cells were plated in twelve-well plates at a density of 6 x 10 4 cells/well or 3 x 10 5 cells/well for ddMEFs and 18 ddHEKs, respectively. After 24 hours, fresh doxycycline was prepared at 1 mg/mL, diluted, and then added to 19 the cells at the appropriate concentrations to induce a range of HA-MBNL1 or GFP-MBNL1 protein expression 20 (0 -10 ng/mL for ddHEKs, 0 -1000 ng/mL for ddMEFs). To induce mOrange-RBFOX1, fresh ponasterone-A 21 (Thermo Fisher Scientific) was dissolved in 100% ethanol and then added to the media (final concentration = 22 10 μM). Cells were treated with 100% ethanol as a vehicle control for -RBFOX1 conditions. Cells were then 23 harvested 24 hours post-drug treatment for experimental analysis.  Table 5. PCR-amplified cDNA samples were visualized and quantified using the Fragment 36 Analyzer DNF-905 dsDNA 905 reagent kit, 1-500bp (Agilent Technologies) and associated ProSize data 37 analysis software. Quantified Ψ values were plotted against relative MBNL levels as determined by immunoblot 38 or log (doxycycline (ng/mL)). Curve fitting was performed using the GraphPad Prism software using non-linear 39 curve fitting (log(agonist) vs. response -variable slope (four parameters)) (Ψ = Ψmin+ ((Ψmax-Ψmin) / (1 + 40 10 ((log(EC50) -log[MBNL1]) * slope) ). Where appropriate the minimum and maximum values were restricted to fall 41 between 0 and 1, respectively. Parameters that correlate to biological data, such as log (EC50) and slope, were 42 then derived from these curves. 43

44
Protein expression and purification 45 The RBFOX2 RRM expression vector (XP_016884187, amino acids 100-194) was gifted by of a GST-SBP tandem affinity tag. The RRM protein was recombinantly expressed using BL21 Star (DE3) cells 48 (Invitrogen) and protein expression induced using 0.5 mM IPTG at an OD600 = 0.6-0.7 for 2 hours at 37°C. 49 Following induction, cells were lysed in B-PER bacterial protein extraction reagent (Thermo Fisher Scientific), 50 supplemented with TURBO DNase I (5 U/mL, Thermo Fisher Scientific) and lysozyme (100 mg/mL) for 30 51 minutes at room temperature. The lysate was then diluted with 1 volume of 1X PBS, incubated on ice for 30 52 minutes, and centrifuged at 15,000 rpm for 15 minutes. The supernatant was isolated, filtered, and protein 53 loaded onto a GSTrap FF 5 mL column (GE Life Sciences, Marlborough, MA, USA) using the NGC 54 Chromatography system (Bio-RAD). The column was washed twice with 5 column volumes (CV) of 1X PBS 55 followed by 5 CVs of 1X PBS + 1M NaCl. The GST affinity tag was then cleaved using the PreScission 56 Protease (GE Life Sciences) for four hours in cleavage buffer (50 mM Tris pH 7.5, 150 mM NaCl, 1 mM EDTA, 57 1 mM DTT). The cleaved protein product was then eluted from the column and dialyzed overnight at 4°C into 58 storage buffer (500 mM NaCl, 25 mM Tris pH 7.5, 5 mM BME, 50% glycerol) using a 7000 MWCO Slide-A-59 Lyzer Dialysis Cassette (Thermo Fisher Scientific). Sample purity was assessed via SDS-PAGE gel and 60 working concentrations determined via Pierce 660nm protein assay reagent and BSA standards (Thermo 61 Fisher Scientific). 62 N-terminal GST-tagged MBNL1 (amino acids 1-260) 59 was recombinantly expressed using BL21 Star 63 (DE3) cells (Invitrogen) and protein expression induced using 0.5 mM IPTG at an OD600 = 0.6-0.7 for 2 hours at 64 lysate was then added to GST-beads (Thermo Fisher Scientific) pre-equilibrated in binding buffer and 69 incubated for 1 hour at 4°C. The beads were then washed 3 times with 20 bead volumes of binding buffer by 70 spinning beads at 750 x g for 3 minutes at 4°C and the supernatant removed. For the third wash, beads were 71 incubated in binding buffer for 15 minutes at 4°C. GST-MBNL1 was eluted from the beads three times by 72 incubating the beads with 2 bead volumes of elution buffer (50mM Tris pH 8.0, 150mM NaCl, pH 8.0, 10mM 73 reduced glutathione) followed by spins as described above. Elutions were combined and dialyzed for 2 hours 74 at 4°C in storage buffer (50mM Tris pH 8.0, 150mM NaCl, 1mM DTT, 50% glycerol) using a 7000 MWCO 75 Slide-A-Lyzer Dialysis Cassette (Thermo Fisher Scientific). Dialyzed GST-MBNL1 was aliquoted, flash frozen 76 in liquid nitrogen, and stored at -80°C. Sample purity was assessed via SDS-PAGE gel and working 77 concentrations determined via Pierce 660nm protein assay reagent and BSA standards (Thermo Fisher 78 Scientific). 79 80 Electrophoretic mobility shift assays 81 Purified GST-MBNL1 and RBFOX2-RRM proteins as described above were used to conduct the EMSA 82 assay. Fluorescent RNA oligonucleotide probes with a 3'-fluorescein tag were ordered from Integrated DNA 83 Technologies (IDT). Equilibration working solution was made fresh by combining 1mL of equilibration buffer 84 stock (0.01% IGEPAL CA630, 10 mM Tris, pH 8.0, 100 mM NaCl) with 131.6 L of 100% glycerol and 0.01 85 mg/mL yeast tRNA (Invitrogen). Serial dilutions of GST-MBNL1 were performed with RNAse/DNase free water 86 on ice. The INSR fluorescent RNA probe was heated for 5 min at 60C and allowed to cool for 2 min. The 87 binding reaction was then assembled in order as follows: 5uL of 10 nM fluorescent RNA probe, 43 μL of 88 equilibration working solution, 2 μL of RBFOX2-RRM (192 nM final concentration), 2 μL of serial dilution GST-89 MBNL1 and let the binding reaction sit covered from light at room temperature for 30 min. After incubation the 90 binding reaction was loaded onto Criterion 26 well, 10% polyacrylamide TBE precast gels (Bio-Rad) and run in 91 gel was transferred onto a low-fluorescence PVDF membrane (Bio-Rad) for 1 hour at 50V in 1X transfer buffer 02 (25mM Tris pH 8.3, 192mM glycine, 20% methanol (v/v)). The membrane was blocked for 1 hour at room 03 temperature using SeaBlock (ThermoFisher Scientific) and incubated overnight at 4°C with a primary antibody 04 targeting RBFOX1 [1:5000 RBFOX1, ab183348 (Abcam, Cambridge, MA, USA)]. Membranes were then 05 incubated for 1 hour at room temperature with a secondary antibody [1:15,000 donkey anti-mouse 800CW (LI-06 COR Biosciences, Lincoln, NE, USA)]. Images were acquired using an Odyssey CLx imager (LI-COR). Blots 07 were then stripped using Restore Fluorescent Western Blot Stripping Buffer (Thermo Fisher Scientific) as per 08 the manufacturer's protocol. Effective antibody removal was assessed via re-imaging of blot. Next, blots were 09 re-blocked in SeaBlock for 1 hour at room temperature and then incubated overnight at 4°C with primary 10 antibodies [1:2,000 MBNL1, 4A8 (Millipore Sigma, Burlington, MA, USA), 1:1,000 GAPDH, 14C10 (Cell 11 Signaling Technology, Danvers, MA, USA)]. Membranes were then incubated for 1 hour at room temperature 12 with secondary antibodies [1:15,000 donkey anti-mouse 800CW (LI-COR), 1:15,000 donkey anti-rabbit 680RD 13 (LI-COR)]. Images were acquired using an Odyssey CLx imager and quantification performed using the 14 associated ImageStudio Lite software (LI-COR). Relative MBNL1 levels in the presence or absence of 15 RBFOX1 co-expression were calculated by first normalizing lanes within the same gel to GAPDH and then by 16 normalizing levels of MBNL1 to the average relative MBNL1 expression level of untreated cells (no MBNL1 or 17 RBFOX1 expression). 18 harvested utilizing the Aurum Total RNA Mini kit (Bio-Rad) and DNase treated on column. RNA quality was 23 assessed using the Fragment Analyzer DNF-471 standard sensitivity RNA analysis kit (Agilent Technologies). 24 RNAseq libraries (500 ng, RQN > 9) were prepared using the NEBNext Ultra II Directional RNA Library Kit for 25 Illumina with NEBNext rRNA depletion (Bio-Rad). Protocols were followed as per the manufacturer's guidelines 26 except that 40X adaptor dilutions were utilized, 4X lower concentrations of index primers were used, all bead 27 incubations were performed at room temperature, and 10 cycles of library amplification were performed. 28 Library size and quality was assessed using the Fragment Analyzer DNF-474 High Sensitivity NGS Fragment performed using the Illumina NextSeq 500 (Illumina, San Diego, CA, USA). 32 33 Global splicing analysis 34 RNA sequencing reads were first cleaned using FASTQC and aligned to the Ensembl reference genomes 35 hg38 (human) and mm10 (mouse) using STAR 60 with additional flags (--outSAMtype BAM SortedByCoordinate 36 --quantMode GeneCounts --alignEndsType EndToEnd). Splicing analysis was conducted using rMATS (--cstat 37 0.01 --tstat 4) 61 . Customs scripts were used to further filter exclusively skipped exon events from rMATS output 38 for significant MBNL1-dependent events and to generate a dataset of common events between rMATS 39 pairwise comparisons. Group binning was done manually using the parameters described in the text via 40 Microsoft Excel. 41 42

Conservation 43
A complete list of human and mouse genes that are one-to-one orthologs was generated using BioMart 44 (release 101). A custom script was used to match one-to-one orthologous genes to events from the 9 splicing 45 events. Custom script was used to remove matches that did not have % translated sequence homology >= of RBFOX1 expression. Additionally, ΔΨ between -RBFOX1 and + RBFOX1 conditions are listed for the lowest (ΔΨLOW) and highest levels (ΔΨHIGH) of HA-MBNL1 expression. These values correspond to splicing dose-response curves presented in Fig 1c-e. Data represented as mean ± standard error where appropriate.
Supplemental Table 2: rMATS output with associate co-regulatory groups. Individual sheets for ddHEK and ddMEF data where columns A-O contain extracted rMATS output data containing both inclusion and exclusion SE events and columns P-AC contain co-regulatory group binning of each event using excel formulas using increasing cutoffs (ΔΨ = 0.1, 0.15, 0.2,

0.3)
Supplemental Table 3: Orthologous Events and Matching co-regulatory groups. List of human and mouse orthologous splicing events based on one-to-one gene orthologs. Events binned into groups are matched by chromosomal coordinates of upstream exon start and downstream exon end. Cells in rows that are blank contain duplicate events that have been already matched based on start and end exonic coordinates.
Supplemental Table 4: Quantitative parameters of MBNL1 dose-response curves for endogenous skipped exon events in presence or absence of RBFOX1 expression in ddHEKs and ddMEFs. Quantitative parameters that describe MBNL1 dose-dependent splicing regulation in presence or absence of induced RBFOX1 expression are listed including slope, log (EC50), and span (ΔΨ). Additionally, ΔΨ between -RBFOX1 and + RBFOX1 conditions has been quantified for lowest (ΔΨLOW) and the highest (ΔΨHIGH) levels of MBNL1 expression ( Fig   S2). The assigned co-regulatory group categorization (A-I) as derived from RNAseq (ΔΨ threshold = 0.10, Fig 3f) is also listed. Data represented as mean ± standard error.
Supplemental Table 5: Primers utilized for endogenous RT-PCR splicing analysis.
Table of primers sequences utilized for evaluation of endogenous exon inclusion levels via RT-PCR analysis in ddHEKs and ddMEFs. PCR annealing temperatures, cycle number, and inclusion/exclusion product sizes (bp) are also listed. Chromosomal coordinates of each skipped exon and flanking exons as identified via rMATS analysis are also provided (hg38 and mm10 for ddHEKs and ddMEFs, respectively).