Cell type-specific CLIP reveals that NOVA regulates cytoskeleton interactions in motoneurons
Genome Biology volume 19, Article number: 117 (2018)
Alternative RNA processing plays an essential role in shaping cell identity and connectivity in the central nervous system. This is believed to involve differential regulation of RNA processing in various cell types. However, in vivo study of cell type-specific post-transcriptional regulation has been a challenge. Here, we describe a sensitive and stringent method combining genetics and CLIP (crosslinking and immunoprecipitation) to globally identify regulatory interactions between NOVA and RNA in the mouse spinal cord motoneurons.
We developed a means of undertaking motoneuron-specific CLIP to explore motoneuron-specific protein–RNA interactions relative to studies of the whole spinal cord in mouse. This allowed us to pinpoint differential RNA regulation specific to motoneurons, revealing a major role for NOVA in regulating cytoskeleton interactions in motoneurons. In particular, NOVA specifically promotes the palmitoylated isoform of the cytoskeleton protein Septin 8 in motoneurons, which enhances dendritic arborization.
Our study demonstrates that cell type-specific RNA regulation is important for fine tuning motoneuron physiology and highlights the value of defining RNA processing regulation at single cell type resolution.
A thorough understanding of the complexities of the mammalian central nervous system (CNS) requires detailed knowledge of its cellular components at the molecular level. RNA regulation has a central role in establishing cell identity and function across the numerous cell types in the CNS [1,2,3,4,5,6,7]. Traditional whole tissue-based methods are particularly limited in their power to delineate cell type-specific RNA regulation in the mammalian CNS due to its vast cellular diversity and architectural complexity. While separation or induction of specific cell types in vitro provides a practical way for cell type-specific analyses [1,2,3, 8, 9], alteration of cellular biology due to loss of physiological contexts presents a significant caveat to this approach.
Recent technological breakthroughs using RiboTag and BAC-TRAP mouse lines have allowed for translational profiling at single cell type resolution [10,11,12,13]. These studies revealed remarkable differences in the population of translating mRNAs across various CNS cell types, highlighting the degree of molecular heterogeneity among neuronal cells. These methods offer important ways to study translated mRNAs in specific cell types, but do not provide a way to define other types of cell type-specific RNA regulation. Here we develop a complementary and more general means to study RNA processing and regulation in a cell type-specific manner.
RNA processing is regulated by RNA-binding proteins. Two of the best-studied are NOVA1 and NOVA2, neuron-specific KH-type RNA-binding proteins that bind to YCAY motifs and regulate alternative splicing and polyadenylation [14,15,16]. Using crosslinking and immunoprecipitation (CLIP), a method that allows stringent purification of protein–RNA complexes captured in vivo, we identified NOVA targets in mouse neocortex [16,17,18,19], and have estimated that NOVA participates in the regulation of ~ 7% of brain-specific alternative splicing events in mouse neocortex . Interestingly, NOVA targets are specifically enriched for transcripts encoding proteins with synaptic functions—a group of transcripts that drives CNS cell type diversity [13, 15].
Indeed, NOVA proteins are essential for the function of multiple neuronal cell types [21,22,23,24]. In particular, we previously uncovered a pivotal role for NOVA in maintaining spinal motoneuron (MN) survival and physiology [5, 21, 25]. Here, to more precisely define NOVA-regulated RNA processing in spinal MNs, we developed a new strategy combining BAC-transgenic mice and CLIP to identify MN-specific RNA regulation. NOVA targets in MNs were especially enriched for genes encoding microtubule-, tubulin-, and cytoskeletal protein-binding proteins. These results led us to uncover a NOVA-mediated RNA processing event differentially regulated in MNs involving a cytoskeleton protein, Septin 8, which controls dendritic complexity. Cell type-specific CLIP revealed a previously unidentified role of NOVA in MNs, highlighting the importance of cell type-specific analysis of RNA regulation.
BAC transgenic mice express epitope-tagged NOVA specifically in motoneurons
Nova1 and Nova2, the two mammalian Nova paralogs, encode proteins harboring three nearly identical KH-type RNA binding domains. Nova genes are widely expressed among various neuronal cell types in the spinal cord, including the MNs. For NOVA2, two isoforms arise from alternative translational start sites, which do not shown detectable RNA binding differences (Yano Y. and Darnell R. B., unpublished data). Our strategy to study MN-specific NOVA regulation is twofold—to express AcGFP (Aequorea coerulescens green fluorescent protein)-tagged NOVA specifically in MNs, followed by CLIP using antibodies against the AcGFP epitope tag.
To test if an N-terminal AcGFP tag alters NOVA RNA binding specificity, we compared global RNA binding profiles of untagged and AcGFP-tagged long NOVA2 isoform. NIH/3T3 cells ectopically expressing NOVA2 or AcGFP-NOVA2 were subjected to HITS-CLIP using antibodies against NOVA and GFP, respectively (Additional file 1: Supplemental material and methods; Additional file 2: Figure S1A). We generated complex NOVA:RNA interactomes for both tagged and untagged NOVA2, as defined by 1,775,101 unique CLIP reads for NOVA2 and 9,778,818 for AcGFP-NOVA2. Both sets of CLIP data showed comparable genomic distributions (Additional file 2: Figure S1B). We identified NOVA2 and AcGFP-NOVA2 CLIP peaks, defined as regions with significantly higher CLIP tag coverage than gene-specific background expected from uniform random distribution [26, 27]. The canonical NOVA recognition motif, YCAY, was enriched around NOVA2 and AcGFP-tagged NOVA2 CLIP peaks to the same degree (Additional file 2: Figure S1C, D; R2 = 0.94), with tetramers that may overlap YCAY by three or more nucleotides (YCAY, NYCA, CAYN) showing the highest enrichment at both NOVA2 and AcGFP-NOVA2 peaks (red dots in Additional file 2: Figure S1D). These indicate that the N-terminal AcGFP tag does not perceptibly alter NOVA2 binding to its cognate RNA motifs.
Having confirmed the RNA binding fidelity of AcGFP-tagged NOVA2, we generated transgenic mice expressing AcGFP-tagged NOVA2 long isoform under the MN-specific choline acetyltransferase (Chat) promoter. The AcGFP-Nova2 cassette was inserted into a bacterial artificial chromosome (BAC) harboring the Chat promoter through homologous recombination (Fig. 1a) . The engineered BAC was injected into fertilized C57BL/6 oocytes and two Chat:GFP-Nova2 BAC transgenic lines, #6 and #17, were selected and maintained from five original founder lines.
Chat:GFP-Nova2 transgenic mice were born at expected Mendelian ratios and were phenotypically indistinguishable from wild-type littermates throughout their lifespan (data not shown). Expression of the GFP-NOVA2 fusion protein was confirmed in the spinal cords of both Chat:GFP-Nova2 lines, but not the wild-type controls (Fig. 1b). To assess the MN-specific expression pattern of the transgene, we performed immunofluorescence staining on spinal cord transverse sections. GFP-NOVA2 expression was confined to CHAT-positive neurons in the transgenic lines (Fig. 1c). Thus, we have successfully established transgenic mouse lines with epitope-tagged NOVA specifically expressed in MNs.
HITS-CLIP generates a robust transcriptome-wide NOVA–RNA interaction map in motoneurons
To define NOVA binding sites in MNs, we undertook HITS-CLIP on Chat:GFP-Nova2 spinal cords using a mixture of two monoclonal antibodies against GFP to maximize avidity and specificity. GFP immunoprecipitation was performed using wild-type or transgenic spinal cords with or without UV crosslinking, followed by 32P-labeling of bound RNA. The presence of labeled GFP-NOVA2–RNA complexes was dependent on both UV crosslinking and the expression from the Chat:GFP-Nova2 transgene, further demonstrating the specificity of GFP-NOVA2 CLIP under these conditions (Fig. 2a, lanes 1–3). Partially digested RNA crosslinked to GFP-NOVA2 was subjected to subsequent library preparation steps (Fig. 2a, lane 4; “Methods”). Importantly, the production of cDNA was dependent on reverse transcriptase (Fig. 2b), confirming the absence of DNA contaminants in the CLIP RNA libraries. cDNA inserts between 30 and 80 nucleotides (nt; Fig. 2b) were selected for next-generation sequencing (NGS) and further analysis.
Four sets of biological replicate HITS-CLIP experiments were performed on each Chat:GFP-Nova2 transgenic line, with each biological replicate consisting of pooled spinal cord samples from five to seven 3-month-old mice. We obtained a total of 2,023,726 unique CLIP reads from these eight biological replicates. Since AcGFP-Nova2 expression is confined to MNs, we refer to this group of CLIP reads as MN NOVA CLIP reads. As a reference to all NOVA binding sites in the whole spinal cord, we additionally performed four standard NOVA CLIP experiments on endogenous NOVA1 and NOVA2 in 3-month-old wild-type mouse whole spinal cords (WSC), generating 17,353,049 unique WSC NOVA CLIP reads (Additional file 2: Figure S2A; Additional file 3A). CLIP reads from MN CLIP showed higher intronic and less exonic distribution compared to those from WSC CLIP (p = 4.2 × 10− 7; Additional file 2: Figure S2B).
NOVA peaks in MN and WSC CLIP datasets were defined separately . A total of 25,681 CLIP peaks (p ≤ 0.01) were identified for MN NOVA CLIP, which harbor CLIP reads originating from at least four of the eight biological replicates (biologic complexity (BC) ≥ 4 out of 8; Additional file 3 B) . In parallel, we identified 218,794 WSC NOVA peaks represented in at least two out of four biological replicates (BC ≥ 2 out of 4; Additional file 3C). The CLIP data from our two transgenic lines was highly correlated (R2 = 0.97), as were WSC CLIP reads (Fig. 2c), underscoring the ability of cell type-specific HITS-CLIP to reproducibly and quantitatively measure in vivo cellular protein–RNA interactions.
MN HITS-CLIP generated a genome-wide binding profile characteristic of endogenous NOVA proteins. The canonical NOVA binding motif YCAY was significantly and similarly enriched around MN and WSC CLIP peaks (Fig. 2d). WSC NOVA CLIP peaks mapped to 12,445 mm9-annotated Entrez genes, including to the alternative spliced regions of 303 out of 335 known NOVA-regulated genes (Additional file 3D); MN CLIP peaks mapped to 4450 Entrez genes, including to 174 known NOVA-regulated alternatively spliced regions (Additional file 3E) . It is of note that when compared to all NOVA-regulated genes in the WSC, this subset of MN NOVA targets are especially enriched for genes encoding microtubule-binding (GO:0008017), tubulin-binding (GO:0015631), and cytoskeletal protein-binding proteins (GO:0008092) (hypergeometric test, false discovery rate (FDR) < 0.05; Fig. 2e and Additional file 3F), suggesting specialized functions for NOVA–RNA regulation in MNs.
Since transcript abundance directly affects CLIP signal intensities, we expected to see higher MN/WSC CLIP signals on MN-enriched transcripts. Therefore, we tested whether transcripts known to be enriched in MN, when compared with WSC, showed enrichment of NOVA binding in MN NOVA CLIP. The numbers of WSC or MN CLIP reads within respective peaks were summed for each gene, followed by analysis using the Bioconductor edgeR package  to identify genes with differential enrichment of CLIP reads in MN or WSC; 474 and 2434 transcripts showed enriched and depleted NOVA binding in MN compared to WSC, respectively (FDR ≤ 0.01; Fig. 2f and Additional file 3G, H). Consistent with the cell type specificity of GFP-NOVA2 CLIP, the known MN markers Chat and Chodl were among the top transcripts with the most significantly enriched GFP-NOVA2 CLIP reads coverage (FDR = 2.61 × 10− 23 and 8.39 × 10− 44, respectively; Fig. 2f) . Conversely, Slc6a5, a glycine transporter gene known as an inhibitory neuron marker , as well as Erbb4, an interneuron-restricted receptor tyrosine kinase [32, 33], were two transcripts with the most significantly depleted MN CLIP reads coverage (FDR = 1.84 × 10− 37 and 2.31 × 10− 31, respectively, Fig. 2f).
We further systematically examined transcripts with well-defined anatomic expression patterns in spinal cord. Spinal cord gray matter exhibits a pattern of lamination consisting of ten laminal layers, with large alpha-MN pools located in lamina IX . Based on the Allen Spinal Cord Atlas generated and curated in situ hybridization data , 166 transcripts are exclusively expressed in lamina IX, of which 55 displayed differential NOVA binding (FDR ≤ 0.01) in MN. Meanwhile, 301 transcripts are restricted in one or more laminae other than lamina IX, of which 104 showed significant differential NOVA binding (FDR ≤ 0.01). Among this group of 159 (55 + 104) transcripts with differential NOVA binding in MN, the enrichment or depletion of NOVA CLIP signals on these transcripts in MN compared to WSC were highly concordant with transcript spatial expression patterns (green and purple dots in Fig. 2f), in that 42 of the 55 (76%) lamina IX enriched transcripts had enriched NOVA binding in MN (p < 2.2 × 10− 16, chi-squared test), and 97 out of the 104 (93%) lamina I–VIII enriched transcripts had lower NOVA binding in MN (p = 0.007, chi-squared test). Taken together, we conclude that HITS-CLIP in Chat:GFP-Nova2 lines captured MN-specific transcriptome-wide NOVA–RNA interactions.
NOVA displays MN-specific binding patterns
We were particularly interested in NOVA binding sites along a given gene that were disproportionally over- or under-represented in MNs compared to WSC (Additional file 2: Figure S3). To identify such NOVA–RNA interactions, we bioinformatically pooled WSC and MN CLIP reads to define peak regions referred as “joint peaks” (JPs), and grouped these JPs by genes (Additional file 2: Figure S3 and “Methods”). For each gene with two or more JPs, we performed pairwise Fisher’s exact test to identify JPs with disproportionally enriched MN or WSC CLIP reads (Additional file 2: Figure S3). This subset of JPs denotes NOVA binding sites over- or under-represented in MN compared to WSC.
Using this method, we identified 121,216 genic JPs after performing our peak-finding algorithm on pooled WSC and MN CLIP reads and filtering for biological reproducibility (Additional file 2: Figure S3). Analysis of intronic and exonic JPs as two separate groups revealed 707 (1.07%) over- and 160 (0.24%) under-represented intronic JPs in MNs (Fig. 3a and Additional file 4A), as well as 1058 (2.34%) over- and 407 (0.90%) under-represented exonic JPs in MNs (FDR ≤ 0.1, fold change ≥ 2; Fig. 3b and Additional file 4B). The AcGFP epitope tag does not significantly account for the observed NOVA binding differences between MN and WSC (Additional file 2: Figure S4A–C and Fig. 3a, b). Examples of disproportionally over- and under-represented NOVA binding sites are shown in Fig. 3c (intronic) and 3d (exonic).
We examined potential RNA sequence signatures around MN-specific NOVA binding sites (FDR ≤ 0.1, fold change ≥ 2) by analyzing neighborhood tetramer distributions. We discovered a marked overrepresentation of polypyrimidine (YYYY) tetramers around intronic NOVA binding sites over-represented in MN (Fig. 3e), as well as a striking enrichment of U-rich (two or more uridines) tetramers around 3′-UTR NOVA binding sites over-represented in MN (Fig. 3f). For both polypyrimidine and U-rich tetramers, differential enrichment is confined to regions 100 nt up- and downstream of over-represented MN peaks (Fig. 3e, f). Interestingly, regions flanking changed MN peaks are evolutionarily more conserved compared to all intronic or 3′-UTR NOVA peaks (Fig. 3e, f), suggesting the potential functional importance of differential NOVA binding sites in MN. Taken together, these data suggest that a large number of binding sites are differentially bound by NOVA independent of potential transcript level variations between WSC and MN, such that NOVA displays MN-specific binding patterns.
MN-specific NOVA binding predicts MN-specific alternative splicing
Since NOVA plays an important role in regulating neuronal transcript splicing, we tested whether MN-specific NOVA binding correlates with MN-specific alternative splicing. We took advantage of a high quality MN RNA-seq dataset, where spinal MNs in 3-month-old wild-type mice were collected by laser capture microdissection (LCM) and subjected to RNA-seq . For the WSC transcriptome, we performed RNA-seq on age-, genotype-, and strain background-matched whole spinal cords. Approximately 36 and 80 million mappable reads were obtained from either biological replicates of MN and WSC RNA-seq, respectively, and alternative splicing analysis was performed as described  (http://zhanglab.c2b2.columbia.edu/index.php/Quantas).
Around 30% (1620 out of 5418) of MN-expressed alternative exons showed biologically consistent differential splicing in MNs compared to WSC (FDR ≤ 0.1, BC = 2 out of 2; Fig. 4a and Additional file 5A). Strikingly, an even higher proportion of MN-expressed Nova targets are differentially spliced between MNs and WSC (281 out of 626, 45%, FDR ≤ 0.1, BC2 out of 2; Fig. 4a). The significant enrichment of Nova targets among exons differentially spliced in MNs (p < 2.2 × 10− 16, hypergeometric test) suggests that Nova contributes in an important way to MN-specific splicing patterns.
We tested whether MN-specific NOVA binding correlated with MN-specific alternative splicing. Although NOVA binding in broader regions may influence splice site choice, high confidence predictions on whether NOVA promotes or represses alternative exons is achieved when NOVA binds within a window of 400 nt around the regulated exons —upstream binding highly correlates with splicing repression, and downstream binding with splicing activation. Transcriptome wide, 13 NOVA binding sites over- or under-represented in MNs (FDR ≤ 0.1) are located in these regulatory “hotspots” around alternative exons with RNA-seq coverage sufficient for analysis. Interestingly, the majority of these exons (9 out of 13, 69%) showed MN-specific splicing patterns (FDR ≤ 0.1). Based on our NOVA–RNA map, differential NOVA binding patterns in MNs compared to WSC correctly predicted increase or decrease of alternative exon inclusion in MNs compared to WSC for eight out of the nine alternative exons (89%) (Additional file 5B; Fig. 4c, d; Additional file 2: Figure S5), consistent with NOVA mediating a direct action to regulate MN-specific alternative splicing.
Intriguingly, the seven genes hosting these nine MN-specific splicing events encode a functionally coherent set of proteins. Six of the seven gene products, i.e., MTSS1, MAP7D1, KIF21A, EPB4.1|3, CLASP1, GPHN, and KCNC3, are known to interact with cytoskeleton components (Fig. 4b) [37,38,39,40,41,42,43]. For example, MTSS1 binds to actin monomers and induces membrane protrusion [37, 44]. Mtss1 harbors two alternatively spliced exons, E12 and E12a, at the 3′ end of its coding sequence, and inclusion of E12a, but not E12, is necessary to promote neuritogenesis [45, 46]. Interestingly, MN CLIP revealed a unique NOVA binding site 109 nt upstream of Mtss1 E12 in MN (threefold increase in relative peak height; FDR = 0.046; Fig. 4c), which would predict based on the NOVA-RNA map that NOVA binding would inhibit E12 in MNs. Indeed, we observed a remarkably lower E12 inclusion rate in MNs compared to WSC (dI (MN-WSC) = − 0.69, FDR = 9.65 × 10− 7; Fig. 4c).
Another example is Kcnc3, which encodes a pan-neuronal voltage-gated potassium channel Kv3.3 [47, 48]. MN CLIP and RNA-seq revealed a new alternative exon E3a which showed a significantly higher inclusion rate in MNs compared to WSC (dI (MN-WSC) = 0.05, FDR = 1.38 × 10− 14). The increased utilization of E3a in MNs positively correlated with a dramatically enhanced NOVA binding site in MN 121 nt downstream of E3a (fourfold increase, FDR = 4.02 × 10− 7; Fig. 4d). Interestingly, inclusion of E3a would lead to a Kv3.3 isoform with an extended C-terminal proline-rich domain, which has been shown to modulate channel inactivation through triggering actin nucleation at the plasma membrane . Taken together, these observations suggest that unique NOVA binding patterns around alternative exons in MNs contribute to MN-specific biology, particularly in shaping the cytoskeleton and regulating cytoskeleton interactions in unique ways within spinal cord MNs.
NOVA differentially regulates Sept8 alternative last exon usage in MNs
Regulation of alternative last exon (ALE) usage involves intricate interplay between splicing and cleavage/polyadenylation, two processes both regulated by NOVA . We therefore investigated the potential role that NOVA played in regulating ALE usage in MNs. Sixty-five ALEs in 34 genes were differentially included in MNs compared to WSC (FDR ≤ 0.1, |dI| ≥ 0.2; Fig. 5a and Additional file 6A), with the top differentially utilized ALEs residing in two functionally related genes, Sept8 and Cdc42 (Fig. 5a) [49,50,51]. We also examined the published Nova2 wild-type and knockout (KO) mouse brain RNA-seq dataset , and identified additional NOVA regulation on 17 ALEs in 9 genes (FDR ≤ 0.1, |dI| ≥ 0.2, Fig. 5B and Additional file 6 B). Interestingly, among all genes with ALEs, Sept8 was the only NOVA target differentially regulated in MNs compared to WSC.
Sept8 encodes a family member of the septin proteins, which are multi-functional components of the cytoskeleton [52, 53]. In neurons, septins regulate dendritic and axon morphology through modulating actin and microtubule dynamics [54,55,56]. Although much is known about other septins, SEPT8 is a more recently described family member and less well characterized. Mouse Sept8 harbors two alternative terminal exons, exon 10a and 10b (Fig. 5c). Four mutually exclusive splice acceptors in Sept8 exons 10a and 10b can directly join downstream of exon 9, with splice acceptors 3 in exon 10a and 4 in exon 10b utilized in > 85% of Sept8 transcripts in adult mouse spinal cords (Fig. 5c). Comparison of ALE usage altered between Nova2 wild-type and KO mouse brains  showed markedly lower inclusion rate for exon 10b in Nova2 KO (dI = 0.42, FDR = 3.32 × 10− 9; Fig. 5c), suggesting that NOVA promotes splice acceptor 4 usage, exon 10a exclusion, and exon 10b inclusion (Fig. 5c). In MNs, RNA-Seq data indicate that the majority of Sept8 transcripts use splice acceptor 4 (65% in MNs vs 22% in WSC), while splice acceptor 3 is preferentially bypassed (20% in MNs vs 72% in WSC; Fig. 5c). This observed differential ALE usage between MNs and WSC coincides with higher NOVA binding in MNs at three binding sites located in exon 10a (site A fold change = 10, FDR = 0.0050; site B fold change = 3.4, FDR = 0.0046; site C fold change = 4.4, FDR = 3.84 × 10− 5; Fig. 5c). NOVA binding in site C is in a highly conserved YCAY-dense region encompassing the putative polyadenylation site at the 3′ end of exon 10a (Fig. 5d and Additional file 2: Figure S6A).
NOVA has been shown to prevent cleavage/polyadenylation through binding in close proximity to polyadenylation sites . We tested whether NOVA association with the 3′ end of Sept8 exon 10a blocked cleavage/polyadenylation, thus allowing for RNA polymerase II (PolII) readthrough and inclusion of the downstream exon 10b. A minigene was constructed using the genomic region of the 3′ end of exon 10a and the 5′ part of intron 10 (Fig. 5d). We co-transfected this minigene reporter with vectors expressing NOVA proteins or controls into COS-1 cells (Additional file 2: Figure S6B), which lack endogenous NOVA, followed by quantification of RNA levels up- and downstream of the E10a polyadenylation site. In the absence of NOVA proteins, transcription terminated efficiently at the E10a polyadenylation site, as measured by the less than 0.5% transcription readthrough rate (Fig. 5e). Exogenous NOVA proteins increased the downstream readthrough dramatically (> 10%, > 25-fold; Fig. 5e), while another RNA-binding protein (RBFOX3) showed no effect on cleavage/polyadenylation at E10a (Fig. 5e), indicating that NOVA proteins efficiently blocked cleavage/polyadenylation and promoted downstream transcription.
To test whether direct NOVA association with the YCAY sites is necessary for the observed regulation, we generated two mutant minigenes by disrupting YCAY sites while preserving the GC content (Fig. 5d). Mutant 1, with the two YCAY sites proximal to the E10a poly(A) site mutated, showed moderately dampened responses (6–11-fold) to NOVA overexpression (Fig. 5e). In contrast, little NOVA regulation (less than twofold; Fig. 5e) was observed for mutant 2, where all YCAY sites within 150 nt of the poly(A) site were disrupted. Taken together, these results suggest that NOVA promotes Sept8 exon 10b inclusion by binding close to the polyadenylation site in exon 10a and boosting readthrough transcription and utilization of exon 10b, and that strengthened NOVA binding around the exon 10a poly(A) site in MNs led to higher exon 10b usage observed in MNs.
The Sept8 exon specifically promoted by NOVA in MNs encodes a palmitoylated filopodia-inducing motif that enhances dendritic arborization
Alternative usage of exons 10a and 10b confers the C-terminal variation between SEPT8 protein isoforms X5 and X1, respectively. SEPT8 undergoes palmitoylation in vivo , which is a reversible post-translational process of attaching a 16-carbon saturated fatty acid to cysteine residues . Although the definitive palmitoylation site(s) in SEPT8 is unknown, the only predicted sites are two exon 10b-encoded cysteine residues (C469, C470) in the NOVA-promoted SEPT8-X1 isoform . To assess whether C469 and C470 mediate SEPT8-X1 palmitoylation, we expressed SEPT8 isoforms in COS-1 cells and employed the acyl-resin-assisted-capture (acyl-RAC) assay for palmitoylation detection [60, 61]. This assay relies on the indirect capture of palmitoylated proteins facilitated by palmitoyl-ester-specific cleavage (Additional file 2: Figure S7A). Whereas acyl-RAC failed to detect any palmitoylated SEPT8-X5, 5–10% of SEPT8-X1 was shown to be palmitoylated, as evidenced by the presence of SEPT8-X1 in the resin captured fraction (Fig. 6a, lane 2). This capture was dependent on palmitoyl-ester-specific cleavage (Fig. 6a, lanes 3 and 4). Mutations in C469 and C470 (SEPT8-X1-mut, C469S/C470S) prevented the detection of SEPT8-X1 palmitoylation (Fig. 6a), indicating that the two cysteines encoded by the NOVA-promoted exon 10b are the sites for palmitoyl addition. The data here collectively suggest that NOVA-regulated alternative RNA processing in MNs mediates isoform-specific SEPT8 palmitoylation.
Interestingly, we discovered that Sept8 exon 10b encodes a potential filopodia-inducing motif (FIM), characterized by two adjacent palmitoylated cysteines (C469, C470) and nearby basic residues (R475, R480) (Fig. 6a, orange box) . It has been shown that the palmitoylated FIM was sufficient to promote dendritic branching and spine formation in neurons . To assay the functional effects of SEPT8-X1 on dendritic morphology, we performed isoform-specific knockdown (KD) using a construct co-expressing GFP and an shRNA targeting Sept8 exon 10b (shX1) in primary mouse hippocampal neurons which express SEPT8-X1 (Additional file 2: Figure S7B). KD efficiency (79%) and isoform specificity of shX1 was demonstrated in COS-1 cells expressing exogenous SEPT8-X1 and X5 (Additional file 2: Figure S7C). Compared to neurons expressing the control scramble shRNA, neurons transfected with shX1 had significantly less complex dendritic arbors, as indicated by an 82% reduction in the number of dendritic branching points (Fig. 6b, e; p = 2.9 × 10− 9). Meanwhile, total dendritic length showed a 65% reduction (p = 9.2 × 10− 11), and sholl analysis revealed a similar 66% reduction in critical value in neurons with SEPT8-X1 KD (Fig. 6e; p = 4.4 × 10− 10). Even more strikingly, SEPT8-X1 KD led to an almost complete absence of dendritic spines—dendrites became smooth with a few focal swellings at distal ends (Fig. 6b, f; p =7.0 × 10− 82). These abnormalities were partially rescued by co-expression of the shRNA-resistant SEPT8-X1, but not SEPT8-X5 or the palmitoylation-deficient SEPT8-X1-mut (Fig. 6c, e, f). Therefore, we conclude that SEPT8-X1 promotes dendritic branching and spine formation through its palmitoylated FIM.
Since NOVA promotes the SEPT8-X1 isoform, we predicted that neurons lacking NOVA would exhibit a similar reduction in dendrite arbor and spine density. We knocked down NOVA2 using an shRNA (shNova2) which efficiently depleted 76% of the endogenous NOVA2 in N2A cells (Additional file 2: Figure S7D). Similar to SEPT8-X1 KD, neurons depleted of NOVA2 displayed greatly decreased dendritic arbors compared to control shRNA transfected neurons, as evidenced by a 55% reduction in the number of dendritic branching points, a 47% reduction in total dendritic length, and a 48% reduction in the sholl analysis critical value (Fig. 6b, e; p = 3.6 × 10− 7, 2.5 × 10− 7, 1.8 × 10− 7, respectively). On the other hand, to our surprise, NOVA2 knockdown did not significantly affect dendritic spine density (Fig. 6b, f). Co-expression of SEPT8-X1, but not the two non-palmitoylated variants, partially restored dendritic arbor complexity in neurons with NOVA2 KD (Fig. 6d, e), suggesting that NOVA2 enhances dendritic arborization through promoting the FIM-containing SEPT8 isoform.
Understanding brain function involves understanding its parts, as demonstrated by the discovery of differences in ribosome-associated transcripts by looking at specific neuronal cell types in the basal ganglia (D1 vs D2 neurons) . Here we develop a general method combining BAC-transgenic engineering and CLIP that is conceptually applicable to the study of any protein–RNA interactions within specific cell types. We apply this strategy to analyze differential NOVA binding in mouse spinal MNs, compare that to interactions visible at the gross level of whole spinal cord analysis, and uncover NOVA-regulated biology specific to MNs. MN-specific CLIP revealed a major role for NOVA in regulating cytoskeleton interactions in MNs, a function obscured in previous whole tissue-based analyses. This led us to discover a consequent defect in dendritic morphology in neurons lacking NOVA, which may help explain in part the severe MN defect seen in NOVA1/2 double KO MNs .
MN-specific CLIP unmasks previously undetected aspects of NOVA function
MNs are among the largest neurons in the central nervous system. Their distinct morphology, characterized by a long axon and intricate dendritic arbor, renders traditional cell purifications based on enzymatic tissue digestion and cell purification methods particularly unsatisfactory for understanding the molecular biology of the whole neuron, especially given abundant evidence for RNA localization within the dendritic arbor. The cell type-specific CLIP strategy developed here allows robust and quantitative identification of NOVA binding sites in MNs in vivo, in both the cell bodies and processes. This generated a transcriptome-wide NOVA–RNA interaction atlas in MNs with NOVA binding sites in over 4000 genes. NOVA binding in MNs reflects MN transcriptome signatures, which further demonstrated the cell type specificity of our assay.
The cell type-specific CLIP strategy developed here allowed delineation of a subset of MN NOVA targets from WSC. These MN NOVA targets are enriched in genes encoding synaptic functions to a similar extent compared to WSC NOVA targets, yet they are especially enriched in genes encoding cytoskeleton interacting proteins. Enrichment of MN versus WSC CLIP signals in this subset of genes is significant even when adjusted to MN/WSC transcript level differences (p < 2.2 × 10− 16; Additional file 2: Figure S8), suggesting that these gene ontology (GO) enriched NOVA target genes contain quantitatively strengthened NOVA binding sites in MNs compared to WSC. Furthermore, through combining cell type-specific CLIP with MN transcriptome profiling, we discovered RNA processing events differentially regulated by NOVA in MNs. Around 2% of NOVA binding sites were over- or under-represented in MNs, which is consistent with our findings that ~ 2% (eight alternatively splicing regions out of 303) of NOVA-regulated alternative splicing showed MN-specific splicing patterns correctly predicted by MN-specific NOVA binding. Interestingly, the vast majority of these events also reside in transcripts encoding cytoskeleton-interacting proteins. These findings reveal a previously undiscovered cell type-specific role for NOVA in regulating MN cell biology.
NOVA plays an important role in regulating MN cytoskeleton interactions
The top NOVA-regulated ALE was in Sept8, which encodes a member of the multifunctional septin family that is capable of regulating neurite outgrowth and branching through interactions with cytoskeleton components [54,55,56]. NOVA binding at the polyadenylation site in Sept8 exon 10a correlated with exon 10b usage, and this was only evident in analysis of MN CLIP, not WSC CLIP alone. Consistent with prior observations of position-dependent effects of NOVA on APA , we demonstrated that NOVA directly inhibits cleavage/polyadenylation by binding close to the exon 10a poly(A) site, presumably promoting exon 10b transcription and splicing. Interestingly, the NOVA-dependent exon 10b encodes an FIM capable of promoting dendritic branching and spine formation. Indeed, we discovered that SEPT8-X1, the NOVA-dependent SEPT8 isoform harboring the FIM, specifically promotes dendritic arborization and spine formation. Moreover, we found that neurons lacking NOVA2 displayed decreased dendritic arbors, which are partially rescued by SEPT8-X1. Taken together, these data indicate that NOVA promotes dendritic arborization through Sept8 regulation.
While SEPT8-X1 promotes dendritic spine formation, we were not able to detect changes in dendritic spine density in NOVA2 KD neurons. This may be related to insufficient effects from NOVA2 KD, or from contributing indirect factors, such as NOVA promotion of protein isoforms with antagonistic effects on spine formation. For example, based on in silico palmitoylation site prediction , we identified a total of eight NOVA targets where a predicted FIM is regulated by NOVA. Of these eight genes, NOVA promotes the FIM-harboring isoform in four (Sept8, Ccp110, Sdccag3, and Ccdc84), while inhibiting the FIM-encoding exon in the other four (Sgce, Clip1, Kcnma1, and Ube2e2). It is possible that changes of various NOVA targets upon NOVA KD mitigate the overall effect on dendritic spine density.
We have previously shown that NOVA proteins play an essential role in MN physiology . MNs in mice lacking both Nova family members were paralyzed and failed to cluster acetyl-choline receptors at the neuromuscular junctions . Successful rescue of the neuromuscular junction defect was evident after Nova-knockout MNs were engineered to constitutively express the Nova-regulated Z+ alternatively spliced isoform of agrin, but surprisingly remained paralyzed. Phrenic nerve stimulation revealed that the axonal synaptic machinery for conducting action potentials and synaptic vesicle release from the axon was functional, leading to the conclusion that dysregulation of additional MN NOVA RNA targets contribute to a proximal physiologic defect in Nova-KO MNs. Here, MN-specific CLIP reveals a major role of NOVA in regulating the MN cytoskeleton, including promoting dendritic complexity. Dendrites are the main information-receiving sites of neurons. Spinal MNs, as the gateway controllers of the CNS motor outputs, have elaborate dendritic structures to meet the highly complex demand of precisely coordinating muscle contractions spatially and temporally [63, 64]. Early and progressive dendritic degeneration has been reported in lower MNs in MN disease mouse models as well as amyotrophic lateral sclerosis (ALS) patients [65, 66], suggesting the importance of dendritic integrity to MN function. Our new findings may provide additional avenues for understanding the role that NOVA and RNA regulation plays in MN function.
Differential NOVA binding sites in MN suggests combinatorial control of multiple RNA-binding proteins
When we compared the NOVA–RNA interactions in MNs and WSC, we uncovered over 2000 sites transcriptome-wide that are differentially bound by NOVA in MNs compared to WSC. In this way cell-specific CLIP offers the possibility of revealing cell type-specific regulatory phenomena that are otherwise obscured from analysis of whole tissues. What underlies this unique NOVA binding profile in MNs? It has been shown that combinatorial control is integral in RNA processing regulation. Cellular RNA-binding protein networks play a pivotal role in determining target selection and binding dynamics of a given RNA-binding protein . For NOVA, interactions with PTBP2, the neuronal polypyrimidine tract binding protein, as well as RBFOX proteins have previously been demonstrated [20, 68]. In particular, we have shown that PTBP2 interacts with and antagonizes NOVA in the regulation of glycine receptor alpha 2 subunit (Glra2) splicing . Interestingly, compared to all NOVA binding sites, sequences around NOVA bindings sites that are over-represented in MNs showed enrichment of pyrimidine-rich motifs and higher evolutionary conservation (Fig. 3e, f). This sequence signature around MN-specific NOVA binding sites along with lower PTBP2 levels in MNs  suggests the hypothesis that lower abundance of PTBP2 in MNs allows for unmasking of certain NOVA binding sites that are otherwise occupied by PTBP2 in other neurons. This kind of intricate cell type-specific interaction network may help determine Nova target selection beyond the specificity determined by sequence and structural constraints.
Cell type-specific CLIP
Cell type-specific CLIP is in general applicable to any cell type and a great variety of RNA-binding proteins. The current study is neuron-specific due to the nature of Nova [25, 69], but a further array of cell types within the CNS or other tissues can be studied in vivo through epitope tagging. Indeed, Schaefer and colleagues used a similar strategy to identify AGO-bound miRNAs in D2 neurons of the mouse striatum , and Camk2a-neurons in the forebrain . And as CLIP is applied to identify functional protein–RNA interactions defining the actions of an increasing number of RNA-binding proteins, including splicing factors , RNA-editing factors , and epigenetic regulators , expanding aspects of post-transcriptional and RNA-related regulation can be investigated in a cell type-specific manner. We recently described a method, PAPERCLIP, in which PABPC1 CLIP was used to identify brain-specific transcripts and alternative 3′ UTR processing within those cells . Combining a cell-specific tagging strategy with PAPERCLIP would allow profiling of cell-specific polyadenylated transcripts in the brain or within any tissue, an approach that would complement the delineation of ribosome-associated transcripts demarcated by BAC-TRAP.
We examined whether GFP-NOVA2 expression in our BAC-transgenic MNs could have affected our observations relative to unperturbed neurons. The BAC-transgenic system has intrinsic expression biases and might, for example, have skewed the MN NOVA-PTBP ratios, impacting the observed NOVA binding increases around pyrimidine-rich regions. However, the high concordance between NOVA binding changes in our transgenic MNs and differential splicing in MNs without exogenous NOVA strongly argues for the physiological relevance of our transgenic model. A more elegant approach of studying cell type-specific RNA–protein interactions would be to generate conditionally epitope-tagged knock-in lines followed by the introduction of a cell type-specific Cre recombinase expression. Epitope-tagged protein would be expressed in the desired cell type from its endogenous promoter, thus preserving protein stoichiometry; this strategy has already shown promise with the ubiquitous RNA-binding protein PABPC1 .
Caution should be taken when comparing cell type-specific RNA binding profiles across different cell types or tissues. One issue is that differences in CLIP library complexities between comparison groups may lead to unequal statistical power in identifying over- versus under-represented binding sites. In our case, although we pooled spinal cords from multiple transgenic mice per biological replicate and performed twice as many biological replicates as the WSC CLIP, we were only able to generate less than one-eighth of CLIP reads from MN-specific CLIP compared to WSC CLIP. This is likely due to low abundance of MN-specific GFP-NOVA2 compared to the WSC NOVA pool (Fig. 1b). We hypothesize that the unevenness of library complexities contributed to the identification of a high percentage (> 72%) of NOVA binding sites significantly over-represented in MNs compared to WSC (Fig. 3a, b). We predict that achieving higher MN CLIP library complexity by including more BC and/or pooling more samples per BC can alleviate this bias.
Another challenge of analyzing cell type-specific CLIP data involves obtaining the matching cell type-specific transcriptomes, which is essential for teasing out functionally important protein–RNA interactions. A variety of technologies, such as fluorescence-activated cell sorting, LCM, and TRAP-seq can be explored for this purpose. We utilized a published RNA-seq dataset on laser captured MNs for the MN-specific transcriptome. Our transcriptome reference for WSC NOVA CLIP, on the other hand, is less than ideal. Since NOVA expression is restricted in neurons, the ideal dataset would be neuron-specific transcriptome profiling generated using technologies such as TRAP-seq or cTag PAPERCLIP.
Cell type-specific CLIP may be relevant for the study of many neurological diseases, such as ALS and ataxias, where defined cell types are pathologically affected during the entire course or the initial stages of disease progression. Interplay between the affected cells and their cellular context plays an important role in disease progression [76,77,78,79]. Our strategy allows for in vivo dissection of both cell autonomous and non-autonomous effects resulting from disease-causing mutations or insults, which had not been possible before. Great insights on cell type contributions to disease pathogenesis will be gained upon application of this strategy to a variety of animal models.
Here we demonstrate the feasibility and physiological relevance of delineating neuronal cell type-specific RNA regulation. Through cell type-specific epitope tagging of the RNA-binding protein NOVA2, we generated a MN-specific NOVA–RNA interaction map using CLIP. Cell type-specific CLIP revealed a major role of NOVA in regulating cytoskeleton interactions in MNs, including promoting the palmitoylated isoform of a cytoskeleton protein, Septin 8, which enhances dendritic arbor complexity. MN-specific NOVA binding predicts MN-specific alternative RNA processing, further supporting the idea that cell type-specific RNA regulation contributes to cell identity. Our study highlights a non-incremental gain of knowledge moving from whole tissue- to single cell type-based RNA regulation analysis. As our strategy for cell type-specific CLIP is highly adaptable, we envision wide application of this method in a variety CNS cell types as well as disease models.
Antibodies and dilutions
The following antibodies were used in CLIP experiments: anti-GFP mouse monoclonal antibodies 19F7 and 19C8 (Memorial Sloan Kettering Monoclonal Antibody Facility), anti-NOVA serum from a paraneoplastic opsoclonus-myoclonus ataxia (POMA) patient. The following antibodies and dilutions were used for immunoblotting: anti-NOVA patient serum (1:2000), mouse anti-GAPDH monoclonal antibody 6C5 (Abcam, 1:20,000), rabbit anti-RBFOX3 (1:500) , rabbit anti-HA monoclonal antibody C29F4 (Cell Signaling Technology, 1:100). The following antibodies and dilutions were used for immunofluorescence: rat anti-GFP monoclonal antibody (nacalai tesque, 1:1000), goat anti-GFP polyclonal antibody (Rockland, 1:500), goat anti-CHAT polyclonal antibody (EMD Millipore, 1:500), rabbit anti-HA monoclonal antibody C29F4 (Cell Signaling Technology, 1:1000), rabbit anti-SEPT8 monoclonal antibody (pan SEPT8, EPR16099, 1:200), mouse anti-SEPT8_X1 monoclonal antibody D-11 (Santa Cruz, 1:50), chicken anti-MAP2 antibody (Thermo Fisher Scientific, 1:1000).
Cell culture and transfection
NIH3T3 and COS-1 cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 10% heat-inactivated fetal bovine serum, 100 U/mL penicillin, and 100 μg/mL streptomycin. NIH3T3 and COS-1 cells were transfected with plasmid constructs using lipofectamine 2000 (Invitrogen), according to the manufacturer’s instructions.
Mouse hippocampal neurons were isolated from 18-day-old CD1 mouse embryos and cultured in 24-well plates according to established protocols. For shRNA knockdown and rescue experiments, 0.4 μg of shRNA vector and 0.6 μg of protein expressing vector or control were transfected at DIV 10 using 3 μL of NeuroMag following the manufacturer’s protocol. Transfected neurons were fixed at DIV 12 for immunofluorescence staining and confocal imaging (see below).
Construction of BAC-transgenic mouse lines expressing GFP-NOVA2 in MNs
In brief, AcGFP-fused Nova2 coding sequence was cloned into a pLD53.SC2 plasmid. Subsequently, a DNA fragment homologous to Chat 5′ UTR region was PCR amplified using primers GCCAGGCATCTGAGAGGC and CCTAGCGATTCTTAATCCAGAGTAGCAGAGCTG and inserted into the plasmid through blunt end ligation at AgeI site. The sequence of the resulting plasmid pLD53.SC2-AcGFP-Nova2 was confirmed by Sanger sequencing and deposited to the Addgene database. Recombinant BAC was generated using RP23-246B12 and pSC2-GFP-Nova2 as described previously  and microinjected into pronuclear oocytes of C57BL/6 mice from Charles River. By PCR genotyping, nine mice from 68 offspring were confirmed to carry the transgene. Two founders, #6 and #17, were bred with C57BL/6 mice from Charles River to establish stable transgenic lines.
Paraformaldehyde-fixed spinal cords were transversely sectioned at 14 μm thickness and kept at − 80 °C. Before immunofluorescence staining, sections were rehydrated in PBS for 10 min, followed by a 1-h block in blocking buffer containing 100 mM Tris-HCl, pH 7.5, 150 mM NaCl, 0.2% Triton X-100, 5% horse serum. Primary antibodies diluted in blocking buffer were subsequently added to the sections followed by overnight incubation at room temperature. Secondary antibody hybridization and TO-PRO-3 staining were performed using Alexa Fluor conjugated antibodies diluted in 100 mM Tris-HCl, pH 7.5, 150 mM NaCl, 0.2% Triton X-100 with 1 μM TO-PRO-3. Confocal images were taken using an inverted LSM 510 laser scanning confocal microscope (Zeiss).
For immunofluorescence staining on hippocampal neurons, cells were fixed in 4% paraformaldehyde at room temperature for 10 minutes, permeabilized with 0.1% Triton-X 100 in PBS at room temperature for 15 minutes, and blocked in PBS containing 0.1% Tween-20 and 10% donkey serum for 1 h at room temperature. The primary antibodies, diluted in PBS containing 1% BSA and 0.1% Tween-20, were added to the cells and incubated overnight at room temperature. Alexa Fluor conjugated secondary antibodies (Jackson Immunoresearch) were used at 1:500 dilution to detect either the tag epitope or proteins of interest. Confocal images were taken using an inverted LSM 880 NLO laser scanning confocal microscope (Zeiss).
HITS-CLIP experiments and analysis
Spinal cords were hydraulically extruded in Hank’s Balanced Salt Solution (HBSS) and triturated seven times with a 21½ gauge needle in 10 mL HBSS. Immediately after trituration, tissue suspensions were crosslinked 400 mJ/cm2 three times at 0 °C. Crosslinked tissues were collected by centrifugation at 4000 g for 2 minutes. HITS-CLIP experiments were performed as previously described [16, 26] with PCR primers specified in Additional file 7. For MN-specific CLIP, 400 μL of Dynabeads protein G (Invitrogen) precoated with 50 μg of each anti-GFP monoclonal antibody (mAb) was used to immunoprecipitate GFP-NOVA2 in UV crosslinked spinal cord lysate pooled from five to seven 3-month-old BAC-transgenic mice. For WSC CLIP, 400 μL of Dynabeads protein G precoated with 80 μL of human anti-NOVA serum was used to immunoprecipitate lysate from one 3-month-old wild-type mouse spinal cord. cDNA libraries were prepared as previously described and sequenced on an Illumina Genome Analyzer IIx, HiSeq 1000 or HiSeq 2000.
Bioinformatic processing and mapping of CLIP NGS reads were performed similarly as previously described with slight modifications [18, 80, 82, 83]. Specifically, we removed the 3′ linker sequence from CLIP reads before mapping the remaining sequences to the mm9 build of the mouse genome with novoalign (http://www.novocraft.com) without iterative trimming using parameters –t 85 -l 25 -s 0 -o Native -r None. Following mapping, pentamer frequencies immediately downstream of CLIP reads were ranked. CLIP reads immediately upstream of the top 20 pentamers with zero or one mismatch to GTGTC were removed. We have recently discovered that under our reverse transcription (RT) conditions, the 3′ end of the RT primer could prime reverse transcription by hybridizing to GUGUC or similar pentamers in the CLIPed RNA, leading to their preferential amplification (Park C., personal communication). Bioinformatically removing such reads circumvents this technical caveat.
NOVA binding peaks were defined as previously described using gene regions compiled from refseq, UCSC known genes, and ESTs plus their downstream 10 kb as transcription units [26, 27, 80]. Gene regions with Entrez IDs were referred to as known genes. Defined peaks were then filtered for biological complexity, requiring CLIP reads from at least half of the biological replicates. Joint peaks (JPs), in particular, were defined by running the peak finding algorithm using CLIP reads pooled from both comparison groups (e.g., MNs and WSC), and requiring CLIP reads from at least half of the biological replicates in either group (i.e., WSC BC 2 out of 4 or MN BC 4 out of 8).
Gene-wise CLIP read enrichment between WSC and MNs was evaluated using the Bioconductor edgeR package with tagwise dispersion model . P values of differential NOVA binding on specific sites were calculated using Fisher’s exact test on the following 2 × 2 matrix:
where m and n stand for the number of MN and WSC CLIP reads in a given JP, respectively, while M and N stand for the total number of intronic/exonic MN and WSC CLIP reads within all intronic/exonic JPs in the corresponding transcript, respectively. Coverage is defined as the smallest number among the 2 × 2 matrix. Benjamini-Hochberg multiple-testing correction was performed using all intronic or exonic JPs with a minimum coverage of 10.
RNA-seq library preparation and data analysis
We performed WSC RNA-seq on two biological replicates of 3-month-old SOD1WT transgenic mice in SJL/B6 mixed background. For each replicate, poly(A) selected RNA from spinal cord was used to prepare RNA-seq library using an Illumina TruSeq RNA sample prep kit. We generated 100-nt paired-end reads using an Illumina Genome Analyzer II. In order to compare WSC and MN RNA-seq, we trimmed our WSC reads to 74 nt to match the read length of the MN RNA-seq dataset by Bandyopadhyay et al. (GSE38820). Both WSC and MN datasets were mapped to the mouse (mm9) genome and exon junctions using OLego with default parameters (15-nt seed with 1 nt overlapping, four or fewer mismatches per read) . Only reads unambiguously mapped to the genome or exon junctions were retained for downstream analysis. For MN RNA-seq, 36,484,398 and 38,293,208 NGS reads were mapped for the two biological replicates, respectively, while 74,444,847 and 81,175,598 NGS reads were mapped for WSC RNA-seq replicates. Transcript level and alternative splicing analyses were performed as previously described. Bioconductor edgeR package was used to evaluate statistical significance of transcript level differences between WSC and MNs . Fisher’s exact test was used to calculate p values of differential alternative splicing, and the FDR was estimated using the Benjamini-Hochberg method. Alternative exon splicing was analyzed as described . Alternative exons with splice junction read coverage over 10  were considered “expressed”, and were included in Benjamini-Hochberg multiple test correction. Differential alternative splicing events were identified by requiring FDR ≤ 0.1 in addition to biological consistency (BC2 out of 2).
Gene ontology analysis
GO analysis was performed using GOrilla by running unranked lists of target and background genes . For background genes, we used all genes with rpkm ≥ 1 in WSC or MNs.
Motif enrichment analysis
For motif enrichment around differential MN intronic NOVA peaks, tetramer occurrences 100 nt around the center of each peak were counted and compared to those 100 nt around all intronic or exonic NOVA peaks using a hypergeometric test.
Sept8 minigene assay
To construct the Sept8 minigene, the C57BL/6 genomic region amplified with primers TACGACTCACTATAGGGCGAATTCGGATCCGCATGAATTCTGACCCCTGTGA and CAATAAACAAGTTCTGCTTTAATAAGATCTCCGTAACCTGGCTACCAGTGA was cloned into BamHI and BglII digested pSG5 vector using HIFI DNA assembly kit (New England Biolabs). Mutant minigenes were constructed by assembling the PCR amplified vector backbone (forward, TTGGGCAGTAGCTTCGCTG; reverse, TAACATAACAGAGAAGCAAGCTGGCT) with synthetic gBlocks (IDT DNA) carrying the desired mutations. Minigene (1 μg) was co-transfected into COS-1 cells with 1 μg of control vector or NOVA/RBFOX3 expression vector. Cells were harvested 48 h post-transfection for immunoblotting or RT-qPCR following standard protocols. Primers CTGACCCCGCATATGTTCCTGTGT and CAAGCTGGCTATCCTGGGCCTCTT were used to amplify an exon 10a region, while primers GTCTGCGATGGTTTTGCAGAGGTG and GGCCACAGGAAATGGAGATGTGAG were used for an intron 10 amplicon. Three independent experiments were performed for statistical comparison.
We seeded 250,000 COS-1 cells per well in six-well plates and transfected them with 2 μg of constructs expressing FLAG-HA-tagged SEPT8 variants 18–24 h later. Acyl-RAC assay was performed using the CAPUREome S-Palmitoylated Protein Kit (Badrilla) according to the manufacturer’s protocol. Treated protein lysates were subsequently immunoblotted with rabbit anti-HA antibody.
Maximum projected confocal images were produced using ImageJ using only linear adjustments. Dendrites were semi-automatically traced using the Simple Neurite Tracer plugin. Critical value and number of dendritic branching points were deduced by using the “sholl analysis” and “analyze skeleton” plugins. Fifteen to 17 neurons were analyzed per group.
Nelson SB, Sugino K, Hempel CM. The problem of neuronal cell types: a physiological genomics approach. Trends Neurosci. 2006;29:339–45.
Hobert O, Carrera I, Stefanakis N. The molecular and gene regulatory signature of a neuron. Trends Neurosci. 2010;33:435–45.
Fishell G, Heintz N. The neuron identity problem: form meets function. Neuron. 2013;80:602–12.
Darnell RB. RNA protein interaction in neurons. Annu Rev Neurosci. 2013;36:243–70.
Buckanovich RJ, Posner JB, Darnell RB. Nova, the paraneoplastic Ri antigen, is homologous to an RNA-binding protein and is specifically expressed in the developing motor system. Neuron. 1993;11:657–72.
Taliaferro JM, Vidaki M, Oliveira R, Olson S, Zhan L, Saxena T, et al. Distal alternative last exons localize mRNAs to neural projections. Mol Cell. 2016;61:821–33.
Zhang X, Chen MH, Wu X, Kodani A, Fan J, Doan R, et al. Cell-type-specific alternative splicing governs cell fate in the developing cerebral cortex. Cell. 2016;166:1147–62. e15
Zhang Y, Chen K, Sloan SA, Bennett ML, Scholze AR, O'Keeffe S, et al. An RNA-sequencing transcriptome and splicing database of glia, neurons, and vascular cells of the cerebral cortex. J Neurosci. 2014;34:11929–47.
Cahoy JD, Emery B, Kaushal A, Foo LC, Zamanian JL, Christopherson KS, et al. A transcriptome database for astrocytes, neurons, and oligodendrocytes: a new resource for understanding brain development and function. J Neurosci. 2008;28:264–78.
Sanz E, Yang L, Su T, Morris DR, McKnight GS, Amieux PS. Cell-type-specific isolation of ribosome-associated mRNA from complex tissues. Proc Natl Acad Sci U S A. 2009;106:13939–44.
Heiman M, Kulicke R, Fenster RJ, Greengard P, Heintz N. Cell type-specific mRNA purification by translating ribosome affinity purification (TRAP). Nat Protoc. 2014;9:1282–91.
Heiman M, Schaefer A, Gong S, Peterson JD, Day M, Ramsey KE, et al. A translational profiling approach for the molecular characterization of CNS cell types. Cell. 2008;135:738–48.
Doyle JP, Dougherty JD, Heiman M, Schmidt EF, Stevens TR, Ma G, et al. Application of a translational profiling approach for the comparative analysis of CNS cell types. Cell. 2008;135:749–62.
Jensen KB, Musunuru K, Lewis HA, Burley SK, Darnell RB. The tetranucleotide UCAY directs the specific recognition of RNA by the Nova K-homology 3 domain. Proc Natl Acad Sci U S A. 2000;97:5740–5.
Ule J, Ule A, Spencer J, Williams A, Hu J-S, Cline M, et al. Nova regulates brain-specific splicing to shape the synapse. Nat Genet. 2005;37:844–52.
Licatalosi DD, Mele A, Fak JJ, Ule J, Kayikci M, Chi SW, et al. HITS-CLIP yields genome-wide insights into brain alternative RNA processing. Nature. 2008;456:464–9.
Ule J, Stefani G, Mele A, Ruggiu M, Wang X, Taneri B, et al. An RNA map predicting Nova-dependent splicing regulation. Nature. 2006;444:580–6.
Moore MJ, Zhang C, Gantman EC, Mele A, Darnell JC, Darnell RB. Mapping Argonaute and conventional RNA-binding protein interactions with RNA at single-nucleotide resolution using HITS-CLIP and CIMS analysis. Nat Protoc. 2014;9:263–93.
Darnell RB. HITS-CLIP: panoramic views of protein-RNA regulation in living cells. WIREs RNA. 2010;1:266–86.
Zhang C, Frias MA, Mele A, Ruggiu M, Eom T, Marney CB, et al. Integrative modeling defines the Nova splicing-regulatory network and its combinatorial controls. Science. 2010;329:439–43.
Ruggiu M, Herbst R, Kim N, Jevsek M, Fak JJ, Mann MA, et al. Rescuing Z+ agrin splicing in Nova null mice restores synapse formation and unmasks a physiologic defect in motor neuron firing. Proc Natl Acad Sci U S A. 2009;106:3513–8.
Huang CS, Shi S-H, Ule J, Ruggiu M, Barker LA, Darnell RB, et al. Common molecular pathways mediate long-term potentiation of synaptic excitation and slow synaptic inhibition. Cell. 2005;123:105–18.
Yano M, Hayakawa-Yano Y, Mele A, Darnell RB. Nova2 regulates neuronal migration through an RNA switch in disabled-1 signaling. Neuron. 2010;66:848–58.
Saito Y, Miranda-Rottmann S, Ruggiu M, Park CY, Fak JJ, Zhong R, et al. NOVA2-mediated RNA regulation is required for axonal pathfinding during development. elife. 2016;5:e14371.
Yang YY, Yin GL, Darnell RB. The neuronal RNA-binding protein Nova-2 is implicated as the autoantigen targeted in POMA patients with dementia. Proc Natl Acad Sci U S A. 1998;95:13254–9.
Charizanis K, Lee KY, Batra R, Goodwin M, Zhang C, Yuan Y, et al. Muscleblind-like 2-mediated alternative splicing in the developing brain and dysregulation in myotonic dystrophy. Neuron. 2012;75:437–50.
Shah A, Qian Y, Weyn-Vanhentenryck SM, Zhang C. CLIP Tool Kit (CTK): a flexible and robust pipeline to analyze CLIP sequencing data. Bioinformatics. 2017;33:566–7.
Gong S, Doughty M, Harbaugh CR, Cummins A, Hatten ME, Heintz N, et al. Targeting Cre recombinase to specific neuron populations with bacterial artificial chromosome constructs. J Neurosci. 2007;27:9817–23.
Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.
Enjin A, Rabe N, Nakanishi ST, Vallstedt A, Gezelius H, Memic F, et al. Identification of novel spinal cholinergic genetic subtypes disclose Chodl and Pitx2 as markers for fast motor neurons and partition cells. J Comp Neurol. 2010;518:2284–304.
Kodama T, Guerrero S, Shin M, Moghadam S, Faulstich M, Lac du S. Neuronal classification and marker gene identification via single-cell expression profiling of brainstem vestibular neurons subserving cerebellar learning. J Neurosci. 2012;32:7819–31.
Vullhorst D, Neddens J, Karavanova I, Tricoire L, Petralia RS, McBain CJ, et al. Selective expression of ErbB4 in interneurons, but not pyramidal cells, of the rodent hippocampus. J Neurosci. 2009;29:12255–64.
Neddens J, Fish KN, Tricoire L, Vullhorst D, Shamir A, Chung W, et al. Conserved interneuron-specific ErbB4 expression in frontal cortex of rodents, monkeys, and humans: implications for schizophrenia. Biol Psychiatry. 2011;70:636–45.
Rexed B. The cytoarchitectonic organization of the spinal cord in the cat. J Comp Neurol. 1952;96:414–95.
Henry AM, Hohmann JG. High-resolution gene expression atlases for adult and developing mouse brain and spinal cord. Mamm Genome. 2012;23:539–49.
Bandyopadhyay U, Cotney J, Nagy M, Oh S, Leng J, Mahajan M, et al. RNA-Seq profiling of spinal cord motor neurons from a presymptomatic SOD1 ALS mouse. PLoS One. 2013;8:e53575.
Saarikangas J, Kourdougli N, Senju Y, Chazal G, Segerstråle M, Minkeviciene R, et al. MIM-induced membrane bending promotes dendritic spine initiation. Dev Cell. 2015;33:644–59.
Yu N, Signorile L, Basu S, Ottema S, Lebbink JHG, Leslie K, et al. Isolation of functional tubulin dimers and of tubulin-associated proteins from mammalian cells. Curr Biol. 2016;26:1728–36.
Lee K-H, Lee JS, Lee D, Seog D-H, Lytton J, Ho W-K, et al. KIF21A-mediated axonal transport and selective endocytosis underlie the polarized targeting of NCKX2. J Neurosci. 2012;32:4102–17.
Baines AJ, Lu H-C, Bennett PM. The protein 4.1 family: hub proteins in animals for organizing membrane proteins. Biochim Biophys Acta. 2014;1838:605–19.
Al-Bassam J, Chang F. Regulation of microtubule dynamics by TOG-domain proteins XMAP215/Dis1 and CLASP. Trends Cell Biol. 2011;21:604–14.
Galjart N. CLIPs and CLASPs and cellular dynamics. Nat Rev Mol Cell Biol. 2005;6:487–98.
Zhang Y, Zhang X-F, Fleming MR, Amiri A, El-Hassar L, Surguchev AA, et al. Kv3.3 channels bind Hax-1 and Arp2/3 to assemble a stable local actin network that regulates channel gating. Cell. 2016;165:434–48.
Mattila PK, Salminen M, Yamashiro T, Lappalainen P. Mouse MIM, a tissue-specific regulator of cytoskeletal dynamics, interacts with ATP-actin monomers through its C-terminal WH2 domain. J Biol Chem. 2003;278:8452–9.
Glassmann A, Molly S, Surchev L, Nazwar TA, Holst M, Hartmann W, et al. Developmental expression and differentiation-related neuron-specific splicing of metastasis suppressor 1 (Mtss1) in normal and transformed cerebellar cells. BMC Dev Biol. 2007;7:111.
Sistig T, Lang F, Wrobel S, Baader SL, Schilling K, Eiberger B. Mtss1 promotes maturation and maintenance of cerebellar neurons via splice variant-specific effects. Brain Struct Funct. 2017;222(6):2787–2805.
Goldman-Wohl DS, Chan E, Baird D, Heintz N. Kv3.3b: a novel Shaw type potassium channel expressed in terminally differentiated cerebellar Purkinje cells and deep cerebellar nuclei. J. Neurosci. 1994;14:511–22.
Brooke RE, Atkinson L, Edwards I, Parson SH, Deuchars J. Immunohistochemical localisation of the voltage gated potassium ion channel subunit Kv3.3 in the rat medulla oblongata and thoracic spinal cord. Brain Res. 2006;1070:101–15.
Gladfelter AS, Bose I, Zyla TR, Bardes ESG, Lew DJ. Septin ring assembly involves cycles of GTP loading and hydrolysis by Cdc42p. J Cell Biol. 2002;156:315–26.
Sadian Y, Gatsogiannis C, Patasi C, Hofnagel O, Goody RS, Farkasovský M, et al. The role of Cdc42 and Gic1 in the regulation of septin filament formation and dissociation. elife. 2013;2:e01085.
Kinoshita M. Assembly of mammalian septins. J Biochem. 2003;134:491–6.
Weirich CS, Erzberger JP, Barral Y. The septin family of GTPases: architecture and dynamics. Nat Rev Mol Cell Biol. 2008;9:478–89.
Mostowy S, Cossart P. Septins: the fourth component of the cytoskeleton. Nat Rev Mol Cell Biol. 2012;13:183–94.
Tada T, Simonetta A, Batterton M, Kinoshita M, Edbauer D, Sheng M. Role of Septin cytoskeleton in spine morphogenesis and dendrite development in neurons. Curr Biol. 2007;17:1752–8.
Xie Y, Vessey JP, Konecna A, Dahm R, Macchi P, Kiebler MA. The GTP-binding protein Septin 7 is critical for dendrite branching and dendritic-spine morphology. Curr Biol. 2007;17:1746–51.
Hu J, Bai X, Bowen JR, Dolat L, Korobova F, Yu W, et al. Septin-driven coordination of actin and microtubule remodeling regulates the collateral branching of axons. Curr Biol. 2012;22:1109–15.
Kang R, Wan J, Arstikaitis P, Takahashi H, Huang K, Bailey AO, et al. Neural palmitoyl-proteomics reveals dynamic synaptic palmitoylation. Nature. 2008;456:904–9.
Wan J, Savas JN, Roth AF, Sanders SS, Singaraja RR, Hayden MR, et al. Tracking brain palmitoylation change: predominance of glial change in a mouse model of Huntingto’s disease. Chem Biol. 2013;20:1421–34.
Xie Y, Zheng Y, Li H, Luo X, He Z, Cao S, et al. GPS-lipid: a robust tool for the prediction of multiple lipid modification sites. Sci Rep. 2016;6:28249.
Forrester MT, Hess DT, Thompson JW, Hultman R, Moseley MA, Stamler JS, et al. Site-specific analysis of protein S-acylation by resin-assisted capture. J Lipid Res. 2011;52:393–8.
Antinone SE, Ghadge GD, Ostrow LW, Roos RP, Green WN. S-acylation of SOD1, CCS, and a stable SOD1-CCS heterodimer in human spinal cords from ALS and non-ALS subjects. Sci Rep. 2017;7:41141.
Gauthier-Campbell C, Bredt DS, Murphy TH, El-Husseini AE-D. Regulation of dendritic branching and filopodia formation in hippocampal neurons by specific acylated protein motifs. Mol Biol Cell. 2004;15:2205–17.
Abdel-Maguid TE, Bowsher D. Alpha- and gamma-motoneurons in the adult human spinal cord and somatic cranial nerve nuclei: the significance of dendroachitectonics studied by the Golgi method. J Comp Neurol. 1979;186:259–69.
Paxinos G. The human nervous system. Cambridge: Academic Press; 2012.
Martin E, Cazenave W, Cattaert D, Branchereau P. Embryonic alteration of motoneuronal morphology induces hyperexcitability in the mouse model of amyotrophic lateral sclerosis. Neurobiol Dis. 2013;54:116–26.
Kato T, Hirano A, Donnenfeld H. A Golgi study of the large anterior horn cells of the lumbar cords in normal spinal cords and in amyotrophic lateral sclerosis. Acta Neuropathol. 1987;75:34–40.
Fu X-D, Ares M. Context-dependent control of alternative splicing by RNA-binding proteins. Nat Publ Group. 2014;15:689–701.
Polydorides AD, Okano HJ, Yang YY, Stefani G, Darnell RB. A brain-enriched polypyrimidine tract-binding protein antagonizes the ability of Nova to regulate neuron-specific alternative splicing. Proc Natl Acad Sci U S A. 2000;97:6350–5.
Buckanovich RJ, Yang YY, Darnell RB. The onconeural antigen Nova-1 is a neuron-specific RNA-binding protein, the activity of which is inhibited by paraneoplastic antibodies. J Neurosci. 1996;16:1114–22.
Schaefer A, Im H-I, Venø MT, Fowler CD, Min A, Intrator A, et al. Argonaute 2 in dopamine 2 receptor-expressing neurons regulates cocaine addiction. J Exp Med. 2010;207:1843–51.
Tan CL, Plotkin JL, Venø MT, von Schimmelmann M, Feinberg P, Mann S, et al. MicroRNA-128 governs neuronal excitability and motor behavior in mice. Science. 2013;342:1254–8.
Bahn JH, Ahn J, Lin X, Zhang Q, Lee J-H, Civelek M, et al. Genomic analysis of ADAR1 binding and its involvement in multiple RNA processing pathways. Nat Commun. 2015;6:6355.
Beltran M, Yates CM, Skalska L, Dawson M, Reis FP, Viiri K, et al. The interaction of PRC2 with RNA or chromatin is mutually antagonistic. Genome Res. 2016;26:896–907.
Hwang H-W, Park CY, Goodarzi H, Fak JJ, Mele A, Moore MJ, et al. PAPERCLIP identifies MicroRNA targets and a role of CstF64/64tau in promoting non-canonical poly(A) site usage. Cell Rep. 2016;15:423–35.
Hwang H-W, Saito Y, Park CY, Blachère NE, Tajima Y, Fak JJ, et al. cTag-PAPERCLIP reveals alternative polyadenylation promotes cell-type specific protein diversity and shifts Araf isoforms with microglia activation. Neuron. 2017;95:1334–5.
Pramatarova A, Laganière J, Roussel J, Brisebois K, Rouleau GA. Neuron-specific expression of mutant superoxide dismutase 1 in transgenic mice does not lead to motor impairment. J Neurosci. 2001;21:3369–74.
Lino MM, Schneider C, Caroni P. Accumulation of SOD1 mutants in postnatal motoneurons does not cause motoneuron pathology or motoneuron disease. J Neurosci. 2002;22:4825–32.
Boillée S, Vande Velde C, Cleveland DW. ALS: a disease of motor neurons and their nonneuronal neighbors. Neuron. 2006;52:39–59.
Phatnani HP, Guarnieri P, Friedman BA, Carrasco MA, Muratet M, O'Keeffe S, et al. Intricate interplay between astrocytes and motor neurons in ALS. Proc Natl Acad Sci U S A. 2013;110:E756–65.
Weyn-Vanhentenryck SM, Mele A, Yan Q, Sun S, Farny N, Zhang Z, et al. HITS-CLIP and integrative modeling define the Rbfox splicing-regulatory network linked to brain development and autism. Cell Rep. 2014;6:1139–52.
Gong S, Yang XW, Li C, Heintz N. Highly efficient modification of bacterial artificial chromosomes (BACs) using novel shuttle vectors containing the R6Kgamma origin of replication. Genome Res. 2002;12:1992–8.
Zhang C, Darnell RB. Mapping in vivo protein-RNA interactions at single-nucleotide resolution from HITS-CLIP data. Nat Biotechnol. 2011;29:607–14.
Darnell JC, Van Driesche SJ, Zhang C, Hung KYS, Mele A, Fraser CE, et al. FMRP stalls ribosomal translocation on mRNAs linked to synaptic function and autism. Cell. 2011;146:247–61.
Wu J, Anczuków O, Krainer AR, Zhang MQ, Zhang C. OLego: fast and sensitive mapping of spliced mRNA-Seq reads using small seeds. Nucleic Acids Res. 2013;41:5149–63.
Eden E, Navon R, Steinfeld I, Lipson D, Yakhini Z. GOrilla: a tool for discovery and visualization of enriched GO terms in ranked gene lists. BMC Bioinformatics. 2009;10:48.
Yuan Y, Phatnani H, Maniatis T, Darnell RB. Cell type-specific HITS-CLIP reveals differential RNA processing in motor neurons. NCBI GEO. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE71294. Accessed 4 July 2018.
The authors would like to thank Rada Norinsky, Roxana Cubias for oocyte microinjection, Nathaniel Heintz for pLD53.SC2 plasmid, Michael Moore, Christopher Park, Mariko Kobayashi for helpful discussion, and Scott Dewell for high-throughput sequencing. We wish to thank those who reviewed the manuscript for their constructive comments (Additional file 8).
This work was supported by the National Institutes of Health grant 5RC2NS069473–02. R.B.D. is a Howard Hughes Medical Institute Investigator.
Availability of data and materials
Raw .fastq files of RNA-seq and HITS-CLIP data have been deposited to the NCBI Gene Expression Omnibus (GEO) under accession number GSE71294 .
The review history is available as Additional file 8.
Ethics approval and consent to participate
All animal experiments were performed in the Association for Assessment and Accreditation of Laboratory Animal Care (AAALAC)—accredited Animal Resource Center at the Rockefeller University under protocol numbers 07111, 14,678, and 17,013.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplemental material and methods. (DOCX 115 kb)
Figure S1. GFP-Nova2 retains Nova2 RNA binding preference. Figure S2. HITS-CLIP on BAC transgenic spinal cords and WSC. Figure S3 Data analysis method. Figure S4. Bioinformatic analyses of disproportionally different NOVA binding sites in MNs. Figure S5. Additional cases where MN-specific NOVA binding correlates with MN-specific alternative splicing. Figure S6. NOVA regulates Sept8 ALE usage. Figure S7. Illustration of palmitoylation detection assay; SEPT8 expression in hippocampal cultures; SEPT8-X1 and NOVA2 knockdown in cell lines. Figure S8. NOVA binding on genes encoding cytoskeleton interacting proteins are over-represented in MNs. (PDF 30346 kb)
a Sequencing and mapping statistics of MN and WSC NOVA CLIP. b Sequencing and mapping statistics of MN and WSC NOVA CLIP. c mm9 coordinates and genomic annotation of WSC NOVA CLIP peaks. d List of known NOVA targets harboring WSC NOVA CLIP peaks. e List of known NOVA targets harboring MN NOVA CLIP peaks. f GO enrichment of known MN versus WSC NOVA targets. g Transcripts with enriched NOVA CLIP signals in MNs compared to WSC. h Transcripts with depleted NOVA CLIP signals in MNs compared to WSC. (XLSX 13044 kb)
a Differential analysis results for all intronic joint peaks. b Differential analysis results for all intronic joint peaks. (XLSX 14058 kb)
a Transcript regions differentially spliced between MNs and WSC. b For the nine exons with over- or under-represented MN NOVA binding sites in regulatory “hotspots”, correlation between observed and predicted dI between MNs and WSC. (XLSX 365 kb)
a List of alternative last exons differentially utilized in MNs vs WSC. b List of alternative last exons differentially utilized in wild-type vs Nova2 knockout mouse brain. (XLSX 69 kb)
List of RNA linkers and DNA primers used in CLIP. (XLSX 70 kb)
Review history for this manuscript. (DOCX 74 kb)
About this article
Cite this article
Yuan, Y., Xie, S., Darnell, J.C. et al. Cell type-specific CLIP reveals that NOVA regulates cytoskeleton interactions in motoneurons. Genome Biol 19, 117 (2018). https://doi.org/10.1186/s13059-018-1493-2