Regulatory consequences of neuronal ELAV-like protein binding to coding and non-coding RNAs in human brain

1RNA binding proteins (RBPs) associate with RNAs throughout their life cycle, regulating all aspects of RNA metabolism and function. More than 800 RBPs have been described in human cells. The unique structure and function of neurons, and the need to rapidly adapt RNA regulation in the brain both within and at sites distant from the nucleus, are consistent with specialized roles for RBPs in the brain. Indeed, mammalian neurons have developed their own system of RNA regulation, and RBP:mRNA interactions are thought to regulate local protein translation at synapses, perhaps underlying learning and long-term memory.

2Numerous RBPs have been linked to human neurological disorders (reviewed in ). For example, FUS, TDP-43 and ATXN2 mutations have been found in familial amyotrophic lateral sclerosis patients, TDP-43 has additionally been associated with frontotemporal lobar degeneration, Alzheimer’s Disease (AD) and Parkinson’s Disease (PD), STEX has been linked to amyotrophic lateral sclerosis 4, and spinal muscular atrophy can be caused by mutations in SMN.

3The neuronal ELAV-like (ELAVL) and NOVA RBPs are targeted by the immune system in paraneoplastic neurodegenerative disorders. Mammalian ELAVL proteins include the ubiquitously expressed paralog ELAVL1 (also termed HUA or HUR) and the three neuron-specific paralogs, ELAVL2, 3 and 4 (also termed HUB, C, and D, and collectively referred to as nELAVL; ). nELAVL proteins are expressed exclusively in neurons in mice, and they are important for neuronal differentiation and neurite outgrowth in cultured neurons. Redundancy between the three nELAVL isoforms complicates in vivo studies of their individual functions. Nevertheless, even haploinsuffiency of Elavl3 is sufficient to trigger cortical hypersynchronization, and Elavl3 and Elavl4 null mice display defects in motor function and neuronal maturation, respectively.

4ELAVL proteins have been shown to regulate several aspects of RNA metabolism. In vitro and in tissue culture cells, nELAVL proteins have been implicated in the regulation of stabilization and/or translation of specific mRNAs, as well as in the regulation of splicing and polyadenylation of select transcripts [reviewed in ]. A more comprehensive approach was taken by immunoprecipitating an overexpressed isoform of ELAVL4 in mice, although such RNA immunoprecipitation experiments cannot distinguish between direct and indirect targets. Recently, direct binding of nELAVL to target RNAs in mouse brain was demonstrated by high-throughput sequencing of RNA isolated by crosslinking immunoprecipitation (HITS-CLIP; ); these data, coupled with transcriptome profiling of Elavl3/4 KO mice, demonstrated that nELAVL directly regulates neuronal mRNA abundance and alternative splicing by binding to U-rich elements with interspersed purine residues in 3’UTRs and introns in mouse brain.

5While genome-wide approaches have been applied to studying nELAVL proteins in mice, the targets of nELAVL in the human brain remain largely unknown. This is of particular importance, as nELAVL proteins have been implicated in neurological disorders such as AD and PD. Hence, to advance our understanding of the function of nELAVL in humans and its link to human disease, we set out to investigate nELAVL:RNA interactions in the human brain.

6To globally identify transcripts directly bound by nELAVL in human neurons, we generated a genome wide RNA binding map of nELAVL in human brain using CLIP. CLIP allows the identification of functional RNA-protein interactions in vivo by using UV-irradiation of intact tissues to covalently crosslink and then purify RNA-protein complexes present in vivo. This method has been adopted for a variety of RBPs. Here, we systemically identified tens of thousands of reproducible nELAVL binding sites in human brain and showed that nELAVL binds transcripts that are important for neurological function and that have been linked to neurological diseases such as AD. We validated the functional consequences of nELAVL binding in mice and cultured human neuroblastoma cells and showed that the loss of nELAVL affected mRNA abundance and alternative splicing of hundreds of transcripts. We further investigated RNA regulation in AD brains, and found that numerous transcripts were differentially spliced in AD, which correlated with differential nELAVL binding in some cases. Remarkably, we observed the most significant increase in nELAVL binding in AD on a class of non-coding RNAs, Y RNAs. We recapitulated these findings in human neuroblastoma cells, showing that nELAVL binding is linked to Y ribonucleoprotein (RNP) remodeling acutely during UV-induced stress, and chronically in AD.

7To gain insight into nELAVL-mediated RNA regulation in human brain we performed CLIP on postmortem brain samples of eight human subjects (Supplementary file 1A). Tissue samples were derived from BA9, which is part of the dorsolateral prefrontal cortex (Figure 1A), a brain area that is damaged in later stages of AD and that is important for executive functions such as working memory, cognitive flexibility, planning, inhibition, and abstract reasoning. Antibodies that specifically recognize individual ELAVL paralogs and that can be used for CLIP are currently not available. We therefore purified nELAVL-RNP complexes with antiserum reactive to all three nELAVL proteins (Figure 1B). 32P-labeled nELAVL-RNP complexes were not recovered with control serum or in the absence of UV-irradiation (Figure 1B).

8nELAVL-crosslinked RNA tags were sequenced and mapped to the hg18 build of the human genome. We searched for and identified nELAVL RNA binding sites (peaks) that were significant (p<0.01) and that were present in at least five out of 8 individuals (n = 75,592). nELAVL binding at these sites correlated between individuals (Figure 1—figure supplement 1), and these peaks were further investigated. We determined the nELAVL binding motif by analyzing the top 500 nELAVL peaks (+/- 25nt) using MEME ChIP and HOMER. This revealed that nELAVL binds polyU RNA stretches in human brain, particularly when interrupted by a G (Figure 1C and Figure 1—figure supplement 2). This is in excellent agreement with the nELAVL binding motif identified in mouse brain using CLIP and in vitro binding assays. Because nELAVL has been shown to regulate alternative splicing and mRNA abundance in mouse brain by binding to introns and 3’UTRs, respectively, we analyzed the genomic distribution of the peaks defined here. As previously reported for mouse nElavl, nELAVL binding sites are found in 3’UTRs and introns (Figure 1D), with far higher per-nucleotide density in 3’UTR regions (Figure 1—figure supplement 3). Consistently we observed that nELAVL peaks in 3’UTRs were higher than intronic peaks (Figure 1—figure supplement 4). These results suggest that nELAVL could regulate splicing and mRNA abundance in the human brain.

9To relate nELAVL binding to mRNA abundance, we performed RNAseq on the same brain samples used for CLIP analysis (Figure 1—figure supplement 5 and Supplementary file 1A/B). 74,423 nELAVL peaks mapped to 8681 expressed genes (Supplementary file 1C), referred to as nELAVL targets (Supplementary file 1D) hereafter, which are shown in Figure 1E. We observed that nELAVL binding correlated with mRNA abundance (Figure 1E). This was not unexpected, as nELAVL 3’ UTR binding has previously been shown to increase mRNA abundance, due to its role in mRNA stabilization. The correlation between nELAVL binding and mRNA abundance might therefore not only reflect the dependence of nELAVL binding on mRNA abundance, but also a role of nELAVL in mRNA stabilization. Consistently, we observed that intronic nELAVL binding correlated less with mRNA abundance than 3’UTR binding (Figure 1—figure supplement 6).

10To identify genes most likely to be impacted by nELAVL, we defined the top 1000 nELAVL targets (colored in green in Figure 1E; Supplementary file 1E). Top targets were identified based on normalized nELAVL binding (binding sites were normalized for mRNA abundance and summarized per gene). Thirty-seven percent of nELAVL peaks (n = 27,581) mapped to these top targets. We constructed a subnetwork that connects top nELAVL target gene products based on a literature-based network of protein-protein interactions (PPI) created from multiple online databases. Six clusters were identified within the resulting network (Figure 1F) and gene set enrichment analyses were performed for top nELAVL targets found in the different clusters with Enrichr. Each cluster was examined for enrichment of Biological Processes (BP), Molecular Functions (MF), Cellular Components, and Mammalian Phenotypes (MP) terms (Supplementary file 1F). Enriched terms included RNA processing and transcription regulation, signal transduction, synaptic transmission, synaptic proteins, and abnormal neuron morphology and physiology. Three of the six clusters were particularly important for neuronal function. Cluster 1 is especially enriched in members of the TGFbeta/SMAD signaling pathways and the FOX protein family. Cluster 2 contains many actors of the IGF-I axis, which is important for neuronal development including neurogenesis, myelination, synaptogenesis, dendritic branching and neuroprotection after neuronal damage. Finally, cluster 3 is almost exclusively formed of synaptic proteins including many postsynaptic scaffolding proteins, members of the neuroligin/neurexin families and glutamatergic receptors or voltage-gated channels. Taken together, these data demonstrate that nELAVL associates with transcripts encoding proteins involved in key aspects of neuronal physiology.

11To investigate if nELAVL-mediated RNA regulation is conserved, we compared our dataset with previously published nELAVL targets in mice. At the transcript level, we found that more than 90% of mouse nELAVL targets were among human nELAVL targets (Figure 2A). However, only 20% of human targets were bound by nELAVL in mouse brain, which is at least partly due to the 10-fold increased depth of the human dataset (more than 10 million human CLIP tags compare to less than a million mouse CLIP tags). Yet these differences could also reflect an increased functional complexity of nELAVL regulation in the human brain, and/or the fact that mouse targets were identified at different developmental stages.

12We additionally investigated the overlap of nELAVL binding at individual binding sites. Surprisingly, only a small percentage of binding sites overlapped between mouse and human; 3% of human binding sites showed nELAVL binding in mouse, and 17% of mouse binding sites were bound by nELAVL in human brain (Figure 2B). The vast majority of these overlapping binding sites were in 3’UTRs (88%). These results indicate that many nELAVL targets are shared between mouse and human and that nELAVL binding at the transcript level is conserved, whereas individual binding sites have diverged drastically, especially within introns. These results reflect analogous observations of evolutionary conservation of transcriptional regulation at the gene rather than the positional level.

13We further observed that nELAVL binding on entire transcripts correlated between mouse and human (Figure 2—figure supplement 1), which prompted us to overlay our human CLIP dataset with previously published transcriptome profiling of Elavl3/4 double KO mice. 119 transcripts showed significant changes in their steady-state level in Elavl3/4 KO mice, 91 of which were expressed in human brain. 37 of these 91 transcripts were nELAVL 3’UTR targets in human brain (Supplementary file 2A), and the majority of them decreased in the absence of ELAVL3/4 (n = 26), including transcripts important for neuronal transport and excitation such as RAB6B, HCN3, and KCNMB2 (Figure 2C/D). This indicates that nELAVL 3’UTR binding is likely to be important for increasing the abundance of these transcripts in human brain, and likely has conserved functions across species.

14Elavl3/4 KO mice also reported splicing defects and 59 alternative exons showed a significant change in their inclusion rate (delta inclusion rate, ΔI) between wildtype and Elavl3/4 KO mice. 54 of the misregulated exons were conserved in the human genome, and 25 of them were adjacent to intronic nELAVL binding sites in human brain (Supplementary file 2B). We observed both increased and decreased inclusion of alternative exons – independently of the position of the peak relative to the exon, which has previously been observed for nELAVL mediated splicing regulation. Three alternative exons are shown in Figure 2E,F, and whereas nELAVL seems to prevent splicing of DST by binding upstream and downstream of an alternative exon, nELAVL might promote the inclusion of alternative exons of NRXN1 and CELF2 by binding to downstream sequences. Given that nELAVL regulates the splicing of these 25 exons in mice and that we observe intronic nELAVL binding sites in human brain adjacent to these exons, we propose that nELAVL regulates the inclusion of these exons in human brain. Collectively, these analyses show that many confirmed functional nELAVL interactions in mouse brain show evidence for nELAVL binding in human brain.

15To further validate potential nELAVL targets, we analyzed the effect of nELAVL depletion on mRNA abundance and splicing in human neuroblastoma IMR-32 cells. We subjected IMR-32 cells to mock or ELAVL2/3/4 triple RNAi, achieving 70% knockdown of all three neuronal ELAVL proteins (Figure 3—figure supplement 1). These cells were then analyzed by RNAseq (Figure 3A, Supplementary file 1A/2C). The steady-state level of 784 transcripts was significantly changed in nELAVL RNAi treated cells (Figure 3A), with ~45% showing a decrease in mRNA abundance. Among those genes were all of the neuronal ELAVL paralogs (ELAVL2/3/4), while the ubiquitously expressed paralog ELAVL1 was not affected by nELAVL RNAi depletion.

16We then compared this RNAseq dataset with the nELAVL CLIP analysis in human brain. 96% of the IMR-32 expressed transcripts were expressed in human brain (12,242 out of 12,743), and were used for subsequent analyses. Since nELAVL binding to 3’UTRs can mediate mRNA stabilization, we investigated the change in mRNA abundance of 3’UTR targets upon nELAVL RNAi, specifically examining 3’UTR targets that did not display any intronic binding. These transcripts were less abundant in nELAVL RNAi conditions (Figure 3B, left panel). In contrast, the mRNA abundance of intron targets (intron binding but no 3’UTR binding) slightly increased upon nELAVL RNAi (Figure 3B, right panel). This suggests that specifically nELAVL binding to 3’UTRs increases mRNA abundance. We further observed that nELAVL depletion affected top nELAVL 3’UTR targets as well as nELAVL 3’UTR targets in general (Figure 3C).

17Out of 784 genes that changed significantly upon nELAVL RNAi, 743 genes were expressed in human brain, 327 of which decreased while 416 increased. We investigated which of these transcripts were direct targets of nELAVL based on nELAVL 3’UTR binding (Supplementary file 2D). Significantly changing transcripts that are top nELAVL 3’UTR targets are boxed in blue in Figure 3D. We observed that 68% of downregulated transcripts were 3’UTR targets (n = 226; p = 9.6e-13; hypergeometric test), and that 16% of downregulated transcripts were even among top 3’UTR targets (n = 51; p = 1.3e-6; hypergeometric test). In contrast, only 7% of upregulated transcripts were among top 3’UTR targets (p = 0.76; hypergeometric test), further supporting a role of nELAVL 3’UTR binding in positively regulating mRNA abundance. Several 3’UTR targets showed an mRNA abundance change in both mouse and IMR-32 datasets (n = 8), and in all but one case this change correlated positively between the datasets, providing support for the accuracy of target validation applied here. Because the abundance of multiple disease associated genes, including APPBP2, ATXN3, and SHANK2 (Figure 3E,F), is regulated by nELAVL, we propose that nELAVL mediated regulation of mRNA abundance plays an important role in the human brain.

18To validate a role of nELAVL in the splicing regulation of human brain targets, we analyzed the inclusion rate of cassette exons in mock and nELAVL RNAi treated IMR-32 cells. We compared the change in exon inclusion of 7903 expressed cassette exons (Supplementary file 2E) and observed that 473 cassette exons were differentially spliced upon nELAVL depletion (FDR<0.05 and ΔI > 0.1; Supplementary file 2F). Many differentially spliced exons were adjacent (+/- 2.5 kb) to at least one intronic nELAVL binding site in human brain (n = 155; p = 1.3e-7; hypergeometric test; Figure 4A, Supplementary file 2F), indicating that these exons might be directly regulated by nELAVL. For example, downstream binding in BIN1 and PICALM was associated with lower exon inclusion upon nELAVL depletion, and binding in APP was associated with higher inclusion of both upstream and downstream exons upon nELAVL depletion (Figure 4B/C). Overall, three exons that were differentially spliced upon nELAVL RNAi depletion also changed in Elavl3/4 KO mice, and the splicing changes in both datasets changed in the same direction. We generated a map from intronic nELAVL binding sites that flanked the 155 nELAVL regulated exons as previously described, revealing that upstream nELAVL binding can promote both exon inclusion and skipping (Figure 4D). In conclusion, these data indicate that intronic nELAVL binding regulates alternative splicing of numerous transcripts in human brain, including transcripts associated with central nervous system disorders.

19nELAVL has previously been linked to neurological diseases and we observed that nELAVL regulated the mRNA abundance and splicing of multiple disease-associated genes. We examined nELAVL binding in a set of genes with disease associated 3’UTR single nucleotide polymorphisms (SNPs). We found that these genes were enriched among nELAVL 3’UTR targets (n = 200; p = 0.001; hypergeometric test), and that nELAVL binding sites directly overlapped with 45 disease associated SNPs, including SNPs associated with autism, schizophrenia, depression, AD, and PD (Figure 5—figure supplement 1, Supplementary file 3A).

20nELAVL proteins have been implicated in AD, and among the validated nELAVL regulated RNAs were also several AD-related transcripts, which led us to investigate additional AD-linked genes (hereafter termed AD genes; n = 96; Supplementary file 3B). Indeed, we found that the top nELAVL targets were enriched among AD genes (n = 11; p = 0.03; hypergeometric test; contained in Supplementary file 3B) as well as among AD risk loci identified in a genome-wide association study (GWAS) in AD (n = 77; p = 1.7e-14; hypergeometric test; Supplementary file 3C). To investigate if nELAVL mediated regulation of AD related and other transcripts might be affected in AD, we performed nELAVL CLIP and RNAseq on AD subject brains, age-matched to control subjects (Figure 5—figure supplement 2, Supplementary file 1A/B and 3D). Importantly, ELAVL3/4 mRNA levels were similar between control and AD samples and ELAVL2 showed only a slight decrease in transcript abundance in AD brains (Supplementary file 1B), which allowed us to compare nELAVL binding profiles between control and AD brains. We did not detect many significant changes in nELAVL binding nor mRNA abundance (Figure 5A/B, Supplementary file 1B and 3D), probably due to the variation between human samples, the small sample size, and the potential heterogeneity of AD. We did however observe that 150 transcripts were differentially spliced in the 9 AD subjects (FDR<0.05 and ΔI>0.1; Figure 5C, Supplementary file 3E). Two of these transcripts, BIN1 and PTPRD, have previously been linked to AD, suggesting that the differential splicing of these two transcripts as well as other RNAs might be linked to AD.

21As shown above (Figure 4), nELAVL depletion in IMR-32 cells was associated with the reduced inclusion of an alternative exon of BIN1, suggesting that nELAVL binding promotes the inclusion of this exon. Precisely this exon was differentially spliced in AD subjects, with AD subjects showing a reduced exon inclusion rate compared to control subjects (Figure 5D). Along with the differential exon inclusion, we observed that nELAVL peak binding was fourfold decreased in AD subjects (log2 fold change = -2.35; p = 0.16; Figure 5D). These results are consistent with nELAVL-mediated dysregulation of this exon in AD subjects, with decreased binding leading to decreased exon inclusion. In conclusion, while we did not detect global nELAVL binding and mRNA abundance changes in AD subjects, we observed that splicing of 150 transcripts was affected, which in some cases might be linked to nELAVL dysregulation.

22The largest fold changes in nELAVL binding in AD (relative to the age-matched control population) occurred on a specific class of non-coding RNAs, Y RNAs. Y RNAs are 100 nt long structured RNAs usually found in complex with RO60 (also known as TROVE2; Figure 6A; modified from ). RO60 is believed to act as a sensor of RNA quality, targeting defective RNAs for degradation. RO60 was initially identified as an autoantigen targeted in systemic lupus and some subjects with the paraneoplastic encephalopathy syndrome harbor both anti-RO and anti-nELAVL (Hu) autoantibodies. Four canonical Y RNAs, Y1/3/4/5, have been characterized in humans, but numerous slightly divergent copies of these Y RNAs, especially Y1 and Y3, are distributed throughout the human genome.

23Surprisingly, we observed nELAVL binding to a total of 320 Y RNAs, although Y RNA copies other than the canonical four Y1/3/4/5 genes had previously been considered to be non-functional and were labeled 'pseudogenes' (Supplementary file 3F). We found that 237 of the 320 nELAVL bound Y RNAs were Y3-like RNAs (Supplementary file 3F), and that nELAVL bound Y RNAs showed an enrichment of the nELAVL binding motif (202 Y RNAs contained UUUUUU, allowing a G at any one position), which is also present in the canonical hY3 RNA (Figure 6A/B). We examined the 118 nELAVL bound Y RNAs that did not fit this consensus in more detail. 91 of these Y RNAs (77%) contained either a 5mer version of the motif or the motif with an A or C instead of a G, and we found U/G rich stretches in the remaining 27 Y RNAs (Supplementary file 3F). In addition, some Y RNAs with a strong binding motif did not show any evidence of nELAVL binding. In general, these Y RNAs showed a lower expression compared to nELAVL bound Y RNA, which may explain the absence of detectable nELAVL binding (Figure 6—figure supplement 1).

24We next explored nELAVL/Y RNA binding in AD brain. We observed a drastic increase in nELAVL/Y RNA association in AD subjects (Figure 6C), while Y RNA levels remained largely unchanged (Figure 6D). This suggests that Y RNPs undergo nELAVL-dependent remodeling in AD. Interestingly, we did observe a high variability in nELAVL/Y RNA association between AD samples (Figure 6—figure supplement 2), with three of them showing a very strong nELAVL/Y RNA association. Efforts to relate this difference to the expression of stress-related genes, post-mortem interval, age, extent of disease and cause of death were not conclusive, and the cause for the variation in nELAVL binding to Y RNAs among AD subjects remains elusive.

25The observation of increased nELAVL/Y RNA association in AD raised the possibility that Y RNP remodeling is associated with neuronal stress. Y RNP remodeling has previously been linked to UV-induced stress, and both bacterial and mouse cells show an increased sensitivity to UV stress in the absence of RO60. ELAVL binding can be modulated in response to stress in cultured cells, and ELAVL proteins, which shuttle between nucleus and cytoplasm in response to environmental cues, preferentially accumulate in cytoplasmic stress granules upon stress. We therefore examined the effect of acute UV stress on Y RNP remodeling in IMR-32 cells. IMR-32 cells were exposed to a low dose of UV stress (not sufficient to induce RNA:protein crosslinking) and allowed to recover for 24 h before being analyzed by nELAVL CLIP. We found that nELAVL bound 132 Y RNAs in neuroblastoma cells (Supplementary file 3F), that Y RNAs showed an enrichment of the nELAVL binding motif (Figure 7A) or at least contained a degenerate version of it (Supplementary file 3F), and that non-bound Y RNAs with a motif show a very low expression (Figure 7—figure supplement 1). Moreover, nELAVL binding on Y RNAs was dynamic and increased in UV stressed cells compared to non-stressed cells (Figure 7B and Figure 7—figure supplement 2), while their abundance did not change upon UV irradiation (Figure 7C). To assess whether Y RNA levels were affected by nELAVL, we depleted nELAVL by RNAi three days prior to the UV exposure, and analyzed Y RNA levels by RNAseq. Y RNA abundance was not affected by nELAVL depletion in UV stressed IMR-32 cells (Figures 7D). These results indicate that increased nELAVL binding to Y RNAs is not a function of Y RNA levels, and that nELAVL binding during stress is not required for Y RNA stability.

26To determine if UV stress induced localization changes of Y RNP components, we investigated the distribution of nELAVL, RO60 as well as Y RNAs upon UV exposure using cell fractionation followed by western blot and qPCR analysis. The induction of a UV stress response was confirmed by measuring CDKN1A mRNA levels (Figure 7—figure supplement 3A). We did not observe a change in nucleocytoplasmic localization of the investigated RNAs or proteins (Figure 7—figure supplement 3B/C), suggesting that the increased nELAVL/Y RNA association upon UV exposure does not result from a difference in the nucleocytoplasmic distribution of nELAVL or Y RNAs. These results are consistent with previous observations that neuronal ELAVL proteins show a higher cytoplasmic localization than the ubiquitous paralog ELAVL1, and that stress-induced nuclear-cytoplasmic shuttling might be limited to ELAVL1. Nonetheless, these results do not rule out the possibility that there may be changes of nELAVL proteins within the nuclear or cytoplasmic compartments themselves with respect to Y RNA binding and localization.

27We next sought to measure the proportion of nELAVL bound to Y RNAs in stressed and non-stressed conditions. Because Y RNAs are relatively short and have a high degree of similarity, our mapping strategy (reporting only unambiguous mapping events) discarded numerous reads that were assigned to multiple Y RNAs. We therefore re-mapped CLIP tags, allowing multiple alignments, but reporting only the best match, permitting a more accurate estimate of overall Y RNA binding. The fraction of the short CLIP reads mapping to Y RNAs was considerably higher using this strategy, revealing that up to 6% of nELAVL CLIP tags map to Y RNAs in AD and UV stressed cells, compared to less than 0.5% in control brain and ~1% in non-stressed cells (Figure 7—figure supplement 4). The significant increase in nELAVL/Y RNA association and our observation that up to 6% of nELAVL was bound to Y RNAs might in fact lead to a sequestration of nELAVL from its targets.

28To investigate if the increased nELAVL/Y RNA association was linked to decreased intronic and 3’UTR nELAVL binding, we grouped AD subjects based on their Y RNA association and compared the two different AD groups to control subjects. We found that the majority of changing nELAVL binding sites decreased in AD subjects with high nELAVL/Y RNA, while nELAVL binding in AD subjects with low nELAVL/Y RNA association was mostly increasing (Figure 8A). Because nELAVL binding in UV-stressed cells also predominantly decreased (Figure 8A) and most of the decreased binding sites were in introns (85%, assessed by annotation of peak locations), we speculate that nELAVL/Y RNA association leads to a sequestration of nELAVL specifically from its intron targets, which might induce similar splicing changes as nELAVL depletion by RNAi. Of note is our observation that nELAVL binding decreased at only a subset of intronic binding sites.

29We further examined the possibility of Y RNA mediated nELAVL sequestration upon UV stress by subjecting mock and nELAVL RNAi treated IMR-32 cells to UV stress and analyzing these cells by RNAseq. nELAVL mRNA and proteins levels decreased upon nELAVL RNAi in non-stressed and UV stressed cells but were not affected by UV stress (Figure 8—figure supplement 1 and Supplementary file 2C). We analyzed exon inclusion rates as above (Figure 4) and found that 9397 cassette exons were expressed between the four conditions (Supplementary file 2E). Comparing UV-induced splicing changes between mock and nELAVL RNAi treated cells, we identified 260 cassette exons that showed a differential inclusion rate upon UV stress only in the presence of nELAVL (Figure 8B). We intersected these splicing changes with nELAVL RNAi induced splicing changes (n = 553), and found a significant overlap between the two lists (n = 68; p = 9.7e-28; hypergeometric test; Figure 8C and Supplementary file 2E). Importantly, splicing of the vast majority of exons (66 out of 68) changed in the same direction, indicating that changes in nELAVL binding due to UV stress partially recapitulate nELAVL RNAi depletion. This finding is consistent with a model of UV-induced nELAVL sequestration from a subset of targets.

30Due to the wide difference between AD subjects and UV stressed IMR-32 cells, the targets affected by nELAVL sequestration in the two systems are likely to be markedly divergent. We nevertheless explored if any of the 66 affected exons in IMR-32 cells were adjacent (+/- 2.5 kb) to intronic nELAVL binding sites in human brain, and observed that 7 alternatively spliced exons were indeed next to intronic nELAVL binding sites (Figure 8D/ Figure 8—figure supplement 2). Remarkably, nELAVL binding at all of these sites decreased in AD subjects with Y RNA association but not in AD subjects without Y RNA association (Figure 8D), although nELAVL peak binding changed significantly at only one binding site (boxed in Figure 8E, FDR = 0.003). We also investigated if decreased nELAVL peak binding was associated with corresponding splicing changes in AD subjects (Figure 8—figure supplement 2). The inclusion rate of only one of the 7 exons changed significantly in AD patients with high nELAVL/Y RNA association (FDR = 0.04), and the direction of the splicing change was indeed the same as the splicing changes observed upon UV or nELAVL RNAi treatment. This is consistent with our observation that only few splicing changes are shared between AD and UV treatment, which likely reflects the much more complex situation in human brain.

31To directly test the model of Y RNA mediated nELAVL sequestration, we overexpressed canonical Y3wt (wild type) and Y3mut (with a mutated nELAVL binding site) using lentiviral infections of IMR-32 cells. We initially confirmed the overexpression of the infected Y RNAs using qPCR, which was more pronounced in the Y3mut infected cells (Figure 9A). The distribution of neither nELAVL nor RO60 was affected upon infection (Figure 9—figure supplement 1). To evaluate the extent of Y3wt and Y3mut overexpression compared to endogenous Y3 RNA expression we additionally analyzed infected cells by RNAseq (Figure 9B, Figure 9—figure supplement 2, Supplementary file 3F). While we consistently observed an increase in the Y3-like RNA expression upon infection, the magnitude of overexpression was modest relative to the endogenous expression of Y RNA copies. Nonetheless, we observed a small increase in total Y3-like RNAs in Y3wt but not Y3mut infected cells (Figure 9B).

32We next investigated mRNA abundance and splicing changes upon Y RNA overexpression (Supplementary file 2C/E). The mRNA abundance changes upon Y3wt and Y3mut infection compared to non-infected controls were very similar, indicating that most mRNA abundance changes are due to lentiviral infection (Figure 9—figure supplement 3A; 70% of Y3wt changes overlapped with Y3mut changes; p = 1.3e-175; hypergeometric test). To investigate virus-independent changes, we focused on the changes between Y3wt and Y3mut infection. nELAVL 3’UTR targets were not enriched among mRNAs that changed between Y3mut compared to Y3wt infected cells (Figure 9—figure supplement 3B), which is in agreement with our hypothesis that nELAVL sequestration predominantly affects intronic nELAVL binding sites and thus nELAVL mediated splicing.

33In contrast to the mRNA abundance changes, only few splicing changes overlapped between Y3wt and Y3mut infection when compared to non-infected cells (17% of Y3wt induced changes overlapped with Y3mut induced changes). Most of the observed splicing changes are therefore likely to be specific to Y RNA overexpression. Importantly, we observed an enrichment of nELAVL bound exons and of nELAVL RNAi dependent exons among the exons that changed upon Y3wt but not Y3mut overexpression (Figure 9C/D and Figure 9—figure supplement 4, 5). The relatively small enrichment is consistent with the modest increase in total Y3-like Y RNAs. These results suggest that Y RNA overexpression results in nELAVL sequestration from some of its intronic targets and consequent splicing changes, and partially recapitulates the stress induced nELAVL sequestration due to increased nELAVL/Y RNA association seen in AD patients and UV treated IMR-32 cells.

34nELAVL proteins are abundant neuron-specific RNA binding proteins which have been suggested to regulate various neurological processes and have been linked to neurodegenerative disorders including AD and PD. Yet the RNA targets of nELAVL in human brain were completely unknown. Here, we generated a comprehensive genome-wide RNA binding map of nELAVL in human brain, identifying 75,592 significant binding events within 8681 transcripts. We observed a significant overlap between these binding sites and disease-associated 3’UTR SNPs, and the potential disruption of nELAVL-mediated RNA regulation at these sites might contribute to disease manifestation. Most deleterious variants to date have been identified by exome sequencing while as many as 50% of disease-causing mutations are thought to affect splicing. With whole genome sequencing being increasingly available, non-coding variants are also increasingly detected, some of which may be linked to disease. As the majority of nELAVL binding occurs in introns and 3’UTRs, we expect that many binding sites will overlap with prospective disease-associated non-coding variants. The overlap between deleterious variants and nELAVL binding sites, and the observation that nELAVL binding at individual sites diverged between mice and human, underscores the importance of this study and illustrates the caveat of relying solely on mouse models when studying human disease. Considering the widespread nature of nELAVL binding in human brain and that RNA dysregulation has been linked to numerous neurological disorders, we believe that this binding map will be a valuable resource for the scientific community.

35To analyze the functional consequences of nELAVL binding, we used two different loss-of-function models: Elavl3/4 KO mice and nELAVL RNAi depletion in neuroblastoma cells. Due to the incomplete RNAi depletion of nELAVL in neuroblastoma cells, and potential differences in mRNA abundance and therefore nELAVL binding between the different samples, it is likely that we validated only a fraction of nELAVL-regulated transcripts. Despite these technical limitations we demonstrated that nELAVL impacts mRNA abundance and/or splicing of hundreds of targets. Among the nELAVL regulated transcripts were many transcripts implicated in human disease, including AD, which led us to investigate RNA regulation in AD subjects. Due to the relatively small sample size and the heterogeneity between these samples, likely due to both differences between individuals and sample preservation during postmortem collection, we did not detect many reproducible changes in mRNA abundance or nELAVL binding between AD and non-AD subjects. However, we found that 150 transcripts were differentially spliced in AD subjects, which in some cases coincided with differential nELAVL binding. Unexpectedly, the most significant binding change in AD was a dramatic increase in nELAVL binding to a class of non-coding RNAs, termed Y RNAs. This change was evident on a specific subset of Y RNAs harboring the nELAVL binding site. nELAVL/Y RNA binding also increased during UV stress in human neuroblastoma cells, while the abundance of Y RNAs remained constant in AD subjects and upon UV exposure. The increased nELAVL/Y RNA association correlated with decreased nELAVL binding at a subset of intronic binding sites, and was associated with similar splicing changes as induced by nELAVL depletion, suggesting that nELAVL/Y RNP remodeling during acute and chronic stress sequesters nELAVL from its mRNA targets. We provided further evidence for a Y RNA dependent nELAVL sequestration by overexpressing Y3 RNAs harboring either a wild type or mutated nELAVL binding site. Exons that were differentially spliced upon Y RNA overexpression were enriched for nELAVL bound exons, indicating nELAVL sequestration, which was dependent on an intact nELAVL binding site in the Y RNA.

36nELAVL 3’UTR binding has been implicated in increasing mRNA abundance in vivo. We described numerous nELAVL 3’UTR targets in brain, and were able to validate many of these targets, including disease-associated transcripts, indicating that nELAVL 3’UTR binding is important for the regulation of mRNA abundance in human brain. While ELAVL binding is frequently reported to result in an increase in mRNA abundance, we found several cases where nELAVL binding seemed to have an opposing effect. ELAVL proteins can compete or collaborate with miRNAs as well as RBPs like AUF1, CUGBP1 and TIA1 to regulate its targets. The ultimate outcome of nELAVL 3’UTR binding might therefore vary between individual transcripts.

37nELAVL has also been shown to regulate splicing in mouse brain by binding to intronic sequences. We observed many instances of intronic nELAVL binding events adjacent to alternative exons in brain, and confirmed that nELAVL regulates many of these exons in mice and neuroblastoma cells. In contrast to the position-dependent splicing observed for other RBPs, we observed that upstream nELAVL binding was associated with both exon skipping and inclusion. While nELAVL binding was observed within 25-50 nucleotides upstream of skipped exons, coinciding with the branch point sequence, nELAVL binding peaked within the proximal 25 nucleotides upstream of included exons, overlapping the polypyrimidine tract. Binding of auxiliary splicing factors, including nELAVL, to the branch point sequence usually interferes with spliceosome assembly and thus leads to exon skipping. Polypyrimidine tract binding however can lead to both exon inclusion and skipping, presumably depending on the recruitment of splicing enhancers or silencers. Our data indicates that upstream nELAVL binding can both interfere with the assembly of the spliceosome as well as promote splicing, most likely by recruiting splicing enhancers.

38Splicing defects have been associated with many neurological diseases, and among the nELAVL-regulated transcripts we describe here are numerous transcripts related to disease, including AD. For example, intronic nELAVL binding of the gene encoding the amyloid precursor protein, APP, was associated with skipping of exons 7 and 8. Both exons have previously been shown to be alternatively spliced and encode for the Kunitz protease inhibitory (KPI) motif, a domain that has been linked to APP processing. Remarkably, KPI domain containing isoforms of APP have been shown to be increased in AD, indicating that APP splicing might contribute to AD pathogenesis, and that nELAVL binding in human brain might be important to regulate the inclusion of the KPI domain. nELAVL regulates the splicing of two more AD-related transcripts, PICALM and BIN1, by promoting the inclusion of alternative exons 13 and 6a, respectively. Both proteins have been implicated in APP trafficking and both exons lie within domains mediating protein-protein interactions. Moreover, inclusion of the alternative exon 13 in PICALM has been linked to an AD-associated SNP, and we observed in this study that exon 6a of BIN1 shows a higher inclusion rate in controls compared to AD subjects. Since nELAVL binding promotes the inclusion of this exon, and control subjects show higher nELAVL binding, we propose that the altered splicing of BIN1 in AD subjects might be due to differential nELAVL binding. In fact, several nELAVL-regulated exons have been shown to be differentially spliced in AD subjects, further strengthening the link between nELAVL dysregulation and AD.

39While Y RNAs have not been linked to AD before, they have been implicated in various types of stress responses. The RNA binding protein RO60 usually associates with Y RNAs and is required for their stabilization. Besides RO60, Y RNPs contain several other RBPs such as ZBP1, MOV10, and Y-box proteins, and have been found to be remodeled upon stress. Our data suggests that nELAVL becomes increasingly associated with specific Y RNAs during both UV-induced stress and AD. ELAVL proteins can shuttle between nucleus and cytoplasm in response to environmental cues and preferentially accumulate in cytoplasmic stress granules upon cellular stress, and ELAVL binding to the CAT-1 transcript is modulated in response to stress in cultured cells. Interestingly, while we found that nELAVL specifically associates with Y RNAs during AD and acute UV stress, the nucleocytoplasmic distribution of nELAVL, RO60, and Y RNAs was not affected by UV stress. Because Y RNA levels remained constant, we propose that Y RNP complexes are specifically remodeled during AD and acute stress, which is not likely due to a change in nucleocytoplasmic protein/RNA distribution. These results are consistent with previous observations that stress induced shuttling might be limited to ELAVL1. Our observation of Y RNP remodeling in two very different systems of neuronal stress suggests that differential nELAVL/Y RNA association may be a widespread phenomenon and a focus of future studies.

40In addition to the four canonical human Y RNAs, hY1/3/4/5, hundreds of additional Y RNA genes are distributed throughout the human genome. The apparent lack of promoters upstream led to a premature designation of these Y RNAs as pseudogenes. Surprisingly, we found that hundreds of these Y RNA copies are expressed in human brain and neuroblastoma cells, although it remains unclear if these Y RNAs can still associate with RO60, because the RO60 binding site in many Y RNA copies is mutated. We observed that numerous Y RNA copies were more strongly associated with nELAVL in AD brain and acutely stressed cells, yet nELAVL binding did not affect their levels, indicating a function for this interaction other than Y RNA stabilization. While the outcome of nELAVL/Y RNA remains to be elucidated, our work revealed an aspect of nELAVL/Y RNA association related to stoichiometry. Hundreds of Y RNAs are bound by nELAVL in AD and UV-stress, which corresponds to up to 5% of all nELAVL CLIP tags. This shift of nELAVL binding may distort the normal stoichiometry of nELAVL interactions with its mRNA targets. Indeed, non-coding RNAs have previously been shown to affect RBP-RNA stoichiometry and therefore the biological function of other RNAs or RBPs. Our data indicate that the binding of nELAVL to Y RNAs during stress may lead to a redistribution of nELAVL binding and/or competition of nELAVL from other RNAs. Consistently, we found that high nELAVL/Y RNA association was associated with a general decrease in nELAVL binding at a subset of binding sites, especially within introns, and consequential splicing changes were reminiscent of splicing changes provoked by nELAVL depletion. Consistently, splicing changes induced by Y RNA overexpression showed an enrichment of nELAVL binding that was dependent on the presence of the ELAVL binding motif in Y RNAs. Hence we propose that the increased association of nELAVL and Y RNAs during stress causes sequestration of nELAVL from its mRNA targets.

41Taken together, our data indicate that nELAVL becomes strongly associated with Y RNAs in some AD subjects as well as in cells subjected to UV stress, and this is linked to a sequestration of nELAVL from some of its intronic targets, partially recapitulating splicing changes induced by nELAVL depletion. Our results are consistent with a hypothesis that a relatively subtle and perhaps long-term effect of Y RNA binding on normal nELAVL stoichiometry may underlie subtle and long-term changes in nELAVL biology. Perhaps analogously, the sequestration of the RBP, TDP-43, has previously been linked to neurodegenerative disorders. While the underlying mechanisms of TDP-43 and nELAVL sequestration are distinct, relatively subtle and long-term rearrangement of RNA:protein stoichiometry and interactions might be a recurrent theme of neurodegeneration.

42In conclusion, we have determined a robust and reproducible map of nELAVL binding sites on both mRNA and Y RNA targets in human brain. This database has both common and human-specific features that confirm and enhance previous work in mice and underscore its value as a resource for scientific inquiry. We have linked the data to genome-wide measurements of both mRNA levels and splicing suggesting specific functions for these important interactions. Moreover, we have uncovered a stress-modulated interaction of nELAVL proteins with non-coding Y RNAs. Y RNPs are remodeled during both UV-induced and chronic AD-related stress, which may be causally related to the pathophysiology of stress by causing a redistribution of nELAVL RNA target binding. Data implicating nELAVL proteins in several neurologic diseases, as well as the overlap of nELAVL binding sites we describe here with SNPs linked to the same and additional human diseases underscore the importance of this resource for the ongoing study of the molecular basis of human neurologic disease.

43Frozen brain tissue was obtained from the Mount Sinai Brain Bank. Subjects were classified according to CERAD criteria. Summaries of cognitive performance (assessed by clinical dementia rating, CDR), neurofibrillary tangles (assessed by Braak staging, BB) and plaque pathology (assessed by plaque counts) are shown in Supplementary file 1A. Brains were dissected along Brodmann areas and the two hemispheres were either stored as crushed powder at -80°C for biochemical analyses, or processed for stereological analyses, respectively. Eight control (no neurofibrillary tangles or plaque pathology), and nine advanced stage AD (CDR between 4 and 5) subjects matched for age and gender with short post mortem intervals (PMI) were selected. Crushed brain powder derived from the dorsolateral prefrontal cortex was either subjected to Trizol extraction or UV irradiation (3x 400 mJ/cm; see below) followed by CLIP analysis.

44Human HEK293T and neuroblastoma IMR-32 cells were obtained from ATCC (Manassas, VA) and maintained in DMEM (Fisher Scientific, Pittsburgh, PA), containing 10% (v/v) fetal bovine serum (Thermo Scientific, Waltham, MA), and 100 unit/ml penicillin/streptomycin (Life Technologies, Carlsbad, CA). 0.25% trypsin and 1% EDTA (Invitrogen) were used to passage cells every three days. Prior to UV treatment, cells were transferred to poly-D-lysine-coated 10-cm culture dishes (BD Biosciences, San Jose, CA) and allowed to adhere for at least 24 hr. Cells were washed in PBS (Invitrogen), layered with 5 ml PBS, and exposed to UVC (254 nm; 0.2 and 0.5 mJ/cm2) with 15 cm distance from the UV source using a Stratalinker. Cells were allowed to recover for 24 hr in fresh media before being harvested. For RNAseq, cells were washed in 5 ml PBS, and frozen in 1 ml Trizol (Invitrogen) at -80°C. For CLIP, cells were washed in 5 ml PBS and irradiated at 400 mJ/cm2 (see above). Cell pellets were collected by centrifugation at 2500 rpm for 3 min at 4°C and frozen at -80°C.

45For nELAVL RNAi experiments, cells were to transferred to poly-D-lysine-coated 6-well plates (BD Biosciences, San Jose, CA) and allowed to adhere for at least 24 hr. Accell Non-targeting pool (Cat#D-001910-10-20), and a mixture of Accell siRNA pools against Elavl2/3 and 4 (Cat#E019801-00-0010, #E011264-00-0010, #E016006-00-0010) from Dharmacon (LaFayette, CO) were added to Accell siRNA Delivery Media (Cat#B-001910-10), supplemented with 2% FBS for a final concentration of 1 µM per pool. Growth medium was removed and 1 ml Accell siRNA and media mixture was added. The Accell siRNA and media mixture was replaced with growth medium after 48h, cells were UV treated at 72 hr (see above, cells were layered with 1 ml PBS during UVC exposure), and cells were harvested at 96 hr for RNAseq analysis. Cells were washed in 1 ml PBS and frozen in 1 ml Trizol (75%) or lysed (25%) in 50 µl CLIP Wash Buffer 1 (see below).

46Canonical hY3 and hY3mut (sequences are shown below and mutated nucleotides are indicated) were cloned with AgeI and EcoRI into the transfer plasmid plKO.1 (; Addgene [Cambridge, MA] #10878). The puromycin selection cassette was replaced with GFP to monitor infection efficiency. HEK293T cells were transfected with transfer and packaging plasmids (of a three-plasmid lentivirus packaging system) at 80% confluency, using X-tremeGENE9 DNA transfection agent (Roche, Indianapolis, IN). Lentivirus containing supernatant was harvested after 24 hr, concentrated, aliquoted and stored at -80°C. IMR-32 cells were transduced with virus in the presence of polybrene (Sigma, St. Louis, MO) with supernatant of one 10 cm plate per 6 well (transductions were performed in triplicates). The infection efficiency was assessed after 48 hr (80% of cells were GFP-positive), cells were expanded and harvested after 72 hr for RNAseq analysis (one 6-well per 1 ml Trizol) and 96 hr for cell fractionation experiments.

47GCTGGTCCGAGTGCAGTGGTGTTTACAACTAATTGATCACAACCAGTTACAGATTTCTTTGTTCCTTCTCCACTCCCACTGCTTCACTTGACTAGCCTTTT

48GCTGGTCCGAGTGCAGTGGTGTCTACAACTAATTGATCACAACCAGTTACAGATCTCTCCGTTCCTTCTCCACTCCCACTGCTTCACTTGACTAGCCTTTT

49Cells grown in 10-cm culture dishes were scraped in 1 ml PBS after media removal and spun at 500 g for 5 min at 4°C. Cells were lysed in 500 µl buffer A (150 ml NaCl; 50 mM Tris, pH 7.4; 0.01% Saponin; 1x Protease Inhibitor Cocktail [Thermo Fisher Scientific, Waltham, MA]), and incubated on ice for 10 min. Cytoplasm and nuclei were separated with a 10 min spin at 3000 g at 4°C, and both fractions were washed in 500 µl buffer A before an additional 10 min spin at 3000 g at 4°C. Nuclei were resuspended in 500 µl buffer A and sonicated 3x for 5 s. RNA was Trizol (Life Technologies) extracted from 100 µl lysate and reverse transcribed using iScript (Bio-Rad, Hercules, CA). qPCR was performed with FastStart SYBR Green Master (Roche, Indianapolis, IN), requiring at least one primer in each mRNA primer pair to be specific for an exon junction. qPCR results were normalized as indicated. Western blot analysis was performed using 15 µl cell lysate per lane, and the following antibodies: α-GAPDH (1:25:000; Abcam [Cambridge, MA] ab8245), Hu subject antiserum (1:1000 dilution), α-HSP90 (1:1000; Cell Signaling Technology [Danvers, MA] 4877S), α-H3 (1: 2000; Abcam ab1791), α-RNA PolII (1:1000; Millipore [Billerica, MA] 05–623), and α-RO60 (1:50; gift from S. Wolin). Immunoreactive bands were analyzed with the Odyssey Infrared Imaging System (LI-COR, Lincoln, NE) and normalized as indicated.

50nELAVL HITS CLIP was performed as previously described with the following modifications. nELAVL-RNP complexes were immunoprecipitated using paraneoplastic Hu-antiserum, which recognizes all three neuronal ELAVL isoforms. 80 mg human brain powder or one 10-cm 70% confluent plate of IMR-32 cells were used per immunoprecipitation. For controls (no UV irradiation, control serum, and overdigested control), 40 mg of human brain powder or half a mouse brain were used, respectively. Prior to phosphatase treatment, beads were washed for 3 min each in a series of wash buffers (see below). nELAVL bound RNA fragments from two IMR-32 samples were ligated to an indexed degenerate 5’ RNA linker (see Supplementary file 1A); RNA fragments from the remaining samples were ligated to a degenerate 5’ RNA linker. The two-step PCR amplification was performed with Accuprime Pfx (Invitrogen). cDNA from IMR-32 and subject samples was initially amplified with DP3/5. Brain cDNA from subjects 1–3, and 9–11 was then amplified with DSFP3/5, samples 4–8, and 12–17 were amplified with MSFP3 and indexed MSFP5. Brain samples were sequenced on an Illumina (San Diego, CA) GAIIX at the Rockefeller University Genomics Resource Center with the standard Illumina primer (samples 4–8, and 12–17) or the custom primer SSP1 (samples 1–3, and 9–11). IMR-32 cDNA was either amplified with SP3/5-PE (samples with indexed 5’linker), and sequenced on the MiSeq system with the custom primer SSP1, or with MSFP3 and indexed MSFP5 and sequenced on the MiSeq system with the standard Illumina primer. See Supplementary file 1A for used indexes.

51Wash buffer 1: 1xPBS; 0.1%SDS; 0.5%Na-DOC; 0.5% NP-40

52Wash buffer 2: 5xPBS; 0.1%SDS; 0.5%NA-DOC; 0.5% NP-40

53Wash buffer 3: 15 mM Tris, pH7.4; 5 mM EDTA; 2.5 mM EGTA; 1% Triton X-100; 1% Na-DOC; 0.1% SDS; 120 mM NaCl; 25 mM KCl

54Wash buffer 4: 15 mM Tris, pH7.4; 5 mM EDTA; 2.5 mM EGTA; 1% Triton X-100; 1% Na-DOC; 0.1% SDS; 1 mM NaCl

55Wash buffer 5: 15 mM Tris, pH7.4; 5 mM EDTA

56Wash buffer 6: 50 mM Tris, pH7.4; 10 mM MgCl2; 0.5% NP-40

573’ RNA linker: 5’P – GUG UCA GUC ACU UCC AGC GG 3’ – puromycin

585’ RNA linker: 5’OH – AGG GAG GAC GAU GCG GNN NNG 3’ – OH

595’ RNA linker indexed: 5’OH – AGG GAG GAC GAU GCG GXX NNN NG 3’ – OH

60DP5: 5’ – AGG GAG GAC GAT GCG G – 3’

61DP3: 5’ – CCG CTG GAA GTG ACT GAC AC – 3’

62DSFP5: AATGATACGGCGACCACCGACTATGGATACTTAGTCAGGGAGGACGATGCGG

63DSFP3: CAAGCAGAAGACGGCATACGACCGCTGGAAGTGACTGACAC

64SP5-PE: AATGATACGGCGACCACCGAGATCTACACCTATGGATACTTAGTCAGGGAGGACGATGCGG

65SP3-PE:CAAGCAGAAGACGGCATACGAGATCTCGGCATTCCTGCCGCTGGAAGTGACTGACAC

66MSFP3: CAAGCAGAAGACGGCATACGAGATCCGCTGGAAGTGACTGACAC

67MSFP5: AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCTXXXXAGGGAGGACGATGCGG

68SSP1: CTATGGATACTTAGTCAGGGAGGACGATGCGG

69Analyses were carried out using the Galaxy suite of bioinformatics tools (http://main.g2.bx.psu.edu/), in addition to publicly available in-house tools. Data was visualized with UCSC genome browser (http://genome.ucsc.edu/). To reduce mis-alignments due to sequencing errors, reads were initially filtered based on quality score (≥20 in the degenerate linker region; average of ≥20 in the remaining read). Exact sequences were collapsed to remove PCR duplicates. The degenerate barcode (and index if applicable) were removed and the 3’ linker was trimmed. Using FASTA files as input, reads were subsequently mapped to the hg18 build of the human genome by novoalign (www.novocraft.com), requiring unambiguous mapping and a maximum of two mismatches. To identify unique CLIP tags, we applied stringent filtering, and collapsed reads with the same starting genomic positions, and only unique tags (Supplementary file 1A) were used for subsequent analyses. Peaks with significant nELAVL binding compared to background (p-value<0.01) were identified utilizing a similar approach from previous studies. Specifically, we used a scan statistics to compute p-values in which each observed PeakHeight (PH, CLIP tags within a peak) was compared to the PH one would expect by random chance, assuming a background of uniformly distributed CLIP tags in each gene. Bonferroni correction was applied to the resulting p-values to correct for multiple hypothesis testing. All scripts used in the analysis including the peak finding algorithm and more information can be publically obtained at http://zhanglab.c2b2.columbia.edu/index.php/Resources.

70RNA from human brain and IMR-32 neuroblastoma cell lines was Trizol (Invitrogen)-extracted, Ribo-Zero-selected (Epicentre, Madison, WI), DNase-treated (Roche), and prepared for sequencing, following the Illumina High-throughput TruSeq RNA Sample Preparation Guidelines. The libraries from subjects 6–8 and 15–17, as well Y RNA overexpression samples, were sequenced on an Illumina HiSeq 2500 at the New York Genome Center, generating 125-bp paired end reads. The remaining libraries were sequenced on an Illumina HiSeq 2000 at the Rockefeller University Genomics Resource Center, yielding 100-bp paired end reads. Reads were aligned to the hg18 build of the human genome using TopHat, allowing a maximum of two mismatches. Only unambiguously mapped reads were kept for analysis. Additional exon junctions not observed due to the gap between paired end mates were inferred using a Bayesian method, and a set of non-redundant constitutive exons with relative high inclusion rate (according to ESTs) was used to estimate gene expression. Exon and exon junction reads were inferred as previously described. 8163 cassette exons in brain, 9629 cassette exons in the RNAi/UV IMR-32 samples, and 10,432 cassette exons in the Y RNA infected IMR-32 samples remained after filtering on exon junction coverage to reduce multiple testing (coverage was normalized for library size; normalized coverage of each isoform and each condition ≥5).

71CLIP and RNAseq data were deposited in the GEO database with accession number GSE53699. Analyses of unique CLIP tags and RNAseq data were carried out in R (www.r-project.org;) using Bioconductor and the packages Biostrings, edgeR, GenomicRanges, ggplot2, Hmisc, plyr, qvalue, reshape, scales, and VennDiagram. mRNA abundance and nELAVL peak binding (PeakHeight, PH) changes were assessed by differential analysis of raw sequencing counts in edgeR using the TMM methodology for normalization and a negative binomial generalized linear model. Only expressed genes and robustly bound peaks were analyzed to reduce multiple testing. Expressed genes were defined as genes that had more than one count per million (cpm) in at least 5 control brain samples (out of 8), in 6 AD brain samples (out of 9), or in 4 IMR-32 samples (out of 12 RNAi/UV samples or 9 Y RNA infection samples, respectively). Robustly bound peaks were defined as peaks that had more than one cpm in at least 5 control or AD brain samples (out of 8 and 9, respectively), or in 2 control or 4 UV IMR-32 samples (out of 3 and 6, respectively). We controlled for batch effects of brain samples (see Supplementary file 1A) and estimated the false discovery rate (FDR) with an optimized FDR approach (q-value methodology) to correct for multiple hypothesis testing. PeakHeight (PH, CLIP tags within peaks) was normalized for library size, and, to account for differences in gene expression level, normalized PH was determined by dividing PH by rpkm of the corresponding gene. nELAVL binding was defined as the summary of PH per gene, whereas normalized nELAVL binding was defined as the summary of normalized PH per gene. We similarly defined 3’UTR and intronic binding by summarizing only 3’UTR or intronic peaks, respectively. Top targets were defined as the top 1000 genes according to normalized nELAVL binding. mRNA abundance was defined as cpm RNAseq tags per gene. We did not normalize nELAVL binding nor mRNA abundance for transcript length as we performed differential analysis of the datasets. A pseudocount of 1 was added before log2 transformation. Cross-correlation plots show log2 raw counts of PH or RNAseq reads per gene. Splicing changes were determined using the observed inclusion and exclusion junction read counts and by fitting a generalized linear model with a logit link function. For each exon, GLM likelihood ratio test was conducted to test if there was a significant difference in the fraction of exon inclusion (delta Inclusion, ΔI) between conditions, and controlling for batch variables of brain samples. ΔI (equivalent to delta PSI) was calculated by subtracting fraction exon inclusion of Sample 2 from fraction exon inclusion of Sample 1. FDR was calculated as described before using the q-value method. High confidence alternative splicing events were identified by requiring a stringent criteria of FDR<0.05 and ΔI≥0.1. P-values to assess the statistical significance of the overlap between gene lists were calculated with a hypergeometric test (Figure 2A/9D, Figure 9—figure supplement 5). A paired Wilcoxon rank sum test was used to evaluate differences in Y RNA binding (Figure 6C/7B), and a t-test was used to evaluate changes in mRNA/Y RNA/protein abundance (Figure 3B/9A,B, Figure 3—figure supplement 1, Figure 6—figure supplement 1, Figure 7—figure supplement 1, 3,Figure 8—figure supplement 1A, Figure 9—figure supplement 1, 2). Differences in cumulative density curves were evaluated with a one-sided KS (Kolmogorov-Smirnov) test (Figure 3B). P-values to assess motif enrichment were calculated with a Fisher’s exact test (Figure 6B/7A).

72Y RNAs were derived from GENCODE (v19 release). nELAVL binding on Y RNAs was defined as CLIP tags on Y RNAs. To discriminate between individual Y RNA copies we used CLIP reads, aligned with novoalign. While this mapping strategy allowed us to discrimate between individual Y RNA copies with high confidence, CLIP reads that mapped to multiple Y RNAs were discarded. Exclusively for Figure 7—figure supplement 4, CLIP tags were therefore additionally aligned with Bowtie2, allowing multiple alignments and reporting one. This allowed us to estimate the overall amount of CLIP tags on all Y RNAs (shown in Figure 7—figure supplement 4). Y RNA abundance in brain, and UV- and RNAi-treated IMR-32 cells was defined as RNAseq tags on Y RNAs, using unambiguously mapped tags aligned with TopHat. Similar to above, reads that mapped to multiple Y RNAs were discarded using this strategy. To estimate overall Y RNA overexpression upon lentiviral infection, we mapped RNAseq reads of infected samples with Bowtie2, allowing multiple alignments and reporting the best matched alignment (Figure 9B). Y RNA overexpression was determined with two additional strategies. We performed qPCR with primers that recognized infected Y RNAs (Y3wt and Y3mut) but also an unknown number of endogenous Y RNAs (Figure 9A). Additionally, we directly searched for 40- and 68-nt sequences of Y3wt and Y3mut sequences in the Illumina sequenced library reads (both 40- and 68-nt sequences spanned the mutated residues). The 68-nt sequence was only present in the infected Y RNAs and the canonical hY3 RNA, and searching for this sequencing length showed that both Y3wt as well as Y3mut were overexpressed. To more accurately estimate the extent of overexpression, we additionally searched for a 40-nt Y RNA sequence. The Y3wt sequence is present in a large number of endogenous Y RNAs, whereas the Y3mut sequence only detects infected Y3mut. By comparing Y3wt and Y3mut reads in Y3mut infected samples, we observed that infected Y RNA reads (Y3mut reads) correspond to 10% of total Y3 RNA reads (Y3wt reads).

73The MEME-ChIP Suite and HOMER were used for motif discovery. Coordinates of the top 500 peaks in healthy brain +/- 25 nt were used as input. The given strand was searched for 6–10-nt-long motifs for up to 10 motifs, allowing any number of repetitions with MEME-ChIP. We additionally searched the given strand with HOMER v4.7 using default parameters and 50 nt flanking sequences of the 500 peaks as background sequences. The respectively top enriched motif is shown.

74The 1000 nELAVL top target gene products were seeded in a large protein-protein interaction (PPI) network containing 14,191 genes and more than 197,000 interactions compiled from various online databases: BioGrid, MINT, KEGG, PPID, HPRD, DIP, BIND, IntAct, InnateDB, and SNAVI. To reduce the false positive interactions in the background network, we excluded interactions reported using high-throughput approaches. Direct interactions with other seed genes were kept for each seed gene, and a connected subnetwork was created using these seed genes as nodes. Organic clustering of the obtained subnetwork was performed using the network visualization software yEd (http://www.yworks.com/en/products_yed_about.html). In a second network, six AD genes (APP, BACE1, MAPT, PICALM, PSEN1 and PSEN2) were concatenated with the input seed gene lists to investigate their relationship with nELAVL target genes. The composition and the clustering of the PPI networks with and without those six genes were identical. Only the network including those additional seeds is shown.

75The clusters of nELAVL target genes in the PPI subnetwork represent protein clusters based on direct physical interactions. Gene products in clusters with at least 10 nodes were examined for an enrichment of functional annotations with Enrichr using the Fisher’s exact test and computing an adjusted p-value for multiple hypothesis testing using the Benjamini-Hochberg correction. The six additional AD associated genes used for the subnetwork clustering were not taken into account for the enrichment analysis.