|Home | About | Journals | Submit | Contact Us | Français|
RNA splicing plays a critical role in the programming of neuronal differentiation and, consequently, normal human neurodevelopment, and its disruption may underlie neurodevelopmental and neuropsychiatric disorders. The RNA-binding protein, fox-1 homolog (RBFOX1; also termed A2BP1 or FOX1), is a neuron-specific splicing factor predicted to regulate neuronal splicing networks clinically implicated in neurodevelopmental disease, including autism spectrum disorder (ASD), but only a few targets have been experimentally identified. We used RNA sequencing to identify the RBFOX1 splicing network at a genome-wide level in primary human neural stem cells during differentiation. We observe that RBFOX1 regulates a wide range of alternative splicing events implicated in neuronal development and maturation, including transcription factors, other splicing factors and synaptic proteins. Downstream alterations in gene expression define an additional transcriptional network regulated by RBFOX1 involved in neurodevelopmental pathways remarkably parallel to those affected by splicing. Several of these differentially expressed genes are further implicated in ASD and related neurodevelopmental diseases. Weighted gene co-expression network analysis demonstrates a high degree of connectivity among these disease-related genes, highlighting RBFOX1 as a key factor coordinating the regulation of both neurodevelopmentally important alternative splicing events and clinically relevant neuronal transcriptional programs in the development of human neurons.
The genetic programs that regulate neuronal differentiation and development are of fundamental importance. In addition to transcriptional levels of control, RNA splicing plays a critical role in establishing these genome-wide programs in a tissue-specific manner (1–3). Neuron-specific RNA splicing programs are therefore of considerable interest in understanding cellular differentiation and neurodevelopment (4). Examples include the polypyrimidine tract-binding protein (PTBP1), which plays a key role in regulating the switch between splicing programs leading to neuronal and non-neuronal cell differentiation (5), and the NOVA1 splicing factor, which directly regulates hundreds of neuron-specific splicing events involving such critical developmental processes as synaptogenesis and neuronal migration (6,7). As virtually all multi-exon genes undergo alternative splicing, with the majority occurring in a tissue-specific manner (2,8), a more detailed understanding of the regulation of such events is essential if we are to better understand neuronal differentiation.
RBFOX1 (also known as A2BP1 or FOX1) is a neuron-specific splicing factor that exerts both positive and negative regulatory effects on alternative splicing (9) and has been previously implicated in several neurodevelopmental and neuropsychiatric disorders including autism spectrum disorder (ASD), mental retardation and epilepsy (10–15), attention-deficit hyperactivity disorder (16), bipolar disorder, schizoaffective disorder and schizophrenia (17–19). Previous work using bioinformatic predictions has suggested that RBFOX1 affects a large tissue-specific splicing regulatory network (2,4,20,21) and this has been further supported by RNA-protein-binding studies of its homolog RBFOX2 (RBM9 or FOX2) in human embryonic stem cells (22), and studies in knockout mice which show altered synaptic transmission, increased membrane excitability and a predisposition to seizures (23). Together, these studies suggest that RBFOX1 is a high-level regulatory factor in the neuronal development cascade acting through the control of RNA splicing patterns and subsequent alterations in gene expression.
To address whether RBFOX1 regulates early genetic programs important for proper neuronal development and integrative coordinate functioning of neurons, we sought to identify the specific splicing events regulated by RBFOX1 by performing a genome-wide analysis using RNA sequencing in differentiated primary human neural progenitor (PHNP) cells. Here we identify the splicing regulatory network downstream of RBFOX1 and show that it regulates the alternative splicing of a large set of genes important for neuronal development, maintenance and proliferation. We further demonstrate that disruption of RBFOX1 function also leads to widespread transcriptional changes involving additional genes important to these pathways. Furthermore, a subset of these transcriptionally altered genes show coordinate expression and are implicated in the molecular pathogenesis of ASD and other neurodevelopmental conditions, suggesting that RBFOX1 function is also important for the correct regulation of a network of clinically relevant neuronal transcriptional programs in human neurodevelopment.
To study the role of RBFOX1 in human brain development, we used cultures of fetal PHNP cells (24–26) to avoid evolutionary or species-specific differences in RBFOX1 splicing regulatory programs (23). As RBFOX1 is itself highly alternatively spliced (9) (Fig. 1A), we performed quantitative real time-polymerase chain reaction (qRT-PCR) using primers directed to the first two constitutive exons of all RBFOX1 isoforms (Fig. 1B) to measure its expression in progenitors and differentiating neurons. We determined that undifferentiated PHNP cells do not produce significant amounts of RBFOX1. However, upon differentiation into neurons, there was a substantial increase in gene expression, including isoforms with an active RNA-binding domain (27) (Fig. 1B). Consistent with this differentiation-induced increase in expression, in situ hybridization of RBFOX1 in a 19-week human fetal brain shows the majority of expression to be in regions harboring post-mitotic neurons, with a reduced signal in the germinal zones relative to the basal ganglia and cortical plate (Fig. 1C). This is similar to what is seen in developmental mapping within the mouse brain, suggestive of expression in post-mitotic projection neurons and interneurons (28), implicating RBFOX1 in the development and maturation of human neurons.
To further validate the use of the PHNPs, we assessed whether the RBFOX1 isoform profile was comparable with that of the fetal or adult human brain (Fig. 1D). Similar isoform expression, including both nuclear and cytoplasmic forms (29), was noted among the human fetal brain, adult brain and the undifferentiated and differentiated PHNPs, indicating that RBFOX1 isoform expression in this primary cell culture system is similar to that of the human brain in vivo. Furthermore, PHNPs recapitulate brain-specific inclusion of exon 16 and exclusion of muscle-specific exon 17 (30) (Fig. 1D).
To better study its role in neuronal development, we reduced all RBFOX1 isoforms in PHNPs through RNA interference by introducing a short hairpin RNA (shRNA) against exon 9 (shRBFOX1) (Supplementary Material, Fig. S1A). PHNPs were then differentiated over a period of 4 weeks as RBFOX1 mRNA levels appear to be maximal at this point (Fig. 1B), suggesting a functionally important time point. Transduction into PHNP cells resulted in approximately a 50% reduction of mRNA levels (data not shown) and a 60–70% reduction of RBFOX1 protein after 4 weeks of differentiation compared with the non-targeting hairpin (Supplementary Material, Fig. S1B). This system recapitulates the peripheral blood levels of RBFOX1 found in clinically haploinsufficient patients with autism and developmental delay (10,12) and the reduction observed in brains from patients with ASD (15).
To assess genome-wide effects of RBFOX1 knockdown on maturing human neurons, RNAs from five biological replicates transduced either with shRBFOX1 or shGFP were analyzed using next-generation sequencing. Analysis of alternative splicing was performed according to published methods (8,31). We identified 996 significant alternative splicing events involving 603 unique genes (Supplementary Material, File S1).
As our methods will detect alternative splicing of genes regulated by RBFOX1 both directly and indirectly, we examined the intronic sequence flanking the most significant alternatively spliced exons for the presence of the RBFOX1-binding sequence (U)GCAUG (4,9,32,33) (Fig. 2B, Supplementary Material, File S2). Overall, the canonical RBFOX1 site was detected in 56% of the 996 total significant alternative splicing events, suggesting that other splicing factors may contribute to the observed alternative splicing changes or that additional novel RBFOX1-binding sites may exist. The presence of the RBFOX1 site did not correlate with the magnitude of the splicing effect (Supplementary Material, File S2). We also found that the RBFOX1-binding site was enriched above baseline only in the downstream introns (Fig. 2B, Supplementary Material, File S2), similar to what has previously been seen with genes showing altered alternative splicing in autistic brains with reduced RBFOX1 levels (15).
We next examined whether the locations of the observed RBFOX1 sites were consistent with the alternative splicing changes detected by RNA sequencing (Fig. 2B, Supplementary Material, File S2). After excluding alternative exons that possessed RBFOX1 sites both upstream and downstream, only 51% of the remaining sites were positioned consistently with the observed splicing change, suggesting this canonical model (9,33) may not be universal to all alternative exons. Interestingly, 63% of all correctly positioned RBFOX1 sites were found downstream of the alternative exon (Fig. 2B, Supplementary Material, File S2), suggesting that the majority of downstream RBFOX1 sites likely do act as enhancers, whereas upstream sites may act as silencers or enhancers depending on local context.
In mice, Rbfox1 and Nova1 coordinate alternative splicing during neurodevelopment (6), so we examined the intronic sequences surrounding all the alternatively spliced exons in our data set for the enrichment of the NOVA1-binding site, YCAY (4,6). We found a bias against the site in the upstream introns and no significant enrichment in the downstream introns where RBFOX1 sites are enriched (Fig. 2B). In contrast, binding sites for the unrelated splicing factor, PTBP1 (4), were found to be enriched in both upstream and downstream introns (Fig. 2B), likely consistent with its role in the repression of neuronal exons in non-neuronal tissues (5). On more detailed inspection, we did find upstream enrichment of alternate sequences that NOVA1 has been reported to bind (34) (Supplementary Material, File S2), but in vivo cross-linking and immunoprecipitation experiments do not suggest these to be of functional significance (35). Therefore, these data suggest the splicing events observed here are predominantly regulated independent of the NOVA1 pathway. To test this prediction, we chose a list of Nova1 splicing targets derived from the mouse brain (6), where an estimated 15% may be under coordinate regulation with Rbfox1/2 (6), and found high concordance (Fig. 2D, Supplementary Material, File S3) indicating that, despite a lack of enrichment for the NOVA1 core binding site (Fig. 2B), many of these genes are regulated by NOVA1, likely through splicing events distinct from those observed here. This observation highlights the independent yet coordinative effects of RBFOX1 and NOVA1 on human neuronal development.
Because we did not observe RBFOX1- or NOVA1-binding sites near many of the genes we found to be alternatively spliced, we explored the hypothesis that RBFOX1 might regulate the alternative splicing or differential expression of other splicing factors that could in turn affect downstream splicing events. We found that three splicing/RNA-processing factors (34) were alternatively spliced in our data set of predicted direct RBFOX1 targets: heterogeneous nuclear ribonucleoprotein D (HNRNPD), known to regulate gene expression in the developing brain (36), and HNRNPH1 and HNRNPA1, proteins that can modulate splice site selection both independently and in collaboration (37) (Supplementary Material, File S3). In our data set, HNRNPD sites were found to be enriched downstream of 29% of the alternative exons (285 of 996 events), whereas HNRNPH1 was enriched upstream of 21% (205 of 996 events) (Supplementary Material, File S2). We also found that the splicing factor ELAVL2, thought to play a role in neuronal differentiation (38), was differentially expressed. ELAVL2-binding sequences were enriched both upstream and downstream of alternatively spliced exons in our data set, with bindings sites present in 36% (355) and 31% (305) of introns, respectively (Supplementary Material, File S2). Together, binding sites for either differentially spliced or differentially expressed splicing factors were found in the flanking introns of 67% of alternative exons (666 of 996) with 480 upstream events (48%) and 411 downstream events (41%) (Supplementary Material, File S2). More than half of these (56%, 374 of 666) also had a flanking RBFOX1 site but 29% overall (292 of 996) did not, indicating that although the majority of events are likely direct, regulatory factors downstream of RBFOX1 likely contribute to some of the splicing changes seen following RBFOX1 knockdown.
To validate these RBFOX1-dependent splicing results, we performed quantitative RT-PCR (qRT-PCR) and obtained a supportive confirmation rate of 68% for a set of 25 splicing events (Fig. 2C). The presence of residual RBFOX1 in these cell lines (Supplementary Material, Fig. S1B) likely influenced these results, as many of the observed RNA sequencing changes were less than 2-fold (Fig. 2C) and this level has been previously noted to reduce splicing validation by this method (39). Therefore, to independently validate our results, we compared the RBFOX1-dependent gene set with various relevant published gene lists, including genes bioinformatically predicted to be regulated by RBFOX1 (20) or regulated by its homolog RBFOX2 (22), as well as genes showing altered alternative splicing in autistic brains with reduced RBFOX1 levels (15). Our set showed highly significant concordance rates with all these lists, and no enrichment when compared with control lists of genes from cellular organelles (Fig. 2D, Supplementary Material, File S3).
We hypothesized that the genes within this splicing network were likely involved in key neuronal pathways which might be regulated as part of the cellular role of RBFOX1. Gene ontology analysis revealed that genes with increased exon inclusion are involved in biological processes related to gene expression, such as cellular protein metabolism, post-transcriptional modification, translation and RNA processing, including RNA splicing (Table 1, Supplementary Material, File S4). Gene ontology analysis of the corresponding subset of genes with increased exon skipping revealed involvement in neuronal development and cytoskeletal organization (Table 1, Supplementary Material, File S4), processes related to cell differentiation. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis revealed similar results (Supplementary Material, File S4), as did restricting the gene ontology analysis to those genes most likely to be direct targets of RBFOX1 (Supplementary Material, File S4).
Since previous analysis of gene expression in ASD brains identified RBFOX1 as a hub in a gene co-expression module enriched for genes involved in synaptic transmission (15), synaptic transmission is altered in Rbfox1 knockout mice (23) and a number of autism-related genes are involved in synaptic function (40) (SFARI database, http://genes.sfari.org), we examined whether genes whose alternative splicing was affected by RBFOX1 knockdown were also enriched in the human post-synaptic density (PSD). We compared proteomic analysis of the murine PSD (41) with the RBFOX1 splicing targets identified here and found that nearly a quarter of the RBFOX1 target set is enriched in PSD proteins (136 genes, 23% of total targets) (Fig. 2D; Supplementary Material, File S3), indicating a role of RBFOX1 in the splicing regulation of genes involved in synaptic maturation and function (see also Supplementary Material, Analysis).
Given that altered splicing in the RBFOX1-depleted PHNP cells affects factors related to gene expression (Table 1, Supplementary Material, File S4) and evidence from autistic brains with reduced RBFOX1 levels suggests that transcriptional programs play an important role in neurodevelopmental disease (15), we next examined how gene expression, irrespective of splicing, was affected by RBFOX1 knockdown. We identified 981 significant differentially expressed genes in the RBFOX1 knockdown cell line (Supplementary Material, File S5). Validation was performed using qRT-PCR and a confirmation rate of 86% was obtained for a sample of 44 differentially expressed genes (Fig. 3A).
Given that RNA transcription and processing are coupled and alterations in the splicing of a transcript can influence its expression (42,43), we compared the RBFOX1-dependent transcriptional and splicing data sets and found that the differentially expressed genes regulated downstream of RBFOX1 were not alternatively spliced (Fig. 3B, Supplementary Material, File S6). Gene ontology analysis revealed that genes with increased expression upon RBFOX1 knockdown are related to cell migration, signaling and proliferation (Table 2, Supplementary Material, File S7), whereas genes with decreased expression are involved in neuronal development, synapse function and cell differentiation (Table 2, Supplementary Material, File S7), remarkably similar pathways to those seen for the genes whose alternative splicing was dependent on RBFOX1. KEGG pathway enrichment was also similar (Supplementary Material, File S7). Despite the overlap in functional pathways, the transcriptionally regulated gene set did not overlap with RBFOX1-dependant splicing targets identified here or elsewhere (Fig. 3B, Supplementary Material, File S6). Together this implies that similar functional pathways and programs are being regulated by these two processes, splicing and transcription, and that although these RBFOX1-regulated transcriptional and splicing networks are distinct from one another, they converge on common biological processes important to neuronal differentiation.
As RBFOX1 has been implicated in several neurodevelopmental and neuropsychiatric disorders (10–19), we hypothesized the mechanism could be via the missplicing of other disease genes, as has been previously suggested (15). We compared our differential splicing list with lists of genes associated with spinocerebellar ataxia, epilepsy and autism but found these represented no more than 1–2% of the genes within the set (Fig. 2D, Supplementary Material, File S3). Despite haploinsufficiency of RBFOX1 previously being shown to cause autism and ASD (10,12), we notably did not find a significant enrichment for alternatively spliced genes currently predicted to be related to autism (40) (SFARI database, http://genes.sfari.org) (Fig. 2D, Supplementary Material, File S3). This lack of strong overlap indicates that the relationship between RBFOX1 function and neurodevelopmental disease is unlikely to occur primarily due to the direct RBFOX1-controlled missplicing of known disease genes, as was suggested by previous analysis of the ASD brain (15).
In contrast to the differentially spliced genes, those with differential expression were 3-fold enriched in known disease genes (Figs 2D and and3B,3B, Supplementary Material, Files S3 and S6), suggesting altered transcriptional regulation could play a role in RBFOX1's association with neuropsychiatric disease. Among these differentially expressed genes were several well-supported ASD candidates including DHCR7, GABRB3, ITGB3, NRXN1 and MET (40) (Supplementary Material, File S6). This subset of transcriptionally regulated genes implicated in ASD suggests an important role for RBFOX1 in the regulation of critical transcriptional programs required for proper human neurodevelopment that can increase risk for common neurodevelopmental disorders when disrupted.
To better understand the systems level organization of the transcriptional network altered by RBFOX1 knockdown, we performed weighted gene co-expression network analysis (WGCNA) (15,44–46). To identify the modules most related to RBFOX1 function, we assessed the module eigengene relationship with RBFOX1 knockdown status (47). Two modules were identified that showed the largest difference between RBFOX1 knockdown and control, designated as the blue and yellow modules, respectively (Fig. 4, Supplementary Material, File S8).
Since the differentially expressed data set was enriched in genes with established relationships to neurodevelopmental disease, we assessed the top RBFOX1 modules for the enrichment of ASD genes. The blue module consisted of 737 genes (Fig. 4A, Supplementary Material, File S8) and its members include the transcription factor FOXP2, which functions in the regulation of genes important to CNS development and language (45), as well as several other genes linked to ASD (CADPS2, DMD, HLA-DRB1, INPP1, ITGB3, LRFN5, MDGA2, NTRK3, OPRM1, PARK2 and SDC2) (48).
The yellow module contained 312 genes (Fig. 4B, Supplementary Material, File S8) and its hub genes include DLG1, a protein important in modulating cytoskeletal interactions to establish cell polarity (49), and implicated as a potential candidate gene in both ASD and schizophrenia (50,51). The yellow module contained several other ASD candidate genes including CNTN3, DPYD, IMMP2L, MACROD2, MET, RAPGEF4, REEP3, SLC4A10 and TPH2 (48). Of these ASD genes, six are differentially expressed following RBFOX1 knockdown (DLG1, MET, IMMP2L, REEP3, SLC4A10 and TPH2) and represent a key subset of the yellow module with a statistically significant (P ≤ 1 × 10−4) greater average connectivity and shorter mean path length relative to all other module members (Supplementary Material, File S8). These observations show that RBFOX1 functionally contributes to the transcriptional regulation of important neurodevelopmental cellular programs that are clinically relevant to autism. This confirms previous indirect studies in human brain, suggesting a role for RBFOX1 in regulating ASD-related transcription (15).
Here we demonstrate that the neuronal splicing factor RBFOX1 regulates the alternative splicing of an extensive network of genes involved in neuronal differentiation and maintenance. We examined these alternatively spliced genes and observed enrichment for pathways involved in cellular proliferation, cell signaling, cytoskeletal organization and cell adhesion, gene expression and regulation and neuronal development (Table 1, Supplementary Material, File S4), suggesting that RBFOX1 may represent a molecular switch in neuronal progenitors which promotes correct differentiation into functional neurons (Fig. 5). In humans, RBFOX1 haploinsufficiency (similar to the level of knockdown achieved here) (Supplementary Material, Fig. S1B) causes neurodevelopmental disease such as ASD (10,12). Based on previous observations in autistic brains (15), we initially had predicted that genes undergoing differential alternative splicing would be related to ASD or similar clinical disorders. However, this was not observed. Instead, it was the network of genes that were differentially expressed in cells lacking RBFOX1 that were associated with neurodevelopmental disease. Interestingly, these differentially expressed genes were also related to similar biological pathways as those that were alternative spliced (Tables 1 and and2,2, Supplementary Material, Files S4 and S7). A subset of these differentially expressed genes were also highly connected in a module of co-expressed genes that were enriched in ASD candidates (48) (Fig. 4B, Supplementary Material, Files S6 and S8). The RBFOX1 protein is therefore involved in regulating both neuronal alternative splicing and coordinative transcriptional programs, both of which contribute to normal neurodevelopment (Fig. 5).
It is interesting to compare these data with recent studies, using microarray analysis in the brains of adult knockout mice lacking Rbfox1, identifying differentially spliced genes which impacted synaptic transmission and membrane excitability, leading to an increased susceptibility for seizure events (23). Similar to the mouse, a number of the RBFOX1-dependent alternatively spliced genes in our data set appear to be involved in synapse formation and function (Fig. 2D, Supplementary Material, File S3), areas of particular interest in the study of the pathogenesis of neurodevelopmental disease and autism (52,53). Despite this, we did not observe significant differences in the RBFOX1 knockout cell line with regard to basic cellular morphology (data not shown). We attribute this to our cells not being fully devoid of RBFOX1 (Supplementary Material, Fig. S1B), with residual activity potentially preventing observation of the more extreme cellular phenotypes, although still reflecting the underlying molecular pathogenesis affecting splicing and transcription. In this sense, the current model approximates human haploinsufficiency observed in patients, as suggested by the highly significant overlap between our alternative splicing gene data set and the genes with altered splicing identified in autistic brains with reduced RBFOX1 (15) (Fig. 2D, Supplementary Material, File S3). To extend this comparison further, we evaluated our list of RBFOX1-dependent splicing events for genes related to epilepsy and found altered splicing of FLNA, a cytoskeletal gene involved in neural migration and heterotopia (54), and SLC1A3, a glutamate transport protein associated with episodic ataxia and epilepsy (55) (Fig. 2D, Supplementary Material, File S3). Additional genes associated with epilepsy were detected in the differentially expressed RBFOX1-dependent gene set including DCX, GABRB3, GAD2, KCNQ2, SLC12A5, SV2B and SYN1 (Fig. 3B, Supplementary Material, File S6). This observation is similar to the increase in ASD genes in this transcriptional data set and consistent with the seizures observed in the knockout mice (23), as well as the clinical finding of seizure disorders in patients with RBFOX1 mutations (10,12).
In addition to HNRNPD, HNRNPH1 and HNRNPA1, 15 additional RNA-binding proteins (56) were identified in the alternatively spliced data set (Supplementary Material, File S3), which may contribute to the neurodevelopmental phenotypes associated with RBFOX1. These include DDX17, a developmentally essential factor involved in miRNA processing and transcriptional co-regulation in cell proliferation (57,58), and NPM1, implicated in neuronal cell proliferation (59). Several of these genes also show altered expression in neurodevelopmental disorders, including NONO and FUBP1 (60), RPL3 (61) and RPL7 and RPS17 (62). In addition to ELAVL2, another six RNA-binding proteins (56) were also differentially expressed (Supplementary Material, File S6).
We also found 17 transcription factors that were alternatively spliced, including TEAD1, whose binding site was enriched in the promoters of the differentially expressed gene set (Supplementary Material, Analysis and File S9). Several of these transcription factors are known to be involved in neurodevelopment, including ARID1B, which causes an ASD-associated syndrome when haploinsufficient (63,64); ARNT2, a gene critical to hypothalamic development and associated with autistic traits (65); RPL7, differentially expressed in mothers of ASD children (62); GTF2I, involved in the Williams–Beuren neurodevelopmental syndrome (66); MEIS2, involved in forebrain development (67); and NOTCH3, causative for the inherited leukoencephalopathy CADASIL (68) (Supplementary Material, File S3). Forty-two transcription factors (69) were also differentially expressed (Supplementary Material, File S6), including ETS1, whose binding site was enriched in genes found in the blue module (Supplementary Material, Analysis and File S9). Binding sites for two transcription factors with some circumstantial evidence for involvement in ASD (70,71), SP1 and RORA, were enriched among the differentially expressed genes (Supplementary Material, File S9), and transcriptional regulation may also occur among the differentially spliced genes (Supplementary Material, Analysis).
One ASD candidate gene which may be particularly relevant is the cell polarity gene DLG1 (49–51), which is differentially expressed with RBFOX1 knockdown (Supplementary Material, File S5) and is a hub in the yellow WGCNA module that is enriched for ASD genes (Fig. 4B, Supplementary Material, File S8). DLG1 is also differentially spliced in the RBFOX1 knockdown line (Supplementary Material, File S3) and in autistic brains (15) (Supplementary Material, File S3), and is also a target of Nova1 in the mouse brain (6) (Supplementary Material, File S3). Further study of DLG1 may provide valuable insights into the downstream effects of disrupting the RBFOX1 networks. NLGN1, another interesting ASD candidate gene (72), was also differentially expressed.
Because RBFOX1 was originally identified as an ataxin-2-binding protein (73,74), its role in cerebellar ataxia is another potential area of interest. Here we observed differential alternative splicing of the ataxia genes ATXN2, ATXN10, FMR1 and SLC1A3 (Fig. 2D, Supplementary Material, File S3). Furthermore, we identified RBFOX1-dependent differential expression of several other ataxia genes (Fig. 3B, Supplementary Material, File S6), including CA8 and VLDLR, both of which cause a developmental syndrome characterized by cerebellar hypoplasia, congenital ataxia and mental retardation (75,76). CA8 further influences cell signaling via the function of the inositol triphosphate receptor, ITPR1. Mutations in ITPR1 also cause spinocerebellar ataxia, type 16 (76,77), and ITPR1 was also found to be differentially expressed in the RBFOX1-deficient cells (Fig. 3B, Supplementary Material, File S6). Interestingly, CA8 was found in the blue WGCNA module, and ITPR1 was in the yellow module (Fig. 4, Supplementary Material, File S8), indicating a close relationship between RBFOX1 and these cerebellar ataxia genes. Whether certain variants in RBFOX1 could lead to clinical ataxic syndromes independent of neurodevelopmental phenotypes has yet to be determined.
Our data reflect both direct and indirect regulatory events resulting from the loss of RBFOX1 protein. The finding that several RNA-processing factors and transcription factors are indeed alternatively spliced (Supplementary Material, File S3), and that other factor-binding motifs are enriched in the introns and promoters of the differentially spliced and expressed genes, respectively (Supplementary Material, Files S2 and S9), points to a complex regulatory system, likely involving multiple cis- and trans-acting elements, that is disrupted by the loss of RBFOX1 (Fig. 5). One key question will be whether novel intronic sequences found to be enriched in the differentially spliced data set (data not shown) represent specific RBFOX1-associated regulatory elements. The observation of enriched motifs in the 3′-untranslated regions and promoter regions of the RBFOX1 splicing network including some, like SP1, shared with the RBFOX1 transcriptional network (Supplementary Material, File S9) suggests that certain regulatory factors and/or signaling pathways may converge to coordinately regulate both networks in the parallel execution of the RBFOX1-mediated developmental program (Fig. 5). Elucidation of the factors mediating these processes in both cell lines and in the human brain (15), and how the loss of RBFOX1 impacts their function, will have important relevance to our understanding of neurodevelopment and the pathogenesis of neurodevelopmental diseases.
The RBFOX1 short-RNA hairpin targeting exon 9 was originally derived from clone ID TRCN0000073260 (Open Biosystems, Thermo Fisher Scientific, Rockford, IL, USA) expressed in the pLKO.1 vector. As use of this vector in our PHNP cultures can lead to vector-related cell death regardless of the RNA hairpin expressed, this hairpin was redesigned in the shRNAmir-30 context and subcloned in the pGIPZ vector (Open Biosystems, Thermo Fisher Scientific) at the BamHI and XhoI sites as per the manufacturer's recommendations. For viral production, the hairpin was transferred to the pLCR vector (78) at the EcoRI and XhoI sites. The RBFOX1 overexpression vector was generated by shuttling RBFOX1 transcript 4 (NM_018723.3) from pENTR-A2BP1 (Addgene, Cambridge, MA, USA) to pPELTR5.5.2-GA, using Gateway cloning (Invitrogen, Life Technologies, Carlsbad, CA, USA) to create pPELTR5.5.2-GA-RBFOX1. The pPELTR5.5.2-GA vector is a second generation lentiviral vector where the expression of the insert is under the control of an EF1-HTLV fusion promoter and includes the addition of an N-terminal HA epitope tag. A red fluorescent reporter protein is also expressed in cis with the insert with an intervening T2A sequence to allow for co-translational cleavage of the two proteins.
Normal human tissues were obtained from the NICHD Brain and Tissue Bank for Developmental Disorders at the University of Maryland School of Medicine (Baltimore, MD, USA; NICHD contract numbers N01-HD-4-3368 and N01-HD-4-3383) (fetal heart age 19 weeks and fetal brains age 18 and 19 weeks). The role of the NICHD Brain and Tissue Bank is to distribute tissue and therefore cannot endorse the studies performed or the interpretation of results. Additional tissues were obtained from the UCLA Center for AIDS Research Gene and Cellular Therapy Core (fetal brain age 17 weeks) and the UCLA Alzheimer's Disease Center Neuropathology Core (adult brain age 97 years).
PHNP cells were cultured from a 17-week human fetal brain, proliferated and differentiated as previously described (26). Genotyping was performed to ensure cells did not harbor major structural aberrations (24). Biological replicates (7–10 per condition) were grown on individual 10 cm plates and all transductions and other sample preparations were handled independently. For shRNA transduction, high-titer virus production was performed by the UCLA Jonsson Comprehensive Cancer Center Vector Core (Los Angeles, CA, USA), utilizing the specified lentiviral vectors, and viral transduction was performed using standard protocols.
HeLa cells (ATCC, Manassas, VA, USA) were grown under recommended conditions. Cells were transfected with the specified vectors using the FuGene 6 transfection reagent (Roche, Indianapolis, IN, USA) according to the manufacturer's recommended protocol. After 48–72 h, cells were harvested and protein was extracted for western blotting.
Cells were harvested and RNA was extracted using the miRNeasy Mini Kit (Qiagen, Valencia, CA, USA) as per the manufacturer's recommended instructions. For PHNP cells, protein was isolated during the RNA extraction procedure according to the supplemental protocols available at the manufacturer's website (www.qiagen.com). For HeLa cells, whole-cell protein lysates were prepared by directly lysing cells in 0.5% Nonidet P-40 buffer containing protease and phosphatase inhibitors. For tissues, RNA and protein were prepared from 30 mg tissue using the AllPrep DNA/RNA/Protein Mini Kit (Qiagen) according to the manufacturer's instructions. All protein concentrations were determined using the Bradford (Bio-Rad, Hercules, CA, USA) or BCA (Pierce, Thermo Fisher Scientific) protein assays.
RNA was initially treated with amplification-grade DNase I (Invitrogen, Life Technologies) and then converted to cDNA, using the SuperScript III First-Strand Synthesis System Kit (Invitrogen, Life Technologies). The resulting cDNA (25–50 ng/reaction) was utilized for qRT-PCR, using Bio-Rad iTaq SYBR Green Supermix with ROX (Bio-Rad) according to the manufacturer's recommendations. Reactions (10 µl volume) were performed in triplicate using either an ABI 7900HT (Applied Biosystems, Life Technologies, Carlsbad, CA, USA) or a LightCycler 480 (Roche) qPCR machine. Data analysis was performed using the comparative Ct method normalized to β-actin. Reactions where the replicates deviated by more than 0.2 Ct from the mean were excluded from further analysis. Each data point represents a minimum of two individual reactions. For confirmation of differential gene expression and splicing from transduced PHNP cell lines, an equal mixture of RNA from five biological replicates was used as a template to minimize bias from an individual sample. All primers were designed to cross an exon–exon junction and, for the detection of splicing events, this would span the junction between the alternative exon and one of its neighboring exons.
For semi-qRT-PCR analysis, 10 ng of cDNA was amplified within the linear range, using Accuprime Pfx Supermix (Invitrogen, Life Technologies) with primer pairs designed to detect products with and without the alternatively spliced exon. The program used was 94°C for 15 s, 55°C for 15 s and 72°C for 15 s for an average of 37 cycles. Products were separated on a 3.5% agarose gel, stained with ethidium bromide, photographed and analyzed by densitometry, using the publically available software ImageJ (http://imagej.nih.gov/ij/).
RBFOX1 genomic organization nomenclature used in this report is as published in Underwood et al (9). For isoform analysis, RNA was harvested from the indicated cultured cell lines and human tissues, converted to cDNA as above and RT-PCR was performed using PCR Supermix (Invitrogen, Life Technologies) with primers spanning exons 15–20, the region of highest alternative splicing diversity. The resulting products were subcloned using the TOPO TA cloning kit (Invitrogen) and 30–50 individual inserts were sequenced.
Sagittal and coronal sections of human fetal brain, age 19 weeks, were probed using an S35-labeled antisense riboprobe directed against exons 8–13 of RBFOX1 as previously described (79). The sense riboprobe was used as the control.
Protein (50–100 µg) was loaded onto 10–12% SDS–polyacrylamide gels, separated by electrophoresis, transferred onto Immun-Blot PVDF membranes (Bio-Rad), blocked in 5% milk, probed with primary antibody diluted into TBST for 1 h to overnight at 4°C and then probed with a secondary antibody diluted into TBS for 1 h to overnight at 4°C. Detection was performed using SuperSignal West Pico chemiluminescent substrate (Pierce, Thermo Fisher Scientific) according to the manufacturer's instructions and membranes were exposed to film. Band quantitation was performed from scanned images, using publically available densitometry software and all samples were normalized to β-actin levels.
Primary antibodies were against HA-tag (1:2000, Covance, Princeton, NJ, USA), RBFOX1 (rabbit polyclonal, 1:500, Aviva Systems Biology, San Diego, CA, USA) or β-actin (mouse monoclonal, 1:250 000, Sigma-Aldrich, St Louis, MO, USA). Secondary antibodies used were goat anti-rabbit horseradish peroxidase (1:2000, Cell Signaling Technology, Danvers, MA, USA) and goat anti-mouse horseradish peroxidase (1:4000, Millipore, Billerica, MA, USA).
RNA sequencing at a 75 bp single-end read scale was performed on an Illumina Platform G2 sequencer (Yale Center for Genome Analysis, West Haven, CT, USA), using polyA-enriched RNA from five biological replicates of the shGFP and shRBFOX1 PHNP cell lines differentiated for 4 weeks. Unique sequence reads were filtered for repetitive and ribosomal RNA sequences and then aligned to the human genome (hg19, NCBI 37), using the Burrows-Wheeler alignment tool (BWA) (80), allowing a maximum of three mismatches to generate reads mapping to exons (54 779 644 ± 1 258 444 for shGFP and 55 146 141 ± 2 027 402 for shRBFOX1) and exon junctions as defined by Refseq. Biological replicate counts were pooled for both gene expression analysis and alternative splicing analysis. The annotation used for mapping potential alternatively spliced reads consisted of a previously published genome-wide alternative gene expression database source (39). To be counted as a junction, a read had to overlap two adjacent exons by a minimum of 5 nucleotides and total counts had to equal or exceed 20. Splicing events were defined as C1–A1A2–C2, where the alternative exon is A1A2 (A1, 5′ 70 nucleotides; A2 3′ 70 nucleotides), and C1 and C2 represent the adjacent 70 nucleotides of the upstream and downstream exons, respectively. The corresponding junctions were therefore defined as C1A1 and A2C2 for the alternative isoform and C1C2 for the constitutive isoform. Events were excluded if the counts for either [C1A1] or [A2C2] exceeded the other by >2-fold. Percent inclusion was calculated as (([C1A1] + [A2C2])/2)/([C1C2] + ([C1A1] + [A2C2])/2). Statistical comparisons for both differential splicing and differential gene expression were made using a negative binomial distribution model with the DESeq software package (81). A cutoff threshold of P ≤ 0.05 was used for significance. Examination of the most variant gene expression indicated a reproducible difference between the shRBFOX1 and the shGFP lines, but also between transduced and wild-type cells as well (Fig. 2A). To correct for variation due to the transduction process, only the shGFP and shRBFOX1 lines were used for further comparisons. Gene ontologies and KEGG pathways were determined using the DAVID Bioinformatic resources (NIAID, NIH). Statistical comparisons to other gene lists were made using the hypergeometric probability. WGCNA was performed as previously described (45,47). Connectivity and shortest path calculation were performed using the igraph software package (82).
All RNA sequencing data are deposited in the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo/) (accession number GSE36710).
Intronic sequences flanking each significant alternative exon were interrogated for 400 nucleotides in either direction for the presence of splicing and RNA-processing factor binding sites. To calculate enrichment for sites assessed individually, an equivalent number of random intronic sequences were assessed from a pool of all genes in the human genome and this process was permuted 100 times to generate an average and standard deviation. Significance was calculated using the normal distribution. In addition, we also utilized a database of experimentally derived target RNA sequences (34) to generate position-weight matrices that were used to scan the intronic sequence using the Clover algorithm (83). For these searches, enrichment was calculated based on three different background data sets comprising all introns in the human genome, human CpG islands, and human chromosome 20 sequence.
For transcription-factor binding site analysis, 1000 base pairs upstream of the transcription start site of each gene were examined using experimentally defined position-weight matrices from the JASPAR database (84) with the Clover algorithm (83). The MEME algorithm (85) was employed for predicting site occurrence. Enrichment was calculated based on three different background data sets comprising all human gene promoters, human CpG islands, and human chromosome 20 sequence.
This work was supported by the National Institute of Mental Health (grants K08MH86297 to B.L.F., R37MH060233 to D.H.G., R01MH081754 to D.H.G., K08MH074362 to E.W., R00MH090238 to G.K. and T32MH073526 to N.P.); the Eunice Kennedy Shriver National Institute of Child Health and Human Development (grant p30HD004612 to F.G. and D.H.G.); the UCLA Caltech Medical Scientist Training Program (to N.P.); the Dr Miriam and Sheldon G. Adelson Medical Research Foundation (to C.V.); and the John Douglas French Alzheimer's Research Foundation (to E.W.).
The authors wish to thank Irina Voineagu, Benjamin Blencowe and Xinchen Wang for helpful comments regarding the splicing analysis; Daning Lu for technical assistance with cell culture; Hongmei Dong for technical assistance with in situ hybridization; and Emily Batton, Sonia Kantak and Emily Matthews for additional technical assistance.
Conflict of Interest statement. None declared.