Concentration-dependent Splicing through Suboptimal Motifs Enables Waves of Gene Regulation in Neuronal Development by Bridget E. Begg B.A., Biological Sciences and English (2013) Wellesley College SUBMITTED TO THE DEPARTMENT OF BIOLOGY IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY AT THE MASSACHUSETTS INSTITUTE OF TECHNOLOGY FEBRURARY 2022 © 2022 Massachusetts Institute of Technology All rights reserved Signature of author…………………………………………………………………………………. Bridget E. Begg Department of Biology December 10, 2021 Certified by………………………………………………………………………………………… Christopher B. Burge Professor of Biology Thesis Supervisor Accepted by………………………………………………………………………………………... Amy E. Keating Professor of Biology and Biological Engineering Co-Director, Biology Graduate Committee 1 2 Concentration-dependent Splicing through Suboptimal Motifs Enables Waves of Gene Regulation in Neuronal Development by Bridget E. Begg Submitted to the Department of Biology on December 10, 2021 In Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy Abstract Alternative splicing, which occurs in over 95% of human genes, is the process by which exons are differentially included in transcripts produced from the same gene to produce a variety of transcript isoforms. This mode of post-transcriptional gene regulation endows the 20,000 protein-coding genes in the human genome with tunable expression, manifold protein–protein interactions, and an additional layer by which to regulate protein localization, activity, and signal response. Increasingly, alternative splicing is understood to be a significant contributor to organismal complexity in animals. The Rbfox family of splicing factors regulates alternative splicing during animal development and in disease, impacting thousands of exons in the maturing brain, heart, and muscle. Although it is well established that Rbfox binds to the RNA sequence GCAUG with high affinity and specificity, this motif is responsible for only half of the Rbfox binding sites observed in cellular and neuronal contexts. We incubated recombinant RBFOX2 with over 60,000 mouse and human transcriptomic sequences to reveal substantial binding to several moderate-affinity, non-GCAYG sites at a physiologically relevant range of RBFOX2 concentrations. We find that these “secondary motifs” bind Rbfox robustly in cells and that several together can exert regulation comparable to GCAUG in a trichromatic splicing reporter assay. In the brain, secondary motifs regulate RNA splicing in neuronal development, enabling a second wave of splicing changes as Rbfox levels increase. Additionally, secondary motifs are activated in neuronal subtypes according to cellular gene expression levels of Rbfox, contributing to tissue diversity in the mammalian brain. This work presents the first observation of spatiotemporal regulation through suboptimal motifs in RNA-binding proteins, a phenomenon that may be a widespread. Furthermore, the characterization of Rbfox secondary motifs reveals new regulatory targets of an essential splicing factor, which may contribute to mammalian brain development, synaptic plasticity, and human disease. Thesis Supervisor: Christopher B. Burge Title: Professor 3 4 Table of Contents Abstract………………………………………………………………………………………..…3 Chapter 1. Introduction…………………………………………………………………………7 Alternative splicing in eukaryotes……………………………………………………………….7 Regulation of alternative splicing by RNA-binding proteins…………………………………..10 Diverse molecular modes of RBP–RNA recognition………………..…………………………15 Quantitative, high-throughput biochemistry to describe specificity and affinity………………20 Suboptimal motifs in gene regulation…………………………………………………………..24 Chapter 2. Concentration-dependent splicing is enabled by Rbfox motifs of intermediate affinity……………………………………………………………….29 Natural sequence RBNS recovers known features of RBFOX2 binding………………………33 RBFOX2 binds a set of secondary motifs of intermediate affinity in vitro……….…………....36 Biochemical characterization of Rbfox secondary motifs……………………………………...38 RBFOX2 binds to specific secondary motifs in vivo…………………………………………..42 Secondary motifs regulate splicing in an Rbfox-dependent manner…………………………...46 Secondary motifs regulate splicing at specific stages of neuronal differentiation……………...49 Rbfox expression-dependent splicing through secondary motifs across neuronal subtypes…...54 A quantitative model for Rbfox expression-dependent, differential motif activity…………….57 Insights into concentration-dependent regulation via RBP suboptimal motifs………………...58 Acknowledgements…………………………………………………………………….…….....62 Materials and methods………………………………………………………………………….63 Supplementary figures………………………………………………………………………….71 Chapter 3. Future Directions.……………………………………………………...……….….86 The biological impact of Rbfox secondary motifs……………………………………………...86 Suboptimal motifs in neuronal development…………………………………………………...89 References……………………………………………………………………………………….92 Appendix A. Exon-mediated activation of transcription starts…………………………….106 Appendix B. Loss of LUC7L2 and U1 snRNP subunits shifts energy metabolism from glycolysis to OXPHOS………………………………………………………….139 5 6 Chapter 1. Introduction Alternative splicing in animals Alternative splicing of precursor messenger RNAs (mRNAs) greatly expands the coding capacity of the genome, contributing to variation in mRNA expression and affording the proteome incredible diversity (Ellis et al., 2012; Nilsen & Graveley, 2010; Yang et al., 2016). In humans, over 95% of genes are alternatively spliced, enriching the information encoded across the entire genome (Wang et al., 2008). Indeed, alternative splicing is correlated with measures of organism complexity across metazoans, and likely contributes significantly to the phenotypic innovations of mammals (Bush et al., 2017; Kim et al., 2007; Yang et al., 2021). Such innovation is evident in the abundant utilization of complex alternative splicing modes in the mammalian brain. For instance: Neurexins are cell adhesion proteins in presynaptic neurons that modulate signaling at the neuronal synapse by binding to the postsynaptic protein neuroligin (Scheiffele et al., 2000). In mammals, neurexins are produced by three different genes, each with two alternative promoters and alternative splicing at five different exons. Combinatorally, neurexins can thus express up to two thousand different isoforms, probably contributing to the remarkable synaptic specificity required for synaptic plasticity and contributing to human processes as profound as learning and memory (Gomez et al., 2021). In mammalian cells, splicing proceeds through cis-sequence elements in nascent and pre- mRNAs by recruiting a modular, macromolecular RNA–protein complex referred to as the spliceosome. More than 150 proteins are associated with the mammalian spliceosome, generally functioning to stabilize the catalytic steps mediated by five small nuclear RNAs (snRNAs), U1, U2, U4/U6, and U5, over the course of the splicing reaction. Each snRNA associates with a defined cohort of requisite protein cofactors, together with the snRNA referred to as small 7 nuclear ribonucleoproteins (snRNPs), and snRNPs associate and dissociate with the spliceosome at distinct molecular steps during splicing. However, some snRNP-associated proteins are not obviously required for successful splicing and may facilitate alternative splicing under specific cellular conditions (Jourdain et al., 2021; Papasaikas et al., 2015; Ule & Blencowe, 2019). The U1 snRNP initiates splicing through base pairing of the U1 snRNA with a 5′ splice site (5′ ss) cis-sequence at an exon–intron boundary within the pre-mRNA. U2 snRNP and U2 auxiliary factor (U2AF) also recognize the 3′ splice site (3′ ss) and the polypyrimidine tract (PPT) through cis-sequence elements at a neighboring intron–exon boundary, and physical interactions between these two snRNPs ultimately define the exon to be spliced into or the intron to be removed from the RNA. Following this definition, the U2 snRNP relocates to the branch point sequence (BPS) upstream of the PPT, and subsequent recruitment of U4/U6 and U5 forms a tri-snRNP complex that executes the catalytic steps of splicing: two trans-esterification reactions that enable lariat formation, intron removal, and exon ligation. The biochemical steps of mammalian splicing have been clearly explicated by recent cryo-EM studies of yeast and mammalian spliceosomes (reviewed in Zhang et al., 2019). These structures explained key features of the splicing reaction, including that it proceeds through an RNA-based catalytic mechanism, confirming that the spliceosome is a ribozyme. RNA cis-sequence elements direct both constitutive and alternative RNA splicing. Since the sequencing of the human genome, much work has been done to disentangle the “splicing code” to completely predict RNA splicing given only genomic sequence. Initial efforts relied on the four absolutely required sequence elements that determine exon/intron recognition and splicing for both constitutive and alternative events: the 5′ ss, 3′ ss, BPS, and PPT. Early computational analyses found that intronic sequence elements meaningfully contribute to more accurate 8 splicing predictions, especially in mammals (Lim & Burge, 2001; Yeo & Burge, 2004). High- throughput splicing reporter technology enabled the screening of libraries of diverse genomic or random sequences for putative splicing regulatory elements (SREs) acting as exonic splicing enhancers (ESEs), exonic splicing silencers (ESSs), intronic splicing enhancers (ISEs), or intronic splicing silencers (ISSs) (Fairbrother & Chasin, 2000; Hiller et al., 2007; Wang et al., 2004, 2006; Yeo & Burge, 2004; Yeo et al., 2007). Many of these SREs can, in fact, be bound by multiple distinct splicing factors individually, cooperatively, or competitively, creating context- dependent and tissue-specific splicing determined by the milieu of nuclear RBPs (Fu & Ares, 2014; Wang et al., 2012). Furthermore, studies of individual RBPs would show that the position of an SRE also frequently determines the inclusion or exclusion activity of a splicing factor (Ule et al., 2006; Van Nostrand et al., 2020). These results provide a molecular rationale for why the same RNA sequence element may have opposing effects on splicing, while also highlighting that a major challenge that persists in the quantitative prediction of splicing is its poorly understood dependence on the expression level and binding preferences of hundreds of RBPs. To this end, the application of machine learning methods to the prediction of splicing outcomes through cis- sequence elements has reinvigorated efforts to determine the “splicing code,” including predicting the effects of human variants on splicing and disease outcomes (Barash et al., 2010; Xiong et al., 2015). Most recently, a neural network architecture confirmed that most SREs operate locally, as most splicing could be accurately predicted from sequences within 400 nucleotides (nt) of an exon (Jaganathan et al., 2019). However, neural networks are inherently limited in the biological information that they can infer, and tissue-specific splicing events were often scored ambiguously by their model, demonstrating the continued importance of 9 considering the cellular RBP environment to the full understanding of alternative splicing. Regulation of alternative splicing by RNA-binding proteins The regulation of alternative splicing proceeds largely through RBPs recruited by cis- sequence RNA elements that appear in an alternative exon and its neighboring introns. In addition to spliceosomal RBPs, there are other distinct cohorts of pre-mRNA-associating proteins that function to determine splice site selection early in the splicing reaction by directly facilitating or obstructing the interactions of U1 and U2 snRNP with the pre-mRNA (Mayrand et al., 1986; Wang et al., 1995). Hundreds of proteins have now been implicated in the regulation of splicing (Gerstberger et al., 2014; Ule & Blencowe, 2019), but nuclear RBPs were first identified in rat liver cells and subsequently characterized in the lampbrush chromosome of newt oocytes, where they were distinguished from proteins binding cytoplasmic messenger or ribosomal RNA (Moulé & Chauveau, 1968; Sommerville, 1973). Hypothesized to be proteins associated with the primary transcripts of eukaryotic RNA, nuclei from 3H-uridine-labeled human cells were separated by velocity centrifugation, and the nuclear ribonucleoproteins (nRNPs) were found to exist in two classes: a few stable species in large amounts, possibly associated with the chromatin–RNA fraction, and a large number of smaller species that are rapidly turned over, appearing to associate with RNA after its release from chromatin (Augenlicht & Lipkin, 1976). This was the first recognition of two types of splicing factors: (1) ubiquitous, often constitutive factors like the spliceosomal proteins, Serine/arginine-rich (SR) proteins, and heterogeneous nuclear ribonucleoproteins (hnRNPs), and (2) tissue- and context-specific alternative splicing factors, including highly-studied RBP families such as Nova and Rbfox. 10 SR proteins and hnRNPs contribute to both constitutive and alternative mammalian RNA splicing, with both classes of RBPs comprised of a structurally diverse cohort of splicing factors that appear to have co-evolved with the complexity of alternative splicing in eukaryotes (Busch & Hertel, 2012). SR proteins are typically comprised of one or two RNA recognition motifs (RRMs) and an arginine/serine (RS) repeat domain of variable length; the RS domain interacts directly with the spliceosome and its length is positively correlated with its activity. SR proteins are generally understood to enhance spliceosome recognition of the splice site, resulting in greater exon inclusion. It has now been observed that the effect of an individual SR protein is dependent not only on its binding location, but also on that of other SR proteins in the vicinity of the putative splicing event (Pandit et al., 2013; Shen & Mattox, 2012). hnRNPs, a larger and more diverse class compared to SR proteins, are structurally comprised of modular RRMs and K-homology (KH) RNA-binding domains that recognize SREs as well as variable arginine/glycine/glycine (RGG) boxes and other repeat domains that engage in protein–protein interactions to execute their functions (Busch & Hertel, 2012). Opposite the function of SR proteins, hnRNP proteins are broadly considered to be splicing silencers, but in fact their ultimate effects on splicing have also shown a strong positional dependence relative to the exon. It is thought that increased functional diversification in hnRNP proteins relative to SR proteins may have evolved due to the degeneration of splice-site specificity in mammals that requires precise distinguishing of alternative exons from abundant pseudoexons that appear throughout the genome (Busch & Hertel, 2012). Tissue-specific splicing is most conserved in brain, heart, and muscle, and likely contributes to the highly specialized cells within those organs (Merkin et al., 2012; Barbosa-Morais et al., 2012). Furthermore, alternative splicing can be specific to a given cellular context, determined 11 by differential RBP or cofactor expression, substrate levels, or subcellular localization varying across cell state, cell type, and in response to stimuli. The varying expression of the context- and tissue-specific RBPs in the cell provides the basis for the regulation of alternative splicing (Fu & Ares, 2014; Gerstberger et al., 2014); it is thus critical to understand their individual functions and specificities. Nova was the first mammalian tissue-specific splicing factor discovered, and many common properties of RBPs are reflected in its biology. Nova-1 and Nova-2 were identified as antigens in sera for a motor disorder, paraneoplastic opsoclonus myoclonus ataxia (POMA), and Nova-1 was immediately recognized as potential splicing factor due to its homology with hnRNP K and the yeast splicing factor MEF1 (Buckanovich et al., 1993, 1996; Yang et al., 1998). The specificities of these neuronal splicing factors were found to differ only slightly, with both NOVA1 and NOVA2 binding the core motif YCAY while preferring [UCAU(N)0–2]3 and GAGUCAU, respectively, though both paralogs demonstrate affinity for both motifs (Buckanovich & Darnell, 1997; Jensen et al., 2000; Yang et al., 1998). Nova was the first protein analyzed by protein– RNA crosslinking and immunoprecipitation (CLIP) (paired with microarray analysis), which found that it binds a biologically-coherent set of transcripts in the mouse brain, establishing a framework for conceiving of RBP–RNA interactions as networks of gene regulation (Ule et al., 2003). Subsequent CLIP studies focusing on the core Nova YCAY motif produced “RNA maps” relating Nova binding with its effects on splicing (Ule et al., 2005, 2006). Specifically, they found that binding sites within or just upstream of the alternative exon tended to induce silencing, and binding sites that enhanced inclusion were found in the downstream and upstream introns. For both enhancement and silencing of splicing, clustering of Nova binding sites is frequently observed. 12 Another family of context-specific splicing factors essential for mammalian brain development, Rbfox was first identified as the sex-determining genomic locus feminizing on X (fox-1) in Caenorhabditis elegans (Hodgkin et al., 1994). In mammals, the Rbfox family consists of three paralogs highly expressed in heart, skeletal muscle, and brain tissues. The observation that peak Rbfox expression coincides with late neuronal development suggests that this protein is important for proper brain function, and is supported by the finding that seizures and defects in cerebellar function occur upon targeted central nervous system (CNS) knockout in mice (Gehman et al., 2011, 2012). Mutations in RBFOX1 and RBFOX3 in humans are associated with autism and epilepsy, and Rbfox triple-knockout mouse embryonic stem cells (mESCs) fail to develop into mature ventral spinal neurons (Jacko et al. 2018). To this end, studies of Rbfox in vivo binding and regulation using CLIP have identified thousands of putative regulatory targets responsible for its critical role in development (Jacko et al., 2018; Weyn-Vanhentenryck et al., 2014; Zhang et al., 2008). Furthermore, studies of Rbfox binding in vitro using surface plasmon resonance (Stoltz, 2015) and high-throughput affinity assays have shown that Rbfox proteins primarily recognize the RNA sequence element (U)GCAYG (Jin et al., 2003; Lambert et al., 2014). The identification of this motif provided an important link between the molecular binding preferences of this RBP and its cellular role (Jin et al., 2003; Minovitsky et al., 2005), as (U)GCAUG had already been discovered as a sequence motif enriched in 25 brain-specific alternative exons and their neighboring introns (Brudno, 2001), and Rbfox proteins have since been shown to regulate numerous alternative exons through this motif (Minovitsky et al., 2005; Sun et al., 2012; Weyn-Vanhentenryck et al., 2014; Yeo et al., 2009). However, only about half of CLIP peaks can be accounted for with a GCAYG (Yeo et al., 2009), suggesting the existence 13 of additional determinants of Rbfox binding and its corresponding gene regulatory outcomes (Lambert et al., 2014; Lee et al., 2016). An alternative to single-gene studies of splicing factors, high-throughput discovery of non- spliceosomal RBPs have been facilitated by advances in CLIP and protein identification by mass spectrometry. Early work to characterize the RNA-bound proteome of mammalian cells identified ~850 RBPs, 200–300 more than previously identified by low-throughput and computational methods (Anantharaman et al., 2002; Baltz et al., 2012; Butter et al., 2009; Castello et al., 2012). A complete, manually-curated estimate of human RBPs established a count of 1,542 RBPs, accounting for 7.5% of all genes and ~20% of the expressed proteome (Gerstberger et al., 2014). While the majority of RBPs (~70%) are ubiquitously expressed, up to one third of human RBP families have members that are tissue-specific. More recent efforts to catalogue RBPs have focused on the detection of RBPs without classical RNA-binding domains and have proposed RNA-binding functions for metabolic enzymes (Beckmann et al., 2015; Perez-Perri et al., 2018). This result, as well as the finding that a U1 snRNP component is responsive to glucose starvation, highlight the emerging connection between post-transcriptional gene regulation and energy metabolism (Christofk et al., 2008; Gao et al., 2011; Jourdain et al., 2021). In addition to accounting for unconventional RBPs, it will likely be important for future studies to consider dynamic expression and isoform-level changes in the RBP proteome, as well as post-translational modifications of RBPs and post-transcriptional modifications of their RNA substrates, as each of these features has been shown to modulate the quantitative binding patterns, and therefore transcriptome-wide regulatory activity, of individual RBPs (Hentze et al., 2018). 14 Diverse molecular modes of RBP–RNA recognition Virtually every step in the life cycle of an mRNA—transcription, splicing, nuclear export, modification, translation, and degradation—is established though direct RBP binding to RNA sequences. The temporal and spatial specificity with which individual RBPs bind suggests that a diverse repertoire of sequence and structural determinants are operating within both RBPs and their cognate RNAs. In principle, single-gene structural studies of RBP–RNA recognition should clarify molecular determinants of RBP specificity; however, as high-resolution structures of RBPs in complex with RNA are developed, it becomes clearer that these structures are not generalizable, and the answers to how most RBPs discriminate among millions of cellular RNAs are unknown or incomplete. Most RBPs contain multiple globular RNA binding domains (RBDs) interspersed with unstructured “linker” amino acid sequences to form a highly modular structure (Lunde et al., 2007). Individual RBDs typically bind 3–6 nucleotides (nt), and RBPs predominately bind single-stranded RNA, though the tolerance for RNA secondary structure at individual binding sites is variable, typically increasing in proportion to motif affinity (Dominguez, et al., 2018; Jens et al. in preparation). Catalytically active RBPs contain a variety of enzymatic RBDs mediating the activity of the protein—this group includes nucleotidyltransferases, ribonucleases, RNA-modifying enzymes, GTPases, and helicases. However, more than half of mRNA-binding RBPs consist of RBDs that participate near- exclusively in RNA–protein interactions without an enzymatic function: the most abundant globular domains of this type include the aforementioned RRMs and KH domains as well as Zinc-finger (ZF) domains (Gerstberger et al., 2014). Despite the abundance and order of globular RBDs, it is difficult to infer RBP specificity from the aa sequences of RBDs due to the incredible range of their recognition mechanisms. 15 Structures of RBPs bound to their cognate motifs have revealed an abundance of non-canonical aa interactions, altered tertiary structures within individual RBDs, combinatorial specificity and recognition for multidomain RBPs, and unpredictable roles for disordered regions in binding. For instance, binding by an ~90 aa RRM is canonically configured by aromatic side-chains of the beta-sheet, often via two short (6–8 aa), conserved amino-acid stretches called RNP1 and RNP2. Despite this stated consistency, NMR and crystallographic structures have shown that RRMs recognize RNA in a variety of ways, with differential contributions of each RNP motif, β-sheet extension, and variation in the protein “loops” between individual alpha-helices and β-strands (Cléry et al., 2008). An NMR structure of RBFOX1 bound to its cognate RNA UGCAUGU revealed non-canonical RRM binding via a phenylalanine in the β-1/α-1 loop to coordinate binding to the UGC and induce base-pairing of the G2 and A4 resides; this non-canonical mode of recognition results in exceptionally tight binding to the RNA for a single RRM (Auweter et al., 2006) (Figure 1). KH domains, which have a core protein secondary structure of β-α-α-β also frequently deviate from their “canonical” binding modes. All known KH domains interact with the RNA backbone through a GXXG loop (where X denotes any aa) between α-1 and α-2. Canonically, this positions four nucleobases of the motif YMMY (where Y denotes C or U and M denotes A or C) into a hydrophobic groove in the protein that achieves RNA specificity through both hydrophobic interactions and hydrogen-bonding. The KH1 and KH3 domains of Nova-1 recognize RNA of the motif YCAY in this manner (Lewis et al., 1999; Nicastro et al., 2015; Teplova et al., 2011). Thus, all KH domains were initially thought to be constrained to the recognition of the core dinucleotide MM, but two important counterexamples have been explicated by NMR structures. The KH3 domain of KSRP has a wider hydrophobic groove and key amino acid substitutions that allow it to recognize a GGGU sequence in the precursor of the 16 a b 0.1 90° 3' UTRs0.0 Introns G1 C2 A3 U4 G5 e Fcigure 1. Ribbon structure of RBFOX1 (teal) bound to UGCAUG (gold). Generated with data ● from Auwe2t.e0r et al.(2006) (PDB 2ERR). Protein–RNA hydrogen bondsG aCrAe UinGdicated in black. ● GCACG 1° UU GAAUG 1.5 UG● GUUUG GUGUG GUAUG 2° UC 1.0 ● GCUUG ● UA GCCUG ● ● 0.5 ●● ● GUAAG● GU ● ● AAAAA● ● ● ● ●● ●● GUCCG GG● 0.0 ●●● ● ● ●●● ● ●● ● ●● GCAGG ● ● ● CCCCC GC 4 43 365 1100 GA CU d 121 nM 365 nM CG 1.00 No. motifs CC 1° 1 2 6 CA 0.75 2° 1No2. 6 AU 0.50 motifs AG 17 12 0. 6 AC 25 AA 1,2 2,3 3,4 4,5 1,5 1,3 1,4 2,4 2,5 3,5 0.00 Substitution positions: GCAUG UG CG ary UG CG ary G G ary Substitution pos1iti2o3nU 4s5 CA CA on d CA CA ndo CA CA C nd Enrichment,G G ec G G ec G G ec o 0.9 1.1 1.3 1.5 S S S Figure 2. Fraction sequences bound 5mer enrichment, log2 Nucleotide pair substituted RRRReeeellllaaaattttiiiivvvveeee PPPPhhhhyyyyllllooooPPPP ssssccccoooorrrreeee Controls miRNA let-7 (Biswas et al., 2019; Nicastro et al., 2017). Later, a structure of the ZBP1 KH4 domain revealed a similar mechanism of recognition for CGGA and also found that variation in the paralogous KH domains of the ZBP/IMP family bind a core GG dinucleotide, with nucleotides in the first and fourth position of the motif determined by the variable unstructured loops around the domain (Biswas et al., 2019; Nicastro et al., 2017). GG-dinucleotide recognition could thus be a relatively common mechanism of recognition for KH domains. The CCCH ZF domain, the most common RNA-binding ZF domain, displays large variations in RNA-binding mechanisms. ZF domains contain the core amino acid sequence CX8CX5CX3H, coordinating a zinc ion in compact α-helices (Tanaka Hall, 2005). In the case of TTP, two CCCH Zinc-finger domains coordinate binding to UAUUUAUU, with each domain recognizing UAUU and utilizing substantial aromatic side-chain stacking interactions with the nucleobases to establish strong specificity; the TTP ZF domains appear to be unique compared to other discovered structures (Hudson et al., 2004). Muscleblind (MBNL) family proteins also bind RNA through CCCH ZF domains. MBNL proteins contain four ZF domains that resemble those of TTP. MBNL ZFs 1–2 and 3–4 fold with each other to adopt an unusually compact globular structure to each bind to YGCY; MBNL ZFs 1 and 3 do not even directly contribute to RNA- binding (Park et al., 2017). Evident from the atomic resolution of RBD structures in complex with RNA, RBP recognition of RNA is both diverse and complex, rendering it unpredictable by computational methods alone. Furthermore, available RBP structures tend to be in complex only with their known highest-affinity motifs, and cannot easily reveal possible recognition modes for low- and even moderate-affinity interactions. Further complicating a structural understanding of RBP binding is the abundance of intrinsically disordered regions (IDRs), including RGG boxes and RS domains within their 18 secondary structures. While RBP IDRs have traditionally been treated as linkers between modular domains that mediate the actual function of the RBP, there is increasing evidence that a significant fraction of IDRs participates in biologically relevant RNA binding (Ottoz & Berchowitz, 2020). Disordered linkers between RBDs are highly versatile, with varying length determining the effective concentration of RBDs and difficult-to-predict disorder-to-order transitions of RBP IDRs contributing to specificity and cooperativity. IDRs and repetitive regions, like globular RBDs, can also be subject to post-translational modifications that induce folding and RNA-binding. In addition, alternative splicing often acts upon IDRs to alter protein– protein interactions. Overall, IDRs contribute substantial variation to RBP–RNA interactions that must still be mechanistically clarified (Ottoz & Berchowitz, 2020; Ule & Blencowe, 2019). Compounding the variety of RBDs and their many modes of recognition is the modular structure of RBPs. Most RBPs have more than one RBD, and this protein architecture allows even greater combinatorial and emergent binding properties for a given protein (Gerstberger et al., 2014; Lunde et al., 2007). Multiple RBDs within a single protein can bind a longer linear motif (Lu et al., 2009), bind bipartite motifs of variable spacing (often related to linker length) (Nicastro et al., 2017), induce multimerization (Scheiba et al., 2014), bind RNA in different domain conformations (Wang & Tanaka Hall, 2001), and/or facilitate protein-protein interactions (Park et al., 2000; Rengifo-Gonzalez et al., 2021). Single RBDs tend to bind degenerate RNA motifs (where more than one nucleobase can be recognized at a single position within the motif) relatively weakly, but multidomain RBPs can greatly increase RBP affinity and specificity for RNA, leading to unique binding profiles (Stitzinger et al., 2021). Clearly, even highly studied RBDs and RBPs show a diversity of binding mechanisms that cannot be predicted by aa sequence alone. As such, establishing the binding preferences of any given RBP requires direct, 19 empirical measurement. Such measurements not only inform RBP target prediction, thus elucidating gene regulatory networks, but may also contribute to the understanding of the correspondence between aa sequence and RNA recognition. Quantitative, high-throughput biochemistry to describe specificity and affinity The vast diversity of RBPs and the RNA sequences they bind makes the determination of genuine regulatory targets a major challenge in the field of post-transcriptional gene regulation. In vivo techniques to detect global binding sites of a particular RBP often rely on CLIP-seq. Meta-analysis of dozens of splicing factors and their targets shows that RBPs bind functionally coherent sets of transcripts to execute regulatory programs in cells (Hogan et al., 2008; Mukherjee et al., 2019); an RBP with its putative targets can thus be described as a biological network. At its simplest, a network can be comprised of gene targets as described by a single high-throughput method such as CLIP-seq, but can be refined by integrating multiple data types such as affinity/specificity estimates, evolutionary conservation, gene ontology, and experimental measurements. RBP regulatory networks have been described for Nova and Rbfox proteins, among others (Weyn-Vanhentenryck et al., 2014; Zhang et al., 2010), and can be extended to gene programs, responses to stimuli, and tissues (Ule & Blencowe, 2019). Network analysis of RBP–RNA interactions suggests that RBPs form two main architectures: (1) clusters of RBPs that bind the same transcripts, imparting robust transcript regulation through cooperative and competitive interactions with other RBPs in that cluster, and (2) hierarchical chains of RBPs connecting these clusters that convey directional regulatory information in a less redundant manner (Quattrone & Dassi, 2019). Defining these complex networks is a particular challenge of the field, given their redundancy and dynamism: post-transcriptional gene 20 regulatory networks can vary tremendously given developmental stage, cell type, and cellular environment. Although CLIP identifies bona fide target transcripts for a given RBP, CLIP-based assays can also have high levels of non-specific background and are subject to biases arising differential crosslinking efficiency, differential susceptibility to nuclease digestion, and library amplification (Kishore et al., 2011; Sundararaman et al., 2016; Uren et al., 2012; Van Nostrand et al., 2016). In addition, any single CLIP experiment is obligately cell type- and stage-specific, limiting its generalizability. To accurately determine RBP targets requires precise, quantitative methods to most readily describe binding networks across the spectrum of cell types and states and to accurately describe gene regulatory networks (Ye & Jankowsky, 2020). While many factors influence the binding of RBPs to their target RNAs, the most informative measurements for understanding the cellular targets of an RBP are those providing the biochemical affinities of that RBP for a diversity of RNA sequences. Indeed, many of the aspects by which cellular contexts differ from one another, like differential RBP or cofactor expression, differential substrate levels, or subcellular localization, can be quantitatively predicted with known affinity measurements (Jens & Rajewsky, 2015). Although recent progress has been made in the measurement of RBP binding kinetics in vivo (Sharma et al., 2021), systematic in vitro studies of RBP specificity provide the clearest building blocks on which to construct these complex networks: the intrinsic, biochemical properties of RBP specificity. High-throughput in vitro methods characterizing RBP–RNA interactions typically consist of a large library of RNA molecules incubated with a single RBP of interest, coupled with a mechanism to isolate and sequence bound RNA substrates (Ye & Jankowsky, 2020). Systematic efforts to catalogue RBP motifs have not only recapitulated known RBP binding preferences, but 21 have invited a more nuanced understanding of RBP binding, revealing non-canonical linear motifs as well as RBP-specific preferences for secondary structure, interspersed motifs, binding in multiple registers, and biased flanking nucleotide composition as important determinants of RBP specificity (Dominguez, et al., 2018; Lambert et al., 2014; Ray et al., 2013). Intrinsic specificity must also be understood not as a binary binding outcome (‘bound’ or ‘unbound’), but as a continuum of affinity for a variety of RNA sequences, which can be provided by quantitative analyses of high-throughput in vitro RBP binding assays. One such assay, RNA Bind-n-Seq (RBNS), provides simultaneous, quantitative measurement of binding for an RBP incubated with billions of RNA sequences in vitro (Dominguez, et al., 2018; Lambert et al., 2014; McGeary et al., 2019). RBNS produces precise measurements of binding for RBPs with a variety of domain structures and specificities. Additionally, a variation of RBNS utilizing a synthetic library design based on naturally occurring transcript sequences (nsRBNS) has also been performed, providing insight into the relationship between binding affinity and structural accessibility for endogenous intronic sequences (Taliaferro et al., 2016). While this restricts the sequence space available to be queried in the assay, the sequence composition and oligonucleotide length aid in the discovery of motifs that appear frequently in the transcriptome in arrangements compatible with binding, and they should tend to be those with the largest transcriptome-wide regulatory effects. Recent quantitative studies have meaningfully advanced the study of RBP binding preferences and their implications for biology. For instance, a high-throughput sequencing kinetics approach (HITS-KIN) incubating an apparently non-specific Escherichia coli tRNA- binding protein with ~150,000 RNA sequences found that C5 protein could in fact discriminate among RNA sequences—but its physiological targets (all tRNA precursors) consist of moderate- 22 , and not high-, affinity motifs (Guenther et al., 2013). RBNS of 126 RBPs showed that RBPs that bind similar top motifs in vitro bind significantly different sets of transcripts harboring that motif in vivo, revealing substantial unaccounted-for binding variation inexplicable by their highest-affinity primary sequence preferences (Dominguez, et al., 2018). Furthermore, RBNS of six AGO–miRNA complexes analyzed by a biochemical model of miRNA-mediated repression established that each miRNA had non-canonical sites that were specific to that AGO–miRNA complex; furthermore, these non-canonical motifs functionally contributed to cellular repression of their target mRNAs (McGeary et al., 2019). Novel and consequential binding behaviors have been revealed by detailed biochemical studies even for well-studied RBPs. A study of human Pumilio PUM1 and PUM2 using RNA- MaP, which utilizes a modified next-generation sequencing platform to produce thermodynamic and kinetic binding constants, uncovered previously-uncharacterized binding determinants for the protein, revealing multi-register binding, more extensive base-flipping, and coupling among nucleotide preferences within the linear motif; all of these biochemical considerations impact cellular motif occupancies (Jarmoskaite et al., 2019). Their thermodynamic binding model correlated well with in vivo binding sites measured using enhanced CLIP (eCLIP) signal enrichment, underscoring the relevance of these biochemical properties and further predicting that about two thirds of cellular PUM2 binds noncanonical motifs at physiological (sub- saturating) concentration, making PUM2 binding particularly sensitive to changes in both RBP expression and RNA binding affinity. When the concentration of an RBP is very high, protein occupancy will increase at moderate- and low-affinity binding sites (Jankowsky & Harris, 2015; Jarmoskaite et al., 2019; Lambert et al., 2014). Thus, RBP specificity is well-described as an 23 affinity landscape sensitive to cellular RBP concentrations, enabling more finely tuned gene regulation. Suboptimal motifs in gene regulation In thermodynamic equilibrium, the concentration of any particular RNA–RBP complex will be directly proportional to the concentration of both the unbound RBP and the unbound RNA sequence. However, the concentration of the unbound RBP will itself be a function of its own affinity landscape and the corresponding cellular pool of available RNA, and the concentration of the unbound local RNA sequence will likewise be a function of its own affinity for the expressed RBP proteome. For a highly-expressed RBP with a specific, high-affinity motif, the concentration of that RNA sequence within the cellular pool may be much lower than the concentration of the RBP, such that the RBP will also substantially bind moderate- and low- affinity RNA sequences. By contrast, RBPs whose highest-affinity binding is characterized by degenerate or shorter motifs with more moderate absolute affinity will be more evenly distributed across the RNA pool, such that the concentration of any individual RNA species would be less influential on the binding distribution (Jens & Rajewsky, 2015). The concentration of the other reactant, the RBP, will also drive the reaction rate; large changes in RBP expression are common in dynamic biological processes like development (Weyn-Vanhentenryck et al., 2018) and can thus drive differential occupancies at intermediate-affinity, high-frequency binding sites. This behavior could also be relevant in regimes of high local protein concentration. RBPs are increasingly associated with phase-separated compartments, which rely on diverse multivalent protein and RNA interactions to form (Gomes & Shorter, 2019; Sanders et al., 2020; 24 Ying et al., 2017); suboptimal motifs may contribute to the formation, maintenance, or interactions of cellular microenvironments (Tsai et al., 2017). There is substantial evidence for the utility of suboptimal motifs for transcription factors (TFs) and other DNA-binding proteins, particularly during early development, when precise spatial and temporal gene expression is essential. Morphogens, discovered in 1988 in the Drosophila melanogaster embryo, introduced the field to the regulatory potential of gradients of molecular activity in post-transcriptional gene regulation (Driever & Nüsslein-Volhard, 1988). The physical concentration gradient established by the passive diffusion of the bicaudal mRNA establishes the anterior–posterior orientation of the embryo. A function for DNA suboptimal motifs in tandem with a morphogen gradient in D. melanogaster would also be observed: the spatiotemporal information encoded in the concentration of the morphogen is mechanistically transduced using the spectrum of target affinities of the downstream TF such that binding to high-, medium- and low-affinity sites occurs in cells closest to the morphogen source, and only the highest-affinity sites are occupied further along the morphogen axis (Jiang & Levine, 1993). A thermodynamic model of gene expression computationally applied to segmentation genes in D. melanogaster confirmed the widespread importance of low-affinity sites for transcriptional control of metazoan body patterning (Segal et al., 2008). Further research has extended this paradigm to other TFs, finding that clusters of low-affinity sites are required for proper Ultrabithorax, Pax2, and Senseless function at enhancers (Crocker et al., 2015; Zandvakili et al., 2018). Concurrent with the characterization of bicaudal, the promoters of a histone gene family in sea urchin were found to contain low- or high-affinity Sp1-like TF motifs differentiating their gene expression between the early- or late-blastoderm stage of development (Lai et al., 1988). 25 Later, it was shown that the Caenhorhabditis elegans TF PHA-4 activates its primary motif targets in early pharyngeal development, while its secondary targets are activated later, when protein levels are higher (Gaudet & Mango, 2002). In mammals, the Mus musculus TF Pax6 operates through both a concentration gradient in the nervous system and regulates Neurogenin2 through a low-affinity enhancer element when Pax6 protein concentration is high (Scardigli et al., 2003). Pax6 expression itself is regulated by a low-affinity enhancer element bound by the TF PREP1, which activates its enhancers sequentially according to their affinity during eye lens development (Rowan et al., 2010). At the genomic level, chromatin immunoprecipitation and high-throughput sequencing has revealed extensive low-affinity TF binding sites in the genome, and evolutionary analysis spanning yeast to human genomes suggests widespread functionality for suboptimal motifs (Jaeger et al., 2010; Tanay, 2006). A high-throughput, biochemical analysis of the yeast starvation-response master TF, GCN4, identified distinct functions for target genes activity by high- versus low-affinity motifs: genes with low-affinity motifs were associated with longer-term starvation response genes, whereas genes with high-affinity motifs were associated with a rapid response to amino acid starvation (Nutiu et al., 2011). Most recently, a Bayesian biophysical model that considered TF binding to both high- and low-affinity sites better predicted gene expression in human B-cells (Wang et al., 2015). Although possible suboptimal motifs have been observed in transcriptome-wide binding studies of RBPs, evidence of RBP suboptimal motifs with clear biological functions has been elusive. The C. elegans GLD-1 RBP represses translation in accordance with the affinity of its binding sites in 3′ UTRs (Wright et al., 2011), and work on hnRNP H and MBNL is suggestive of dose-dependent splicing through cis-elements (Wagner et al., 2016; Xiao et al., 2009). Suboptimal motifs for RBPs may have a variety of functions. They could be additive and, 26 relatedly, a product of the incremental evolution of RBP binding sites. Clusters of suboptimal motifs may enhance RBP specificity for a transcript, enhance binding to nearby primary motifs, or provide biological buffering in post-transcriptional gene regulation. As seen for C. elegans GLD-1 and a variety of TFs, suboptimal motifs can tune levels of activity for subsets of transcripts. Lastly, suboptimal cis-regulatory elements may function only at high levels of a trans factor to narrow the temporal or spatial scope of RBP activity. As high-throughput biochemical studies reveal that many RBPs meaningfully bind suboptimal motifs (Jarmoskaite et al., 2019; McGeary et al., 2019), it will be important to integrate complete motif affinities alongside RBP and motif concentrations–extending consideration to abundant sites of lower affinity when RBP activity is high–to fully understand context-specific splicing. 27 28 Chapter 2. Concentration-dependent splicing is enabled by Rbfox motifs of intermediate affinity Bridget E. Begg1, Marvin Jens1, Peter Y. Wang1,2, Christine M. Minor1, Christopher B. Burge1 1Department of Biology, Massachusetts Institute of Technology, Cambridge MA 02139, USA 2Whitehead Institute for Biomedical Research, Cambridge, MA 02142, USA B.E.B. performed experiments, analyzed data, and wrote the manuscript. M.J. analyzed neuronal subtype data, developed the biochemical model, and wrote the manuscript. P.Y.W. performed flow cytometry experiments. C.M.M. assisted with filter binding experiments. C.B.B. supervised the study and wrote the manuscript. Published as: Begg, B.E., Jens, M., Wang, P.Y., Minor, C.M., and Burge, C.B. (2020). Concentration- dependent splicing is enabled by Rbfox motifs of intermediate affinity. Nature Structural and Molecular Biology 27(10): 901-912. 29 Abstract The Rbfox family of splicing factors regulate alternative splicing during animal development and in disease, impacting thousands of exons in the maturing brain, heart, and muscle. Rbfox proteins have long been known to bind to the RNA sequence GCAUG with high affinity, but just half of Rbfox binding sites contain a GCAUG motif in vivo. We incubated recombinant RBFOX2 with over 60,000 mouse and human transcriptomic sequences to reveal substantial binding to several moderate-affinity, non-GCAYG sites at a physiologically relevant range of RBFOX concentrations. We find that many of these “secondary motifs” bind Rbfox robustly in cells and that several together can exert regulation comparable to GCAUG in a trichromatic splicing reporter assay. Furthermore, secondary motifs regulate RNA splicing in neuronal development and in neuronal subtypes where cellular Rbfox concentrations are highest, enabling a second wave of splicing changes as Rbfox levels increase. Introduction The Rbfox family of RNA-binding proteins (RBPs) have been the subjects of numerous genetic, biochemical, and structural studies since their discovery over twenty-five years ago. Rbfox1 was identified in 1994 as the genomic locus feminizing on X (fox-1) in C. elegans; it was subsequently shown that the encoded protein peaks in expression in larval development and controls sex determination by repressing the dosage compensation factor xol-1 post- transcriptionally (Hodgkin et al., 1994; Skipper et al., 1999). Mammalian Rbfox1 (A2BP1) and its paralogs Rbfox2 (RBM9) and Rbfox3 (NeuN) are highly expressed in heart, skeletal muscle, and in the brain, with a characteristic spike in expression in late neuronal development (Conboy, 2017; Gallagher et al., 2011; Kim et al., 2009; Weyn-Vanhentenryck et al., 2018). While Rbfox 30 proteins are predominantly nuclear and regulate pre-mRNA splicing, some isoforms are also expressed in the cytoplasm, where they regulate RNA stability (Kuroyanagi, 2009). The three mammalian paralogs have high sequence identity (94–99%) in their RNA-binding domains (RBDs), substantially overlapping gene targets, and similar activities, making single knockouts difficult to interpret and regulatory targets challenging to define (Gehman et al., 2011; Jacko et al., 2018). Mouse Rbfox1 central nervous system (CNS) knockouts have seizures and defects in neuronal excitability, while Rbfox2 CNS knockout mice exhibit defects in cerebellar development (Gehman et al., 2011, 2012). Binding targets of RBFOX1, both mRNA splicing- and mRNA stability-related, are enriched for genes related to autism, axon and dendrite formation, and electrophysiology in mice (Hamada et al., 2016; Lee et al., 2016; Vuong et al., 2018; Weyn-Vanhentenryck et al., 2014), and mutations in Rbfox1 and Rbfox3 have been associated with autism and epilepsy in humans (Barnby et al., 2005; Bhalla et al., 2004; Martin et al., 2007; Sebat et al., 2007). Recently, it was shown that Rbfox triple knockout mouse embryonic stem cells fail to develop a mature splicing profile when differentiated into ventral motor neurons, underscoring the roles of these proteins in neuronal maturation (Jacko et al., 2018). In 2003, RBFOX1 was found to bind the pentanucleotide GCAUG with high affinity and regulate alterative exon inclusion through binding to flanking intronic regions (Jin et al., 2003). (U)GCAUG had long been established as a prominent signal in computational analyses of neuronal alternative splicing, and Rbfox proteins have been shown to regulate numerous alternative exons through this motif in a variety of neuronal subtypes (Brudno, 2001; Jacko et al., 2018; Minovitsky et al., 2005; Weyn-Vanhentenryck et al., 2014). Binding to the GCAUG motif is mediated by canonical and non-canonical interactions between the protein’s single RNA 31 recognition motif (RRM) and target RNA (Ying et al., 2017). Generally, binding of Rbfox proteins to motifs in the proximal ~200 nucleotides (nt) of the intron downstream of an alternative exon activates exon inclusion, while binding in the proximal upstream intron or alternative exon promotes exon exclusion (Yeo et al., 2009). However, more distal motifs are also able to regulate splicing23,24. Proteins recruited by the Ala/Tyr/Gly-rich C-terminal domain of Rbfox mediate its effects on splicing, whereas its effects on stability may be related to competition with other RBPs and microRNAs (Sun et al., 2012; Vuong et al., 2018). Although Rbfox binding to canonical GCAUG motifs and rarer GCACG motifs is well characterized, studies of binding sites in vivo using crosslinking-immunoprecipitation (CLIP) have observed that about half of CLIP peaks lack an associated GCAYG (Y = C or U) motif (Yeo et al., 2009), suggesting the existence of additional binding determinants (Dominguez, et al., 2018). Indeed, recent studies have noted GUGUG motifs and motifs differing by one base from the canonical GCAUG in CLIP peaks (Jangi et al., 2014; Lambert et al., 2014; Lee et al., 2016; Stoltz, 2015), and RBFOX1 was shown to compete with MBNL1 at a GCCUG motif (Sellier et al., 2018). It has been proposed that binding to GUGUG is mediated by partner proteins such as SUP12 or members of the large assembly of splicing regulators (LASR) complex (Damianov et al., 2016; Kuroyanagi et al., 2007; Kuwasako et al., 2014). Direct binding of Rbfox proteins to these motifs has not been excluded, but if such binding occurs it is not known to have any functional consequence. Almost all studies of Rbfox binding to date have filtered their datasets for the presence of a GCAUG element, discarding many other CLIP peaks (Gehman et al., 2011, 2012; Jacko et al., 2018; Lee et al., 2016; Vuong et al., 2018; Weyn- Vanhentenryck et al., 2014; Zhang et al., 2008). 32 Studies of Rbfox in vivo binding and regulation using CLIP have identified important regulatory targets and uncovered complex regulatory networks. However, CLIP is known to have substantial rates of false positives and false negatives (Kishore et al., 2011; Uren et al., 2012; Van Nostrand et al., 2016), making it challenging to confidently infer binding sites. Here we employ RNA Bind-n-Seq with natural sequences (nsRBNS) as a biochemical approach to refine our understanding of Rbfox binding across mammalian transcriptomes (Lambert et al., 2014; Taliaferro et al., 2016). Our method detects strong binding to the canonical Rbfox motif GCAYG, but also binding to several additional motif variants. We show that Rbfox proteins bind to these “secondary motifs” in vivo, and exert regulatory activity in a manner dependent on Rbfox concentration. Furthermore, we find evidence of important roles for secondary motifs in splicing programs involved in neuronal differentiation and subtype specification, and we propose a model for Rbfox regulation that incorporates these motifs. Natural sequence RBNS recovers known features of RBFOX2 binding To better understand the transcriptomic RNA binding preferences of RBFOX2 across the transcriptome, we designed an nsRBNS library based on naturally-occurring mammalian 3' untranslated region (UTR) sequences (Figure 1a). To construct the library, ~2,200 well- annotated human and mouse 3' UTRs were selected that matched the base and 5mer composition of the entire 3' UTR transcriptome (Figure 1b). The 3' UTRs of each of these transcripts were divided into overlapping 110-base segments (Supplementary Table 1), which are longer than those used in random sequence RBNS (typically 20–40 nt), enabling analysis of binding of RBFOX2 to motif clusters. We performed nsRBNS with recombinant RBFOX2 at 4 nM, 14 nM, 43 nM, 121 nM (in 33 34 Figure 1. 3' UTR natural sequence nsRBNS (nsRBNS) with RBFOX2 captures variation in binding affinity. a. Schematic of nsRBNS. Recombinant protein is incubated with a designed RNA library to equilibrium and bound RBP:RNA complexes are purified. Oligonucleotides are sequenced and the enrichment (R) value is calculated ((reads per million)input/(reads per million)pulldown). b. nsRBNS 5mer frequencies correlated with 5mer frequencies of the 3' UTR transcriptome (n = 1024 5mers). Pearson correlation. c. Distribution of R for nsRBNS sequences containing 0 (n = 49931), 1 (n = 5586), 2 (n = 392), or 3+ (n = 22) Rbfox GCAUG motifs and 0 (n = 54637) or 1-3 (n = 1294) GCACG motifs *** P < 0.001 between lowest and highest counts (two-sided Wilcoxon Rank-Sum test). d. Distribution of enrichment of RBFOX2 eCLIP reads in HepG2 cells with increased motif count for 0 (n = 8004), 1 (n = 1244), or 2+ (n = 118) GCAUG motifs and 0 (n = 9065) or 1-3 (n = 301) GCACG motifs in transcriptomic regions corresponding to those in nsRBNS library (normalized to IgG control). *** P < 0.001, * P < 0.05 between lowest and highest counts (two-sided Wilcoxon Rank-Sum test). e. An iterative method discovers moderate binding by RBFOX2 to six motifs of the sequence format GHNUG (teal) beyond two known Rbfox motifs (gold). After nine rounds of enrichment analysis, remaining GNNUG 5mers (teal) were also included as secondary motifs, while AU-rich (grey) and shifted (light blue) 5mers were excluded from subsequent analyses. f. Distribution of R for nsRBNS sequences containing 0 (n = 25501), 1 (n = 18682), 2 (n = 7957), 3 (n = 2576), 4-6 (n = 1069), or 7-14 (n = 146) secondary motifs. *** P < 0.001 between lowest and highest counts (two-sided Wilcoxon Rank-Sum test). g. Distribution of enrichment of RBFOX2 eCLIP reads in HepG2 cells at 0 (n = 3922), 1 (n = 3188), 2 (n = 1453), 3 (n = 498), 4-6 (n = 262), or 7-14 (n = 43) secondary motifs at library positions in the transcriptome (normalized to IgG control). *** P < 0.001 between lowest and highest counts (two-sided Wilcoxon Rank-Sum). 35 technical duplicate), 365 nM and 1.1µM, as well as a 0 nM (no protein) control, with a constant 250 nM RNA concentration. This range, spanning 2.5 orders of magnitude, is comparable to the natural range of total Rbfox family expression, which varies at the mRNA level by at least three orders of magnitude across cell types4, and by at least an order of magnitude during neuronal differentiation (McNutt et al., 2013). Enrichment (R) values were calculated for each individual oligonucleotide as the ratio of the frequency in the bound pool of RNA over the frequency in control (no protein) conditions (Supplementary Table 2). R values of bound oligonucleotides at similar concentrations were strongly correlated (r ≥ 0.74) (Extended Data Figure 1a), indicating high reproducibility of oligonucleotide-specific binding. Of the 55,931 sequences (87% of the total) that were detectable at 365 nM, 8,412 (15%) bound with R ≥ 1.1 (10% enrichment) about half of which (4,032) contained a canonical GCAYG motif. The distribution of R values was shifted upwards with increasing counts of each motif (Figure 1c, P < 2.01 x 10-13 and P < 5.64 x 10-46, respectively), consistent with previous work38. Similar trends were observed when considering in vivo binding, assessed by the density of reads mapping to the genomic region of each oligonucleotide, using published eCLIP data (Figure 1d, P < 2.7 x 10-23 and P < 0.10, respectively) (Van Nostrand et al., 2016). The nsRBNS R values of sequences in the library that contained a single GCAUG had a moderate positive correlation with eCLIP read density in the corresponding mRNA regions (r = 0.28, P < 1.1 x 10-25, Extended Data Figure 1b), as observed previously (Taliaferro et al., 2016), confirming that this assay captures some aspects of motif context that are relevant in vivo. RBFOX2 binds a set of secondary motifs of intermediate affinity in vitro 36 We aimed to systematically identify sequence elements beyond the “primary” GCAYG motifs that might help to explain the large dynamic range of RBFOX2 binding and the presence of many oligonucleotides with binding above background that lacked GCAYG motifs. We employed an iterative analysis, identifying the most enriched 5mer (GCAUG), then removing all oligonucleotides containing this 5mer and recalculating motif enrichments from the remaining oligonucleotides (Figure 1e). This procedure, similar to that employed in the SKA algorithm for random-library RBNS (Lambert et al., 2014), is designed to prevent the spurious detection of 5mers that are enriched due to overlap with a higher-affinity 5mer. We performed the iterative enrichment analysis on data from the highest RBFOX2 concentration, 1.1 µM, which should favor binding to lower-affinity sites (Lambert et al., 2014). At each of the first four iterations, the most enriched 5mers matched the pattern GNNYG. At iterations 5–9, 5mers matched either: GHWUG (GUAUG, GCUUG), where H = A, C or U; NGCAU (UGCAU), where N = A, C, G or U; or W5 (AAAAA, UUUUU), where W = A or U. We considered these three different classes of potential binding motifs, which might bind RBFOX2 directly or be enriched for other reasons. To cast a wide net, we also considered other 5mers matching these patterns that had R values at least 2.5 standard deviations above the mean in the remaining oligonucleotides, adding the GNNUG 5mers GCCUG and GUGUG, the NGCAU 5mers CGCAU and AGCAU, and the W5 motifs UAUAU and AUAUA to our list of potential binding motifs. At the oligonucleotide level, the six GHNUG secondary motifs robustly increased the R-value of an oligonucleotide when several motif occurrences were present (Figure 1f, P < 9.3 x 10-25). The NGCAU and W5 5mers did so to a somewhat lesser extent (Extended Data Figure 1c-d; P < 8.8 x 10-28 and P < 0.001, respectively). We also observed increased eCLIP enrichment in 37 regions with increasing numbers of GHNUG motifs (Figure 1g, P < 6.1 x 10-7) (Feingold et al., 2004). However, NGCAU and W5 5mers are not robustly enriched in eCLIP data (Extended Data Figure 1e-f, P < 0.071, P < 0.017, respectively). Because of the highly significant binding signal for the six GHNUG 5mers in both nsRBNS and eCLIP, we chose to pursue only this subclass of secondary motifs for further study. Together, these six non-GCAYG 5mers comprised our set of candidate “secondary” Rbfox motifs. Biochemical characterization of Rbfox secondary motifs Our iterative method suggested that the G1, G5 and U4 positions of the Rbfox motif are the most critical sites of recognition for the protein. This observation is consistent with published structures of RBFOX1 bound to (U)GCAUG(U), in which the primary hydrogen bonding contacts occur at G1, U4, and G5, while the other bases are recognized mostly by shape (Figure 2a) (Auweter et al., 2006). Rbfox tolerating nucleotide variants at positions 2 and 3 is also consistent with the conservation pattern of GCAUG in introns and 3' UTRs, where we observed that the G1 and G5 positions are most strongly conserved, with the intervening positions showing less constraint (Figure 2b). Examining the nsRBNS enrichment of 6mers ending with the six GHNUG 5mers, we observed a preference for C or U at the first position, consistent with the preference of GCAYG motifs for an upstream U or C (Extended Data Figure 2a) (Auweter et al., 2006). To further explore these candidate motifs, we reanalyzed data from previous random-library RBNS and intronic nsRBNS experiments performed with RBFOX2 (Lambert et al., 2014; Taliaferro et al., 2016). In random-library RBNS, which employs much shorter sequences that are unlikely to contain more than a single motif by chance, secondary motifs were only weakly 38 a b 0.1 90° 3' UTRs0.0 Introns G1 C2 A3 U4 G5 c e ● 2.0 GCAUG GCACG 1° UU● GAAUG 1.5 UG● GUUUG GUGUG 2° UC 1.0 ● GUAUGGCUUG ● UA GCCUG ● ● 0.5 ●● ● GUAAG● GU ● ● AAAAA● ● ● ● ●● ●● GUCCG GG● 0.0 ●●● ● ● ●●●● ● ●● ●● GCAGG ● ● ● CCCCC GC 4 43 365 1100 GA CU d 121 nM 365 nM CG 1.00 No. motifs CC 1° 1 2 6 CA 0.75 2° 1No2. 6 AU 0.50 motifs AG1 2 6 AC0.25 AA 1,2 2,3 3,4 4,5 1,5 1,3 1,4 2,4 2,5 3,5 0.00 y y y SubsStitu GCAUG UG CG ar UG G ar G G r ub tisotnit upotisointio pnos:s1iti2o3n4s5 CA CA on d CA A C U a G c C on d C G G G c GC A CA d e e G co n Enrichment, e 0.9 1.1 1.3 1.5S S S Figure 2. 39 Fraction sequences bound 5mer enrichment, log2 Nucleotide pair substituted RRRReeeellllaaaattttiiiivvvveeee PPPPhhhhyyyyllllooooPPPP ssssccccoooorrrreeee Controls Figure 2. Rbfox proteins reproducibly bind a class of secondary motifs with moderate affinity. a. Ribbon structure of RBFOX1 (teal) bound to UGCAUG (gold). Generated with data from Auweter et al. (2006) (PDB 2ERR). Protein–RNA hydrogen bonds are indicated in black. b. Relative per-base conservation, as represented by PhyloP score, for each position of GCAUG at all instances of the motif in 3' UTRs (dark teal, n = 26707) and introns (light teal, n = 91791). c. Primary (gold), secondary (teal), polyA and polyC (light grey), and GCAGG, GUAAG, and GUCCG (dark grey) R values are shown across four concentrations of RBFOX2 nsRBNS experiments. d. Analysis of the fraction of oligonucleotides bound in nsRBNS at three concentrations of RBFOX2 for primary (gold (GCAUG and GCACG), nGCAUG-1 = 5435, nGCAUG- 2 = 384, nGCACG-1 = 1246, nGCACG-2 = 30) and secondary (teal, nSecondary-1 = 15676, nSecondary-2 = 6722, nSecondary-6 = 64) motifs. An oligonucleotide was considered bound if it had an R value of at least 1.1. e. nsRBNS R values at 1.1 µM RBFOX2 concentration for all pentamers diverging from GCAUG or GCACG by 1 or 2 bases. Secondary motifs identified here are outlined in black. 40 enriched (Extended Data Figure 2b). However, in an nsRBNS experiment using highly conserved intronic regions, U-rich secondary motifs GCUUG, GUUUG, and GUGUG were enriched (Extended Data Figure 2c). Thus, length and sequence composition of the library used in nsRBNS may impact detection of specific secondary motifs. Our design with varying concentrations of RBFOX2 enabled analysis of relative RBP occupancies and saturation. As seen previously (Lambert et al., 2014), Rbfox primary motifs had high enrichment that declined at the highest protein concentration (1.1 µM), consistent with these sites reaching saturation around the 2nd-highest concentration, 365 nM (Figure 2c). Control 5mers (e.g., AAAAA, CCCCC) or 5mers differing from GCAUG by one substitution but not identified as secondary motifs (GCAGG, GUAAG, GUCCG) had low R values across concentrations, suggesting the absence of specific binding. By contrast, the secondary motifs identified above exhibited steadily increasing enrichment at successive protein concentrations, consistent with moderate-affinity binding that becomes saturated at or above 1.1 µM [RBFOX2]. We noted that about half of the oligonucleotides containing six or more secondary motifs were detectably bound at 1.1 µM [RBFOX2], comparable to the fraction of oligonucleotides containing a single GCAUG that were bound at lower Rbfox concentrations (Figure 2d). This observation suggests that clusters of secondary motifs may function in regulation similarly to single primary motifs when Rbfox concentrations are high. Notably, in a previous study using surface plasmon resonance, several secondary motifs had Kd values in the 100–600 nM range (Stoltz, 2015b), well above the low nanomolar value observed for GCAUG but comparable to the highest-affinity motifs of many other RNA-binding proteins (Helder et al., 2016). Of the 4,358 sequences bound at R ≥ 1.1 in nsRBNS that lacked a primary motif, 1,885 (43%) contain at least one of the GHNUG secondary motifs. Thus, the six secondary motifs identified appear to 41 explain a plurality of the non-canonical (non-GCAYG) RBFOX2 binding to our nsRBNS library. To further explore whether other related 5mers beyond those identified above are specifically bound by RBFOX2, we examined the binding at 1.1 µM [RBFOX2] of all 5mers with one or two positions substituted relative to the primary motif GCAUG (Figure 2e). The most enriched 5mers consisted predominantly of the two primary and six secondary motifs. Besides these, certain AU- rich 5mers with two differences from GCAUG (UUAUG, AUAUG, GUAUU, and GUAUA) showed modest enrichment. We tested RBFOX2 binding independently using a nitrocellulose and nylon filter binding assay (Extended Data Figure 2d, Supplementary Table 3). This assay showed stronger binding to radio-labeled primary motifs than to any secondary motifs, as expected. The secondary motifs GCUUG, GAAUG, GUUUG and GUAUG bound more strongly than GUGUG or any of the negative controls (GCAGG, GCUAG, GUCCG, U23), most of which exhibited minimal or no binding. RBFOX2 binds to specific secondary motifs in vivo We next asked whether Rbfox binding to individual secondary motifs could also be observed in vivo. We analyzed published binding data generated using individual-nucleotide resolution CLIP (iCLIP) from mouse embryonic stem cells (mESCs) (Jangi et al., 2014). We constructed meta-motif plots, summing read density as a function of distance for every instance of a given 5mer in expressed introns or 3' UTRs (Figure 3a). Primary motifs and all six GHNUG secondary motifs showed a characteristic peak located somewhat 3' of the motif location, as expected (Lambert et al., 2014), with broader peaks observed in a few cases (e.g., GAAUG). Some of the motifs also showed evidence of crosslinking at U4 of the motif (Extended Data Figure 3a). The 42 a b GCAUG GCACG ●1.5 ● 3.5 3.5 3.0 3.0 ● 2.5 2.5 1.0 2.0 2.0 ● ● ● ● ● ● 1.5 1.5 ● ● ● ●● ● ● ● ●●●● ● 1.0 1.0 ●●● ●●●●●●●●●● ● ●● 0 25 50 0 25 50 0.5 ●●●●●●● ●● ● ●●●● ●● ● ● ●●●● ● ●●●●●●●●●● ●●●● ● ●●● ●●●●●● ● ●●●● ●●●●●● ●● ● ●●●●●●●●●●●● ●● ●●● GCUUG GAAUG ●●● ● ●●●● ●● ● ●● ●●●●●●●● ●●●●●●●● ● ●●●●●●●●●● ● Primary●●●●●●●●●●●●● ● 1.5 1.5 0.0 ●● ●●●●● ●●●●●●●●●● ● ● ● ●●●●●●●●●● ●●●● Secondary ● ●●●●●●●●●●● ●●●●●●●●●● ●● ●● GCAUG-overlapping 1.4 1.4 ●●●●●●● ●●● ● ●● ●●●●●●●●● ●●●●●●●● ● ●●●●● Other●●●●●●● ● ●●●●●●●●●●●●● ●●●●●● 1.3 1.3 ●●●●● ● ● ●● ●● ●●● ●●●●●●● ●●●●●●●● ● ● ● ●● ●● ●● ● 1.2 1.2 ●●●●●● ●●●●●●● ●●●●●●●●●● ●●●● ●●●●●●●●●●●●●●● 1.1 1.1 ● ●●●●●●●●●● ●●●● ●●●0● 1 2 ●●●● 1.0 1.0 ● ●●●●●● 0 25 50 0 25 50 c GUUUG GUAUG 3'UTR Intron 1.5 1.5 1.4 1.4 1.3 1.3 1.2 1.2 1.1 1.1 1.0 1.0 0 25 50 0 25 50 GAAUG 2 GUGUG GCCUG GUUUG 1 1.5 1.5 0 1.4 1.4 GUAUG 1.3 1.3 –2 1.2 1.2 GUGUG 1.1 1.1 1.0 1.0 0 25 50 0 25 50 d Intron e WWhole 6 br eayinn 100% 4 75%Primary Secondary 4.4X GCAUG-overlapping 5.4X Other 2.4X 50% 2 12.1X 2.7X 8.0X 3.7X 25% 5.3X 0 0 50 100 150 200 5.4X 7.7X 8.1X 19.1X 5mer rank 0% 1.2x 3' UTR 3' UTR Figure 3. 43 5mer enrichment Figure 3. Rbfox proteins bind secondary motifs in vivo. a. Primary and secondary motif reads peak near 0 in a metaplot centered at the motif in introns (black) and 3' UTRs (grey) in RBFOX2 iCLIP data (Jangi et al., 2014) (GEO GSE54794). iCLIP reads containing the motif of interest (see Methods for read counts) were aligned with position one of the pentamer at 0 and normalized to the minimum read count in an 80-nt window (50-nt window shown). Y-axis range was reduced for secondary motifs. b. Correlation of intronic iCLIP- and nsRBNS-enriched 5mers (n = 1024). Secondary motifs indicated in teal, primary motifs indicated in gold. Grey dots indicate “hitchhiking” motifs that overlap the primary motif GCAUG by at least three bases but do not have intrinsic Rbfox affinity. c. RBFOX1 HiTS-CLIP data(Jacko et al., 2018; Weyn- Vanhentenryck et al., 2014) (SRA SRP128054, SRP035321) enrichments for primary motifs and secondary six motifs in both 3' UTRs (left, n = 2963 and 989) and introns (right, n = 847 and 1431) relative to transcriptomic frequencies in mouse whole brain and ventral spinal neurons, respectively. Enrichment was calculated based on the 5mer composition of 100-base CLIP peak regions centered around the apex of the CLIP peak relative to 5mer composition of the transcriptomic region. d. Four secondary motifs (GCUUG, GAAUG, GUUUG, GUAUG) are indicated among the top 200 highly enriched 5mers derived from intronic HiTS-CLIP peaks from mouse ventral spinal neurons. Primary motifs in gold, secondary motifs in teal. Peaks calculated as above. e. High-confidence CLIP peaks in two different Rbfox1 HiTS-CLIP datasets (SRA SRP128054, SRP035321) in ventral spinal neurons(Jacko et al., 2018b) and mouse whole brain cells (Weyn-Vanhentenryck et al., 2014) attributable to primary (light gold), or four secondary (teal; GCUUG, GAAUG, GUUUG, and GUAUG), or both (dark gold) motifs in both 3' UTRs (n = 15487 and n = 24972, respectively) and introns (n = 4800 and n = 1519, respectively). Fold enrichments above transcriptomic background are indicated. Peaks containing neither primary nor two or more secondary motifs are shown in grey. 44 four AU-rich 5mers noted above lacked robust peaks and were thus excluded from further analyses (Extended Data Figure 3b). In general, the peaks were stronger and more clearly defined in introns than in 3' UTRs, likely reflecting the predominantly nuclear localization of RBFOX2 and its prominent role in splicing. We generated a CLIP enrichment value for each motif by normalizing the read density at the peak apex to the density at positions distal from the peak (Extended Data Figure 3c). Comparing iCLIP 5mer enrichments in 3' UTRs with nsRBNS 5mer R values in introns (Figure 3b) and 3' UTRs (Extended Data Figure 3d) yielded similar results: GCAUG- and GCAUG-overlapping 5mers were the most highly enriched by both measures, with GCACG next, and the six secondary motifs near the top of remaining 5mers by both measures. We further explored in vivo binding to specific 5mers in neuronal cell types, where expression of Rbfox protein is often very high. Using two high-throughput crosslinking and immunoprecipitation (HiTS-CLIP) datasets—in whole mouse brain and in vitro-differentiated ventral spinal neurons (Jacko et al., 2018; Weyn-Vanhentenryck et al., 2014)—we analyzed the enrichment of secondary motifs in stringently-filtered CLIP peaks present in two replicates in both 3' UTR and intronic regions. As expected, primary motifs were strongly (4.1–6.2-fold) and moderately (1.3–1.7-fold) enriched (Figure 3c-d, Extended Data Figure 4). Notably, four of the six secondary motifs (GCUUG, GAAUG, GUUUG, GUAUG) were also enriched in CLIP peaks to extents approaching that observed for GCACG. Because secondary motifs GUGUG and GCCUG lacked robust enrichment in HiTS-CLIP data, we chose to exclude them from further analysis of in vivo function, focusing instead on the remaining four motifs GCUUG, GAAUG, GUUUG, and GUAUG, all of the form GHWUG. 45 We explored the extent to which secondary motifs can account for CLIP peaks not explained by presence of a primary motif. To do so, we looked for instances of primary, secondary, or both motif types near 3' UTR or intronic HiTS-CLIP peaks that appeared in both replicates (Figure 3e) (Jacko et al., 2018; Weyn-Vanhentenryck et al., 2014). In both datasets, ~50% of peaks contained a primary motif or both types of motifs, while another 15–20% of peaks lacked primary motifs but contained two or more secondary motifs. The presence of ≥ 2 secondary motifs was enriched approximately five-fold in the intronic HiTS-CLIP peaks, and about 2.5-fold in 3' UTR peaks. Thus, consideration of secondary motifs helps to explain a substantial fraction of CLIP peaks that lack a primary motif. Secondary motifs regulate splicing in an Rbfox-dependent manner We next assessed the potential of secondary motifs to mediate Rbfox-dependent regulation. We adapted a bichromatic splicing reporter (Orengo et al., 2006) (pRG6) by introducing one copy of GCAUG or six copies of one of three secondary motifs (GCUUGx6, GAAUGx6, or GUUUGx6) into a 250 nt region downstream of an alternative exon (Figure 4a). For the secondary motifs, the six copies were inserted into interspersed positions (Supplementary Table 4), and a control vector was constructed for each by inserting a permuted version of the motif at the same positions (designated as pGCUUG, etc., with “p” for permuted). Exclusion of the alternative exon yields a transcript in which the DsRED ORF is in frame, while inclusion of the exon shifts the frame so that EGFP is produced; exon inclusion can be measured by red versus green fluorescence. We co-transfected the modified RG6 plasmid with or without a Cerulean:RBFOX1 fusion protein (pCERU:rbFOX1) into HEK293T cells to augment low endogenous levels of RBFOX2. 46 a 250 nt of downstream intron 5mer motifs Permuted 5mer motifs 1. GCAUG.1 (position 1) 1. pGACGU.1 (position 1) x 2. GCAUG.2 (position 2) 2. pGACGU.2 (position 2)3. GCUUGx6 (positions 1-6) 3. pGUGCUx6 (positions 1-6) 4. GAAUGx6 (positions 1-6) 4. pGAGUAx6 (positions 1-6) dsRED EGFP pCERU: 5. GUUUGx6 (positions 1-6) 5. pGUGUUx6 (positions 1-6) Exclusion isoform Inclusion isoform rbFOX1 pRG6 b c 1.00 inclusion exclusion 0.75 Intronic content 0.50 Permuted motif Primary motif Secondary motif inclusion 0.25 exclusion 0.00 .1G . 2 6 6 6 U UG UG x UG x Gx CA CA CU AA U U d G G G G GU 15 10 GCAUG.1, RBFOX GCUUGx6, RBFOX GAAUGx6, RBFOX 5 GUUUGx6, RBFOX GCAUG.1, Null GCUUGx6, Null GpACGU.1, RBFOX GpACGU.1, Null 0 5 <2.2 2.2 3.0 3.0 4.5 4.5 6.0 6.0 7.5 >7.5 RBFOX expression (Cerulean), log2 Figure 4. 47 Inclusion:Exclusion ratio (EGFP:dsRED), log2 pGUGCUx6 pGACGU.1 GCUUGx6 GCAUG.1 pGAGUAx6 pGACGU.2 GAAUGx6 GCAUG.2 pGUGUUx6 GUUUGx6 PSI Figure 4. Secondary motifs in downstream introns promote exon inclusion in an Rbfox- dependent manner in a splicing reporter. a. Experimental design of Rbfox1 splicing reporter. One GCAUG primary motif or six copies of a secondary motif (GCUUG, GAAUG, or GUAUG) were cloned in a 250-base window downstream of an alternative exon in the RG6 dual fluorescent splicing reporter. Plasmids were co-transfected in HEK293T cells with a plasmid expressing fluorescently labelled RBFOX1 to monitor cellular protein levels. b. Semi- quantitative PCR with 5' 6-FAM-labelled primer indicates exon inclusion in the presence of both primary and secondary motifs in an Rbfox1-dependent manner. c. Mean percent spliced in (PSI) values of exons containing primary motifs (gold), secondary motifs (teal), or motif permutations (grey) in the downstream intron after expression of Rbfox1. Error bars show SD of technical replicates in triplicate. For GAAUG, the median permutation value was used due to the introduction of a splicing silencer in its permuted form. d. Per-cell inclusion:exclusion (EGFP:dsRED, y-axis) ratio for the RG6 alternative exon as Rbfox expression (Cerulean, x-axis) increases as measured by flow cytometry. Primary motifs (gold), six copies of three indicated secondary motifs (teal), primary (intact and permuted) and secondary motifs without co- transfection of Rbfox (dark grey), and a permuted primary motif with co-expressed Rbfox (light grey) are shown. Representative data from one of two replicates (the other is shown in Extended Data Figure 6). Bin numbers can be found in Supplementary Table 5. The center line represents the median, lower and upper hinges the first and third quartiles, respectively, and whiskers extend to the smallest or largest value (at most 1.5*IQR (interquartile range) of the hinge. Outliers are not shown. Notches extend 1.58*IQR/sqrt(n), giving a roughly 95% confidence interval on the medians. Uncropped gel image is available as source data online. 48 To measure the effect of Rbfox-mediated regulation of splicing in the presence of secondary motifs, we measured exon inclusion by RT-PCR with a fluorescently-labeled primer (Figure 4b). Insertion of either a single GCAUG or six copies of any of the secondary motifs tested drove percent spliced in (PSI) values of the exon to nearly 100% in the presence of exogenous RBFOX1 (Figure 4c). We observed similar effects of Rbfox protein expression by flow cytometry. Using Cerulean fluorescence to quantify exogenous RBFOX1 levels, and the EGFP:DsRED ratio to quantify exon inclusion, we could measure the relationship between RBFOX1 levels and exon inclusion driven by secondary motifs in single cells (Figure 4d, Extended Data Figure 5-6, Supplementary Table 5). Generally, six copies of any of the secondary motifs GCUUG, GAAUG, or GUUUG drove exon inclusion to about the same extent as a single copy of the primary motif GCAUG. At low levels of RBFOX1 expression (bins 1–2), all reporters exhibited exon inclusion similar to the GCAUG construct co-transfected with an empty Cerulean vector (GCAUG.1, Null), matching the RT-PCR results. As RBFOX1 expression increased, exon inclusion increased well above background for constructs containing either primary or secondary motifs, but failed to increase or increased modestly with control constructs. Across two replicates, exon inclusion was significantly more correlated with Rbfox expression when Rbfox primary or secondary motifs were present (Extended Data Figure 6b). Together, these observations demonstrate that secondary motifs can robustly regulate splicing, particularly when levels of Rbfox proteins are high. Secondary motifs regulate splicing at specific stages of neuronal differentiation 49 Given our observations that secondary motifs are most bound at high Rbfox concentrations, we examined neuronal differentiation, a process in which levels of Rbfox rise naturally. In an eight-stage time course of mESCs differentiating into glutamatergic neurons, total Rbfox expression increases five-fold from the radial glia stage (RG) to developmental stage 3 (DS3) (Figure 5a) (McNutt et al., 2013). During this period, cells progress from a bipolar-shaped progenitor neuron (RG) to fate-specified developmental stage 1 (DS1) neurons, then to developmental stage 3 (DS3) over seven-days, and eventually become mature neurons after extensive growth and pruning of dendrites (Figure 5a). We examined the relationship between exon inclusion and presence of primary or secondary motifs in the first 250 nt of the downstream intron. We first examined the interval RG–DS1, where Rbfox increases from low to moderate expression, and then DS1–DS3, where Rbfox levels reach their highest point. For primary motifs GCAUG and GCACG, exon inclusion is significantly correlated with motif frequency at both intervals. However, the secondary motifs are correlated with exon inclusion only at the DS1– DS3 interval, suggesting that they are mostly active later, when Rbfox levels are higher (Figure 5b). In order to examine potential functions of splicing changes in the DS1–DS3 interval, we assessed the functions of genes whose splicing changed. For exons associated with primary Rbfox motifs (n=388), Gene Ontology analysis found enrichment for functions related to membrane and cytoskeletal organization, while exons associated with secondary motifs (n=561) were enriched for functions in dendrite development and signal transduction (Figure 5c). This distinction could reflect a regulatory program in which primary motifs mediate earlier splicing events related to neurite outgrowth and secondary motifs mediate a later wave of splicing changes related to dendrite development and signaling as Rbfox levels increase. 50 51 Figure 5. Secondary motifs enable splicing regulation at distinct concentration Rbfox concentration ranges in neuronal differentiation. a. Total expression of Rbfox1, Rbfox2, and Rbfox3 in transcripts per million (TPM) based on RNA-seq during a neuronal differentiation time course (McNutt et al., 2013a) (SRA PRJNA185305). b. Correlation of Rbfox primary (gold) and secondary (teal) motifs in the downstream intron with delta PSI at both low-moderate (RG- DS1) and moderate-high (DS1-DS3) transitions of Rbfox expression during in vitro neuronal differentiation. Increased color intensity represents 0, 1, or 2+ motifs in the downstream intron. Primary motifs, RG–DS1: 0 (n = 1521), 1 (n = 269), 2+ (n = 90), DS1–DS3: 0, (n = 3611), 1 (n = 596), 2+ (n = 171). Secondary motifs, RG–DS1: 0 (n = 2763), 1 (n = 838), 2+ (n = 159), DS1– DS3: 0, (n = 6461), 1 (n = 1882), 2+ (n = 413). The center line of the boxplot represents the median, lower and upper hinges the first and third quartiles, respectively, and whiskers extend to the smallest or largest value (at most 1.5*IQR (interquartile range) of the hinge. Outliers are not shown. Notches extend 1.58*IQR/sqrt(n), giving a roughly 95% confidence interval on the medians. *** P < 0.001 (two-sided Wilcoxon Rank-Sum). c. Gene Ontology categories of splicing events driven by primary motifs (n = 388) (top) and secondary motifs (n = 561) (bottom) during neuronal differentiation. Events were compared to all expressed genes at DS3 and terms were filtered for FDR < 0.1, B > 99, and b > 9. d. Correlation of Rbfox primary (gold) and secondary (teal) motifs throughout stages of neuronal differentiation subsequent to RG stage. Pearson correlation of secondary motif presence with delta PSI is shown at intervals of neuronal differentiation from radial glia stage to mature 28-day glutamatergic neurons. Events: RG–DS1 (n = 940), RG–DS3 (n = 3436), RG–MAT16 (n = 3860) RG–MAT21 (n = 3942) RG–MAT28 (n = 4119). Size of point indicates correlation coefficient, intensity indicates uncorrected p-value < 0.05. e. Reporter design of CD47 intron 10-containing RG6. All secondary motifs were ablated from the 617-nt intron and sequentially reintroduced into the reporter to examine the effects of individual secondary motifs. f. Mean percent spliced in (PSI) values of exons containing secondary motifs in the downstream intron after expression of Rbfox1. Error bars show SD of triplicate technical replicates. 52 We next examined the correlation between motif count and difference in exon inclusion (delta PSI) across every interval of the time course (Extended Data Figure 7). The DS1–DS3 interval was the only interval in which secondary motif counts were significantly correlated with Rbfox activity, coincident with peaking of Rbfox expression. Examining intervals from RG to subsequent stages of the time course (Figure 5d), secondary motifs GCUUG and GAAUG had signal throughout the rest of the differentiation, while GUUUG had signal only at DS3, and GUAUG lacked detectable signal. We identified mouse Cd47 exon 10 as a candidate target of Rbfox proteins likely mediated by secondary motifs in its downstream intron. CD47 is a plasma membrane protein with four known isoforms that encode variation in its intracytoplasmic tail (Mordue et al., 2017). Isoforms containing this exon are associated with memory retention in rats (Lee et al., 2000), and CD47- deficient neurons have impaired axon and dendrite formation in development (Murata et al., 2006). Over the course of neuronal differentiation, CD47 exon 10 shows a 36% increase in inclusion from developmental stage 1 (DS1) to developmental stage 3 (DS3). The 617-base downstream intron contains 2 GCUUGs, 2 GAAUGs, and 4 GUUUGs. We replaced the RG6 intron with the natural sequence of the CD47 intron, mutating each secondary motif (Figure 5e, Supplementary Table 4). We successively re-introduced secondary motifs into the intron in 5' to 3' order, measuring the effect of each additional secondary motif on the PSI value of the upstream exon. Restoration of the cluster of five motifs at the 5' end of the intron gave a 19% increase in PSI for the upstream exon. Adding back a subsequent GUUUG and GAAUG yielded additional 10% and 11% increases in PSI, respectively, demonstrating that individual secondary motifs can act additively to trigger substantial increases in PSI (Figure 5f). 53 Rbfox expression-dependent splicing through secondary motifs across neuronal subtypes We next asked to what extent Rbfox contributes to diversification between differentiated neuronal cell types. We analyzed data from thirteen cell types that span almost three orders of magnitude of Rbfox gene expression, combining mRNA levels of Rbfox1, Rbfox2 and Rbfox3 (Weyn-Vanhentenryck et al., 2018). Specifically, we compared mean differences in PSI between cell types grouped by their Rbfox expression from lowest (EC, TRCbitter, OSN25wk) to highest (CGN) (Figure 6a). We expected that the number of GCAUG motifs downstream of alternatively spliced exons should generally correlate with exon inclusion, whereas secondary motifs should contribute only in cell types with higher Rbfox expression. Indeed, we observed that downstream GCAUG motif count is significantly correlated with exon inclusion, even when comparing medium to low Rbfox expression (P < 0.0015), and is more strongly correlated when comparing high to medium (P < 1.3 x 10-19) and highest to low (P < 2.4 x 10-23) (Figure 6b). Consistent with our previous results, the strongest secondary motifs together contribute significantly to increased exon inclusion in high (P < 7.7 x 10-4) and highest (P < 1.3 x 10-4) Rbfox-expressing cells compared to low-expressing cells. When comparing highest to low Rbfox-expressing cells, the slope of the regression predicts that, on average, each downstream GCAUG increases exon PSI by 15%, while each secondary motif elicits a 3.6% increase in PSI. Thus, in endogenous loci as in our reporter assays, the regulatory activity attributable to a single secondary motif is comparable to that conferred by a set of ~4 secondary motifs. We examined 864 alternative exons with increased inclusion associated with Rbfox expression (Extended Data Figure 8). Of these, 11% had primary motifs, 26.4% had primary and secondary motifs, and 3.2% had at least 4 secondary motifs but no primary motif (2.2-fold 54 a comparison 1 comparison 2 comparison 3 b 1. 2. GCAUG 3. 1. 2. GCACG 3. 1. GCUUG 2. GAAUG 3. GUUUGGUAUG r -log10(P) c d Neuronal differentiation 100 1° motifs 2° motifs 10-1 10-2 10-3 10-4 10-2 10-1 100 101 102 tifs Rbfox expression ary mo GCAUG GCAUG Secondaries Non-specific Prim Figure 6. 55 Fraction of Rbfox bound to site Comparison: Cellular diversification Secondary motifs Figure 6. Secondary motifs are active in neuronal cell types with high Rbfox expression. a. Differentiated neuronal cell types from Weyn-Vanhentenryck et al.(2018) (SRA SRP055008), arranged by combined Rbfox1, Rbfox2, and Rbfox3 expression (log10 of sum of RPKM values). Cell types were grouped into Lowest, Low, Medium, High, and Highest Rbfox expression categories. Cells analyzed: olfactory sensory neurons (OMP+) (OSN), enterochromaffin cells (EC), taste receptor cells (TRC), rod or cone photoreceptors, jugular or nodose visceral sensory ganglia, dopaminergic neurons (DN), motor neurons (MN), dorsal root ganglia sensory neurons (Nav1.8+ or Avil+) (DRG), trigeminal ganglia, Purkinje neurons (PN), and cerebellar granule neurons (CGN). b. Linear regression of 1,909 alternatively spliced exons comparing the Rbfox expression groups indicated in a. Horizontal bars represent the Pearson r value (left) and uncorrected significance (right) of the correlation between the number of occurrences of primary or secondary (GCUUG, GAAUG, GUUUG, and GUAUG) motifs and ΔPSI values between Medium and Low (1.), between High and Low (2.), or between Highest and Low (3.) Rbfox expression regimes. Grey dotted and dashed lines indicate 0.05 and 0.01 significance thresholds, respectively. c. An equilibrium model for Rbfox binding to intronic sequences in the nucleus at various expression levels of Rbfox (low, grey area; high, teal area). GCAUG indicated by gold, GCACG indicated by yellow, secondary motifs indicated in teal, non-specific 5mers in grey. d. Graphical summary of how Rbfox proteins (golden ellipses) regulate distinct splicing events at different expression levels (teal shading). Secondary motifs are functionally relevant only at higher Rbfox levels, occurring at later stages of neuronal differentiation or in cell types with high Rbfox expression, while primary motifs are functionally relevant at earlier stages and in cell types with medium as well as high Rbfox levels. 56 enrichment, P < 0.0084, Fisher’s exact test). Even exons with one to three secondary motifs (n = 354) were 1.3-fold enriched in this set (P < 0.0012). This analysis supports that secondary motifs contribute to regulation of dozens of exons or more in neuronal subtypes with high Rbfox levels. A quantitative model for Rbfox expression-dependent, differential motif activity In order to better understand the distribution of Rbfox protein binding across the nuclear transcriptome, we built a quantitative equilibrium model (Extended Data Figures 9, 10). To illustrate the plausible range of RNA concentrations in the neuronal nucleus, we considered two scenarios: a large cell with lower RNA turnover and lower nuclear RNA concentration (Figure 6c, Extended Data Figure 10a,b) and a small cell with higher RNA turnover and nuclear RNA concentration (Extended Data Figure 10c,d) (Methods). Using this model, we estimated the aggregate concentration of neuronal intronic RNA 5mers as between ~14 and 56 μM (Methods). Of this total, the concentration of GCAUG will fall between ~16 to 63 nM, with GCACG five-fold lower (3–12 nM), and the four secondary motifs together occurring at several-fold higher concentration (67 to 270 nM) (Extended Data Figure 10a). After assigning approximate dissociation constants to all Rbfox 5mer motif variants by calibrating our random RBNS data for human RBFOX2 and RBFOX3 (Dominguez, et al., 2018) to SPR data (Stoltz, 2015) (Extended Data Figure 9, Supplementary Table 6), we modeled the equilibrium distribution of Rbfox proteins across the pool of nuclear binding sites, analogous to a previous model for miRNAs in the cytoplasm (Jens & Rajewsky, 2015), as a function of nuclear Rbfox concentration (Figure 6c). Next, we estimated the nuclear Rbfox concentration from the range of Rbfox mRNA expression across neuronal cells (Figure 6a), extrapolating from the mRNA expression, using proteins-per-mRNA ratios of 2,800–10,000, in line with reported 57 protein:mRNA ratio for Rbfox in cell lines (Schwanhüusser et al., 2011a; Wiśniewski et al., 2014). The model indicates that, in cells with low to intermediate Rbfox expression, 20% to 40% of Rbfox is recruited to primary motifs, even though they represent a small fraction of the available RNA pool. On the other hand, at very high Rbfox expression, the primary motifs become saturated (because of their high affinity and low abundance) and a larger fraction of Rbfox protein occupies secondary motifs (Figure 6c, Extended Data Figure 10d). We estimate that individual instances of secondary motifs reach ~50% occupancy at nuclear Rbfox concentrations between 10 and 50 μΜ (Extended Data Figure 10b,c). To support this model, we performed filter binding experiments wherein radiolabeled primary (GCAUG) or secondary (GCUUG, GAAUG, GUUUG) motifs in three copies were co- incubated with RBFOX2 and unlabeled, single-copy GCAUG at increasing protein concentrations (Extended Data Figure 10f). As predicted, the fraction bound increased for both primary and secondary motifs, with increases occurring at somewhat higher protein concentrations for secondary motifs, mimicking the modeled scenario in which secondary motifs are preferentially bound after saturation of primary motifs. We conclude that the secondary motifs identified here contribute to regulation under physiological conditions, with greatest activity at high nuclear concentrations of Rbfox proteins (Figure 6d). Insights into concentration-dependent regulation via RBP suboptimal motifs We show that Rbfox family proteins bind a defined, abundant set of secondary motifs in vitro and in vivo. These motifs enable concentration-dependent regulation of exon inclusion by Rbfox in neuronal differentiation and diversification, adding a second wave of regulation in cells that with high Rbfox expression. The mediation of temporally-staged and cell-type-specific activity 58 by secondary motifs has rarely been reported for RBPs (Wagner et al., 2016; Xiao et al., 2009), but may be widespread. Examples of spatial and temporal regulation via secondary motifs have been described for several transcription factors (TFs). The C. elegans TF PHA-4 activates its primary motif target in early pharyngeal development, while its secondary targets are activated later, when protein levels are higher (Gaudet & Mango, 2002). In mice, PREP1 activates its enhancers sequentially, according to their affinity, during eye lens development (Rowan et al., 2010), and suboptimal motifs for the GATA and FGF families of TFs enable tissue-specific expression in tunicates (Farley et al., 2015). At the genomic level, a Bayesian biophysical model that considered binding to both high- and low-affinity sites better predicted gene expression in human cells, supporting the importance of low-affinity sites to regulation (Wang et al., 2015). The Rbfox secondary motifs identified in this paper likely represent a combination of the strongest-binding and most frequently-occurring secondary motifs across a continuum of related sequences that bind RBFOX with varying affinity. One feature of the RBNS natural sequence assay, as opposed to a random sequence RBNS experiment (or SELEX), is that its sequence composition and oligonucleotide length aids in the discovery of motifs that appear frequently in the transcriptome in arrangements compatible with binding. Therefore, while the discovered motifs may not represent all of those with highest affinity, they should tend to be those with largest transcriptome-wide regulatory effects. As an in vitro assay, nsRBNS in its current form does not capture the native protein environment of the cell, with interacting proteins and modifications, for example, but is a straightforward method to determine the spectrum of protein binding to cellular sequences. 59 Though variant motifs have occasionally been noted in other Rbfox studies, many go on to filter out sequences lacking GCAUG from their datasets (Gehman et al., 2011, 2012; Jacko et al., 2018; Lee et al., 2016; Vuong et al., 2018; Weyn-Vanhentenryck et al., 2014; Zhang et al., 2008). In some cases, binding of secondary motifs has been attributed to other proteins, e.g., presence of GUGUG in Rbfox CLIP peaks has been attributed to binding by SUP-12 (Kuwasako et al., 2014) or HNRNP M of the LASR complex (Damianov et al., 2016). The potential for direct binding of GUGUG motifs by Rbfox proteins should also be considered. Our finding that secondary motifs exert regulatory activity only in cell types with high Rbfox levels may explain varying enrichments of primary and secondary motifs among CLIP data derived from different cell lines (Damianov et al., 2016; Jacko et al., 2018; Weyn-Vanhentenryck et al., 2014). Suboptimal motifs may function in different ways. First, they may enhance binding to nearby primary motifs. Second, a cluster of several secondary motifs may create a stretch of sequence with total Rbfox binding comparable to that of a high-affinity motif. Consistent with this idea, we observed increased binding to oligonucleotides with larger numbers of secondary motifs in nsRBNS, and found evidence that 4–6 copies of a secondary motif can enhance exon inclusion to an extent comparable to that of a single GCAUG motif in vivo. Third, these motifs can introduce more subtle changes in exon splicing, as we saw for CD47, or incrementally tune splicing over evolutionary time. Lastly, suboptimal cis-regulatory elements may function only at high levels of a trans-factor to narrow the temporal or spatial scope of activity, a phenomenon we observed for Rbfox regulation in cellular differentiation and diversification. Our study argues that it is time to reconsider the simplifying dichotomy of ‘binding site’ versus ‘non-binding site’ and to start incorporating site affinities as well as motif transcriptomic frequencies, extending consideration to abundant sites of lower affinity when RBP activity is 60 high (Jankowsky & Harris, 2015). For Rbfox family proteins, we demonstrate that secondary motifs contribute significantly to Rbfox-dependent gene regulation, with their modest affinity enabling highly cell- and stage-specific regulation. This behavior could also be relevant in regimes of high local protein concentration. For instance, Rbfox and other RBPs are associated with phase-separated compartments, which rely on multivalent protein and RNA interactions to form (Gomes & Shorter, 2019; Sanders et al., 2020; Ying et al., 2017); secondary motifs will have increasing significance as our understanding of such subcellular environments deepens. 61 Acknowledgements We thank past and present members of the C. Burge Lab, J. Conboy, and I. Jarmoskaite for helpful comments on the manuscript. We gratefully acknowledge the courtesy of the C. Zhang lab of Columbia University, who shared intermediate results from Weyn-Vanhentenryck et al. (2018) used for Figure 6 (gene expression, PSI values, exon coordinates). M. Jens received an EMBO Long Term Fellowship ALTF-1130-2015. All other authors were supported by NIH 5- R01-GM085319. Data availability nsRBNS raw data is available under accession number GSE152510 and processed data is available in Supplementary Table 2. Due to their large volume, FACS data are available from the corresponding author on reasonable request. Data used in other analyses can be found at PDB 2ERR (Figure 2a), GEO GSE54794 (Figure 3a,b), SRA SRP128054, SRP035321 (Figure 3c-e), SRA PRJNA185305 (Figure 5), and SRA SRP055008 (Figure 6a-b). 62 Materials and methods Cloning, expression, and purification of RBFOX2 The RRM domain of RBFOX2 (amino acids 100-194) was cloned into the pGEX6P-1 expression vector (GE Healthcare, 28-9546-48) downstream of a GST-SBP tandem affinity tag. Following 12-hr recombinant expression in Rosetta Competent Cells (Millipore, #70954) at 12ºC with ampicillin and chloramphenicol selection, the protein was expressed by addition of IPTG and purified via the GST tag as described previously (Dominguez, et al., 2018; Lambert et al., 2014; Taliaferro et al., 2016). Library design of natural 3' UTR sequences GENCODE59-annotated human transcripts were evaluated for presence of appropriate start codon and stop codon sequences at the corresponding GENCODE-annotated positions. Coding genes that had a 3' UTR between 100 and 10,000 bp were considered for library inclusion. For the H. sapiens library, 120 transcript pairs were included because of evidence of alternative 3' UTRs, 720 transcripts were selected based on expression level of >10 fragments per kilobase million (FPKM) in both HepG2 and K562 cell lines, and 360 other transcripts were selected at random. Each of the resulting transcripts was assigned a homologous M. musculus transcript by identifying a transcript in the homologous M. musculus gene in which the annotated polyA tail was ±150 nt from the location of the homologous H. sapiens polyA tail site, as identified by Batch Coordinate Conversion (liftOver, UCSC Genome). This procedure generated 1108 H. sapiens 3' UTRs paired with 1104 M. musculus 3' UTRs. Each 3' UTR sequence was then split into overlapping 110 nt segments at 43 nt intervals to achieve approximately 2.5X coverage of each 3' UTR, yielding 64,319 unique sequences. Natural sequence RNA Bind-n-Seq procedure and analysis The 64,319-oligonucleotide library was synthesized by Twist Biosciences and nsRBNS was performed as previously described (Dominguez, et al., 2018; Lambert et al., 2014; Taliaferro et al., 2016). Briefly: Library was amplified with Phusion Polymerase (NEB, #E0553L, primers below, Integrated DNA Technologies (IDT)), in vitro-transcribed, treated with Turbo DNase (Thermofisher, #AM2238), and gel- and phenol:chloroform:isoamyl alcohol-purified. Streptavidin T1 magnetic beads (Invitrogen, #65601) and 250 nM of the prepared library was incubated with recombinant, tagged RBFOX2 at concentrations of 0 (no protein control), 4, 14, 43, 121 (2 replicates), 365, and 1100 nM for 1 hr at 4º C to equilibrium binding in binding buffer (25 mM Tris, 150 mM KCl, 0.1% Tween, 0.5 mg/ml BSA, 3 mM MgCl2, 1 mM DTT, pH 7.5). RBP:RNA:bead complexes were pulled down with a magnet and washed gently twice with wash buffer (25 mM Tris, 150 mM KCl, 0.1% Tween, 0.5 mM EDTA, pH 7.5). RNA was eluted by two separate incubations with 4 mM biotin (pH 7.5) for 30 min at room temperature. RNA eluate was purified with Ampure beads (Beckman Coulter, #A63987), and the resulting RNA was reverse transcribed with Superscript III (Thermofisher, #18080093). The cDNA was amplified for 6-16 cycles with Phusion Polymerase (NEB, #E0553L) and gel-purified with ZymoClean Gel DNA Recovery Kit (Genesee Scientific, #11-300C) to produce the final library. Each library was sequenced single-end on an Illumina HiSeq 2500 instrument. Sequences with at least 100 associated reads in the input sample were considered for further analysis. Read counts in each sample were normalized by the total reads in that sample. To produce enrichment (R) values at the single-oligonucleotide level, normalized reads for an oligonucleotide in a given sample were 63 divided by the normalized reads for that oligonucleotide in the 0 nM control sample. In general, sequences enriched at an R value of 1.1 (10% enrichment) were considered “bound”. Replicates of nsRBNS performed with different in vitro transcriptions and on different days correlate at r > 0.9, although there is variability at finer intervals of R values. To reduce technical noise, we have found that it is important to keep the number of PCR cycles after RNA pulldown as close as possible across both experiments and samples. For example, the 0 nM control for the first 121 nM experiment was cycled 9X, while the second was cycled 16X, which may have caused variation due to different amplification biases, expected to be more pronounced for lowly enriched motifs. Approximately 2,700 oligonucleotides had R ≥ 2.0 in these experiments, with an estimated FDR of 20% at this cutoff. PhyloP scores Plus-strand hg19 46-way alignment phyloP (Pollard et al., 2010) scores were obtained for each genomic position from the UCSC genome database. Hg19 46-way alignment scores are not strand-symmetric (https://genome.cshlp.org/content/suppl/2009/10/27/gr.097857.109.DC1/supplement.pdf, S2.6), and minus-strand scores are thus inappropriate to include in transcriptomic analyses and were excluded. For meta-phyloP scores, all scores for the sequence GCAUG were averaged at each base in 3' UTRs and shallow introns (+250 bases) and normalized to the average phyloP score of the 5mer. Identification of secondary motifs by iterative kmer analysis Initial 5mer enrichments were determined by generating five-base sequences from the 3' UTR regions of each valid oligonucleotide in a 0 nM and pulldown experiment. 5mer enrichment values were then determined by dividing the normalized count of a particular 5mer in the pulldown by the normalized counts of the 5mer in the 0 nM experiment. Of the 1024 5mers, GCAUG had the highest R value. To identify other sequences that influence binding, we took an iterative approach in which all oligonucleotides containing the highest R value 5mer were removed from consideration, and enrichments for all other 5mers were regenerated from this set. Following identification of the primary motifs (GCAUG, GCACG) in the 1.1 µM experiment, the top 5mer in each subsequent iteration was considered a secondary motif (GUUUG, GAAUG, UUUUU, UGCAU, AAAAA, GUAUG, GCUUG). After four iterations, there was insufficient power to continue the iterative method with ~18,000 of ~64,000 sequences remaining. All remaining 5mers of the format GNNUG, NGCAU, or W5 that had R value two standard deviations above the mean were considered secondary motifs (GCCUG, GCCUG, CGCAU, AGCAU, AUAUA, UAUAU). Filter binding 23-base oligonucleotides with three copies of motifs of interest spaced by two random bases (IDT, Supplementary Table 3) were in vitro-transcribed as above and 5' dephosphorylated with Calf Intestinal Phosphatase (NEB, #M0290) according to manufacturer’s protocol and phenol:chloroform:isoamyl alcohol-purified. ~120 ng RNA was radiolabeled with T4 Polynucleotide Kinase (NEB, #M0201L) and 2 mCi [γ-32P]-ATP (PerkinElmer). Unincorporated [γ-32P]-ATP were removed with illustra Microspin G-25 Columns (GE Life Sciences, #27532501). RBFOX2 was purified as above and buffer exchanged using a Zeba Spin Desalting Column (7K MWCO, 0.5mL, Thermofisher, #89882). 20uL reactions were prepared in 64 96-well plates with binding buffer (10% glycerol, 25nM Tris, 150mM KCl, 0.1mg/mL BSA, 1mM MgCl2, 1mM DTT, pH 7.5), 1-5nM radiolabeled RNA, and six concentrations of protein: 25uM, 5uM, 1uM, 200nM, 40nM, and 8nM. After 1h incubation at RT, 10uL of the reactions were applied to pre-soaked (wash buffer, 25nM Tris, 150mM KCl, 1mM MgCl2, pH 7.5) stacked 0.45 µm nitrocellulose (Thermofisher, #77010) and nylon (Amersham Hybond-XL, GE Life Sciences, #RPN303S) in a 96-well Bio-Dot vacuum apparatus (Bio-Rad, #RPN303S) on low vacuum. Reactions were immediately washed with 100µL wash buffer. Blots were exposed on a phosphor-screen (GE Healthcare) and imaged on a Typhoon FLA 9500 (GE Healthcare). Fraction bound was quantified using ImageJ software. For competition filter binding experiments in Extended Data Figure 10, roughly equimolar unlabeled, single copy GCAUG was additionally added to each reaction and assay was otherwise as above. For these experiments, due to an increased background with the higher concentration of RNA, the poly-U control values were subtracted from our measurements. eCLIP peak enrichment and iCLIP metaplot analyses eCLIP enrichment values were produced from significant RBFOX2 eCLIP read peaks in HepG2 cells obtained from the ENCODE Project Consortium (Feingold et al., 2004). RBFOX2- pulldown peak densities were normalized to a no-protein input control to produce enrichment values roughly analogous to nsRBNS R values to facilitate comparisons between the two assays. RBFOX2 individual nucleotide crosslinking and immunoprecipitation (iCLIP data) from mouse embryonic stem cells (Jangi et al., 2014) was analyzed at RBFOX secondary motif sites. Adapters and barcodes were trimmed prior to mapping with STAR to the mm10 genome following standard ENCODE guidelines (http://labshare.cshl.edu/shares/gingeraslab/www- data/dobin/STAR/STAR.posix/doc/STARmanual.pdf, page 7). Duplicate PCR reads were removed from the mapped reads to generate final reads. These reads were aligned in a metaplot centering on all possible 5mers to visualize an iCLIP meta-peak. Peak height was quantified with a CLIP enrichment (CE) score centered on position 1 of the 5mer. For 3' UTR peaks, the CE score was calculated as the sum of the read coverage between positions 10 to 15 divided by the sum of the read coverage between positions –85 to –80. For the intronic peaks, the CE score was calculated by the read coverage between positions 5 and 10 divided by the read coverage at positions -85 to -80. The ranges differed between 3' UTR and intronic peaks due to different maximum peak heights in these regions; the maximum range was chosen for each region. Control scores were produced using an untagged RBFOX2 protein with identical data processing; these background reads were subtracted from the metaplots and CE scores to eliminate iCLIP noise. Total base-read counts per regions were 1) 3' UTR: nGCAUG = 8296148, nGCACG = 914837, nGCUUG = 3876531, nGAAUG = 3529113, nGUUUG = 6172128, nGUAUG = 3482546, nGUGUG = 8735260, nGCCUG = 5432910, nUUAUG = 3954351, nAUAUG = 3142196, nGUAUU = 4395833, nGUAUA = 3092874; 2) Intronic: nGCAUG = 10538122, nGCACG = 1316449, nGCUUG = 3649606, nGAAUG = 2319706, nGUUUG = 4261360, nGUAUG = 2624047, nGUGUG = 7388214, nGCCUG = 5971647, nUUAUG = 2503772, nAUAUG = 2051322, nGUAUU = 2194441, nGUAUA = 1675246. HiTS-CLIP motif enrichment Two RBFOX1 HiTS-CLIP (Jacko et al., 2018; Weyn-Vanhentenryck et al., 2014) datasets, in mouse whole brain and differentiated mature neurons, respectively, were analyzed for the presence of secondary motifs. Both datasets were mapped with STAR after removing duplicate reads following ENCODE guidelines, except for requiring mapped reads to be completely 65 unique (--outFilterMultimapNmax 1). For the Jacko et al. (2018) data, read length was relaxed to accommodate the slightly shorter average HiTS-CLIP read length (--outFilterMatchNminOverLread 0.33). Mapped reads were then processed into CLIP peaks using CLIPper with standard specifications. For each dataset, only peaks shared between two replicates were considered in subsequent analyses. Additionally, only peaks in 3' UTRs and shallow (splice site-proximal) regions of introns were analyzed. For each CLIPper peak, a region ±50 from the reported peak apex was analyzed for the frequencies of all possible 5mers to report a total 5mer frequency for each dataset analyzed. 5mer frequencies in CLIPper regions were normalized to the total 5mer frequencies in either 3' UTRs and shallow introns to generate 5mer enrichments. To analyze the enrichments of peaks containing primary, primary and secondary, or secondary motifs, a peak with one of these motifs in the 100-base region was considered to contain the motif. A minimum of two secondary motifs was required in this analysis to adjust for the higher frequency of these motifs in the transcriptome. To assess signal over background for these groups, artificial 100-base intervals derived from 3' UTRs and shallow introns, respectively, were analyzed identically and frequencies were compared. Splicing reporter assay We cloned primary and secondary Rbfox motifs GCAUG.1, GCAUG.2, GCUUGx6, GAAUGx6, or GUUUGx6 250 bases downstream of the alternative exon of the RG6 splicing reporter (Orengo et al., 2006) (see sequences below, with altered nucleotides in capitals) using custom-designed oligonucleotides (IDT) with InFusion cloning (Takara Bio #638920) in HEK293T cells. The GFP of a pEGFP rbFOX1 plasmid (Addgene #63085) was replaced with Cerulean (Cerulean-N1 Addgene #54742) to produce a Cerulean:Rbfox1 vector. The downstream Rbfox1 was also removed to produce a Cerulean:NULL control plasmid (see Supplementary Table 4). A far-downstream GCAUG endogenous to the plasmid was included in all plasmids. We also introduced variation on a natural intron, mouse CD47 intron 9, into the RG6 plasmid. A wild-type and secondary motif-null construct were synthesized with GeneWiz FragmentGene, and sequential mutations in or restorations of secondary motifs were introduced using the QuickChange Lightning Multi Kit (Agilent, #210515) and custom-designed oligonucleotides (IDT) to produce a series of introns with increasing numbers of secondary motifs in the intron (see Supplementary Table 4, with altered nucleotides in capitals). 100 ng of RG6 construct and 300 ng of Cerulean:Rbfox or Cerulean:NULL were co-transfected into HEK293T cells in a 24-well plate with 1 uL of lipofectamine. Transfected cells were harvested after 24h. The cells were washed twice with 2 mL PBS and RNA was extracted using the Qiagen RNAeasy Kit (#74104). Three replicates of each condition were subjected to fluorescent PCR with a FAM-labelled forward primer (below) and with Phusion polymerase (NEB #M0530S) for 32 cycles. The product was imaged on a Typhoon FLA 9500 (GE Healthcare). Resultant bands were quantified using ImageJ to produce relative Percent Spliced In (PSI) values. Flow cytometry and data processing 400 ng of RG6 construct, 400 ng of Cerulean:Rbfox or Cerulean:NULL was transfected into HEK293T cells in a 24-well plate with 4 uL of lipofectamine. Transfected HEK293T cells were harvested after 48 hours of transfection in 6-well plates, with 106 cells per well. After the media was removed, cells were gently washed in 1 mL PBS, and then resuspended into 1 mL ice-cold 66 PBS with 1% BSA and 2 mM EDTA. Cell suspensions were collected into test tubes through a single-cell strainer (Fisher Scientific; Corning Falcon #352235) on ice. Flow cytometry was carried out with the LSR II flow cytometer (BD Biosciences), with the 405 nm laser and 450/50 nm filter for Cerulean, 488 nm laser and 515/20 nm filter for EGFP, and 561 nm laser and 610/20 nm filter for DsRED. A total of 30,000 events from single, live cells were acquired for each treatment set, and processed using the FlowJo software. Minor channel spillover between EGFP and Cerulean was compensated using single-fluorophore controls. Cerulean signal was normalized by dividing over the median Cerulean signal of the no-plasmid control, to account for background fluorescence. To only include cells with at least one copy of both plasmids transfected, we used the 99th percentile of the signal from the three respective fluorophores in the no-plasmid control as the thresholds, and filtered for events with Cerulean above threshold, and with EGFP and/or DsRED above threshold. Events with EGFP or Cerulean signal above 104.5, or DsRED signal above 103.2, were discarded due to signal anomaly near the detector saturation limit. Events were sorted into bins by log2-transformed normalized Cerulean signal, dividing at 2.2, 3.0, 4.5, 6.0, and 7.5, to obtain six bins with generally similar number of events each. Log2-transformed ratio of EGFP signal to DsRED signal was used as the readout for the splicing ratio. Boxplots were visualized using ggplot2 (geom_boxplot). The center line represents the median, lower and upper hinges the first and third quartiles, respectively, and whiskers extend to the smallest or largest value (at most 1.5*IQR (interquartile range) of the hinge. Outliers are not shown. Notches extend 1.58*IQR/sqrt(n), giving a roughly 95% confidence interval on the medians. Bin numbers in Supplementary Table 5. Cell lines HEK293T-A2 were obtained courtesy from the Eugene Makeyev Lab at Nanyang Technological University, Singapore. Cells were tested for mycoplasma by PCR. Cells were not authenticated. Neuronal differentiation analysis Using a deep transcriptomic sequencing dataset characterizing neuronal differentiation from mouse embryonic stem cells (mESCs) to glutamatergic neurons (McNutt et al., 2013), we examined the use of Rbfox primary and secondary motifs in neuronal differentiation. Each of eight (ESC, NESC, RG, DS1, DS3, MAT16, MAT21, MAT28) time points was analyzed for total Rbfox (Rbfox1, Rbfox2, Rbfox3) expression using kallisto (Bray et al., 2016) with standard parameters. Splicing events were analyzed with rMATS.4.0.2(S. Shen et al., 2014) using standard specifications between the radial glia (RG) stage and all subsequent events (e.g. RG– DS1, RG–DS3, RG–MAT16, etc.) as well as between each interval (e.g. ESC–NESC, NESC– RG, RG–DS1, etc.). For all significantly changing cassette exons (FDR < 0.1), the secondary motif content of the first 250 bases of the downstream intron was computed and correlated with the magnitude of the inclusion of the upstream exon as reported by rMATS. Boxplots were visualized using ggplot2 (geom_boxplot). The center line represents the median, lower and upper hinges the first and third quartiles, respectively, and whiskers extend to the smallest or largest value (at most 1.5*IQR (interquartile range) of the hinge. Outliers are not shown. Notches extend 1.58*IQR/sqrt(n), giving a roughly 95% confidence interval on the medians. Gene Ontology Genes regulated exclusively by primary motifs (primary motif-mediated) or secondary motifs (secondary motif-mediated) in the 250 nt of the intron downstream of an exon increasing in 67 inclusion from DS1 to DS3 were subjected to Gene Ontology analysis with GOrilla (Eden et al., 2007, 2009). Results were then filtered by FDR < 0.1, B > 99, and b > 9. Background genes were all genes expressed >1 transcript per million (TPM) in DS3 as assessed by kallisto (Bray et al., 2016). Linear regression analysis of motif frequency and alternative splicing Starting from a table with PSI values for 1,909 alternative exons regulated in neuronal cells and a list of RBP expression values (courtesy of the Chaolin Zhang lab, Columbia University), underlying data (Weyn-Vanhentenryck et al., 2018), cell types were grouped by the sum of Rbfox1, Rbfox2, and Rbfox3 expression values. For each group of cell types, the arithmetic mean of PSI values was computed, from which followed a ΔPSI value for changes in exon inclusion between Low and Medium, Medium and High, and Medium and Highest Rbfox- expressing cells. Next, for each alternative exon, the intronic sequences 8nt to 250 nt downstream of the exon were scanned for the presence of Rbfox motif 5mers. Then, for each vector of ΔPSI values across all exons, linear regression was performed with the number of Rbfox motif occurrences as explanatory variable (using the python scipy.stats.linregress function of the scipy package (Jones et al., 2001)). The resulting P-values and r values were plotted using matplotlib (Hunter, 2007). Furthermore, PSI values were converted to logit scores by computing logit(PSI) = log2 PSI/(1- PSI) after replacing 1 with 0.999 and 0 with 0.001 to avoid infinities. A regression was then performed to find optimal loadings for each considered primary and secondary motif, as well as scrambled negative control motifs, to explain the PSI values observed at high Rbfox expression, as a function of the PSI values observed at low Rbfox concentration and a linear combination of motif loadings for motifs present in the considered region of the downstream intron (see above), after the mean trend had been subtracted. Briefly, we predicted PSI2' = 1/(1 + 2x) - <ΔPSI> with xi = logit(PSI1)i + Σj fi Mij, where M is the matrix of occurrences of motif j for each alternative exon I, and <ΔPSI>= the mean change in exon inclusion across all alternative exons. The regression then computed optimal motif weight fj to minimize the squared difference between predicted and observed PSI: fjopt = argmin { (PSI2' (fj) – PSI2)2 } . Exons with no primary or secondary motifs in the downstream intron +8 to +250 window were dropped from the regression. We performed this regression on 1,000 bootstraps of the alternative exon set (random sampling with replacement) to assess the robustness of the resulting motif loadings. RBNS calibration to surface plasmon resonance (SPR) database We computed RBFOX2 and RBFOX3 7-mer enrichments from RNA Bind-n-Seq (RBNS) experiments performed with 1.1 and 1.3 μΜ respectively (Dominguez, et al., 2018; Lambert et al., 2014). These R-values should directly track with the occupancy of RBFOX protein, but also contain a contribution from non-specific sequences, either captured through the apparatus or due to the fact that 20 nt and 40 nt random sequences were used rather than 7mers (see Lambert et al., 2014 for a discussion). We therefore estimated the R-value of non-specific 7mers Rns from 68 the bottom percentile of R-values and derived corrected R-values as R' = R + Rns * (R – 1)/(1 – Rns). The corrected R-values display a higher dynamic range and provide a slightly better fit to SPR reference affinities (not shown). The results change only marginally if Rns is estimated from other percentiles (5 or 10, not shown). Next, we compiled a list of dissociation constants for 22 7mer sequences binding to RBFOX1, measured via SPR from Auweter et al. (2006) and Stoltz (2015). We then used the scipy.stats.linregress function to find an optimal linear relationship between log(R’) and log(Kd) for these sequences. We then extended the linear interpolation to all 7mer R’ values to assign approximate dissociation constants. The two experimental replicates (with RBFOX2 and RBFOX3, using 20 nt and 40 nt random sequences) yielded highly correlated results (Extended Data Figure 10c). We then computed average 5mer dissociation constants as Kd5 = exp( ), where < … > denotes the arithmetic mean over all 7mers containing the 5mer of interest, and across both calibrated data sets for RBFOX2 and RBFOX3. Values obtained by subtracting or adding one standard error of the mean were used to estimate the error bars. In the affinity histogram (Extended Data Figure 10d), non-primary or secondary 5mers containing partial primary motifs (GCA, AUG, and ACG) were masked, because RBNS R values are always contaminated by the enrichment of kmers that overlap authentic high affinity motifs but are not necessarily directly bound (see Lambert et al., 2014). Estimation of the intronic sequence content of the nucleus Common estimates for total mRNA copy numbers per cell range from 100,000 to 1,000,000 molecules per mammalian cell (consistent with 0.1 to 1 pg of mRNA per cell) (BioNumbers.org ID 111220). This is in line with previously estimated mRNA copy number from human CA1 pyramidal neurons (~1,000,000) and the smaller cell body size of mouse pyramidal neurons (Benavides-Piccione et al., 2019; Kosik, 2016). To arrive at a rough estimate of how much intronic RNA is made, we therefore convert the estimated average half-life time of mRNA molecules of Th = 5 to 10 hours into a rate kdeg = log(2)/Th and assume that mRNA decay is balanced by new mRNA production (BioNumbers.org ID 106378, 104747), consistent with measured mRNA half-lives for mouse embryonic stem cells treated with RA and LIF withdrawal (Sharova et al., 2009). We thus considered two scenarios that span the anticipated range of RNA concentrations based on these estimates: a small cell (concentrated RNA) scenario with 100,000 mRNAs in a cell volume of 500 μm3 with Th = 5 h; and a large cell (dilute RNA) scenario with 1,000,000 mRNAs in a cell volume of 10,000 μm3 with Th = 10 h, with the cell sizes estimated using a rodent neuron (BioNumbers.org ID 112112, ID 106320). We further assume that the nucleus comprises ~10% of the cell’s volume, and that effectively 70% of that volume is available for the diffusion of intronic RNA and Rbfox proteins (subtracting nucleolus and chromatin). Our findings are not dependent on these exact values. By iterating over the catalog of mouse transcripts expressed in differentiating mouse neurons (McNutt et al., 2013) (Gencode M97) and weighting each encountered intron with the RNA-seq- derived TPM (transcripts per million) expression value for the harboring transcript (using kallisto (Bray et al., 2016)) we thus estimate the amount of intronic sequence being transcribed per minute. Assuming an average intronic half-life time of 1 minute, we conclude that, on average a mouse neuronal cell nucleus may harbor between 1.2 and 6 million intronic 5mer sequences that could represent potential protein binding sites. Using the same TPM-weighting scheme, we then compute the share of this total nuclear, intronic 5mer concentration allotted to each 5mer. 69 With a given estimate of total cellular mRNA copy number, TPM values correspond directly to cellular mRNA copies/cell. To estimate how many protein molecules are present per mRNA copy, we investigated Schwanhüusser et al. (2011) who report ~10,800 proteins/mRNA for Rbm9 (aka Rbfox2) in NIH 3T3 cells and a median of ~2,800 proteins/mRNA across all proteins. Li et al., (2014) (BioNumbers.org ID 110236) suggest that this might be an underestimate and place the median at 9,800 proteins/mRNA. We also investigated Wiśniewski et al. (2014) who report ~74,000 copies of Rbfox protein per A549 cell. Combined with mRNA expression data from the EBI Gene Expression Atlas (E-MTAB-4729) of 37 TPM and an estimated mRNA count of 300,000 mRNAs per cell, this yields an estimated ratio of 6,700 proteins/mRNA. We therefore consider the range between 2,800 (a lower estimate for the median across all proteins) as a reasonable lower bound for Rbfox and 10,000 (reported by Schwanhüusser et al. (2011) for 3T3 cells) as reasonable Rbfox protein to mRNA ratio, while assuming that the bulk of Rbfox protein is nuclear. This range reflects the shaded areas in Figure 6c corresponding to low Rbfox expression (10 TPM) and the highest Rbfox expression (1,900 TPM) observed among the neuronal cell types. Equilibrium model for protein binding to a diverse pool of binding sites We follow the procedure employed in Jens & Rajewsky (2015) for miRNA binding sites. Briefly, in equilibrium, within a well-mixed compartment, every potential binding site interacts with the same pool of free (unbound) protein [Pfree]. Further, each binding site’s occupancy is assumed to be solely dependent on its primary sequence affinity, represented by a 5mer (no cooperativity, no competition with other proteins or RNA structure). Under these assumptions, we can compute the occupancy for each 5mer i using its SPR-calibrated, RBNS-derived Kd (see above) as Oi = [Pfree] / ([Pfree] + Kd) . This is the familiar Michaelis-Menten type relationship between protein concentration and binding probability, also known as the Langmuir isotherm. The distinctive features are approximate linearity for [Pfree] << Kd where Oi ≈ [Pfree]/Kd, Oi = 0.5 for [Pfree] = Kd, and saturation with Oi asymptotically approaching 1 for [Pfree] >> Kd. The relationship between total and free protein can be found from mass-action: [Pfree] = [Ptotal] – [Pbound]. And [Pbound] = Σi Oi * ci (where ci is the estimated concentration of 5mer i in the nucleus). This is equivalent to standard formulations in biophysical chemistry for mixtures of many ligands (Cox, 1981). Because [Pfree] occurs on both sides of this non-linear equation, we numerically find the free RBP concentration 0 < [Pfree] < [Ptotal] to optimally satisfy mass-action. From this then directly follow the expected occupancies of Rbfox motifs and the fraction of total Rbfox protein allotted to the corresponding sites as f = Σj Oj * cj / [Ptotal] (where j are motif 5mers). Non-primary, non-secondary 5mers with partial overlap to primary motifs (139 out of 1,024 5mers) were ignored for this analysis. Statistics Supplementary Table 7 contains details of all statistical tests performed. 70 Supplementary figures. 71 72 Figure S1. RBFOX2 nsRBNS reveals binding to a set of moderate-affinity secondary motifs. a. Correlations among seven natural sequence nsRBNS experiments. Pearson correlations are reported for any sequence with an enrichment (R) value greater than 1. Darker color indicates a higher correlation (R 1.1.463 cor.test function). n = 38467. b. Correlation of nsRBNS R with eCLIP enrichment at oligo-derived regions for all oligonucleotides or sequence regions containing a single GCAUG Rbfox primary motif (n = 2946). c. R value distribution of nsRBNS sequences containing 0 (n = 21596) or 1-3 (n = 2397) NGCAU motifs. d. R value distribution of nsRBNS sequences containing 0 (n = 11077) or 1-3 (n = 12916) AU motifs. e. RBFOX2 eCLIP in HepG2 at library positions in the transcriptome for 0 (n = 7041) or 1-3 (n = 711) NGCAU motifs. RBFOX2 peaks were compared to an IgG control to determine enrichments. f. RBFOX2 eCLIP in HepG2 at library positions in the transcriptome for 0 (n = 4610) or 1-3 (n = 3142) AU motifs. RBFOX2 peaks were compared to an IgG control to determine enrichments. 73 Figure S2. Different nsRBNS libraries emphasize different 5mer binding preferences for RBFOX2. a. R value distribution of nsRBNS sequences containing 1-2 copies of different 6mer classes UGNNUG (n = 7725), CGNNUG (n = 1751), AGNNUG (n = 6260), GGNNUG (n = 4935). b-c. Comparison of random (b) and intronic natural sequence (c) RBNS with 3' UTR nsRBNS 5mer enrichments. Primary and secondary motifs are labelled in red and blue, respectively. Dotted lines show 2.5 standard deviations above the mean. d. Filter binding with radiolabeled oligonucleotides containing three copies of the indicated sequence brought to equilibrium with six concentrations of RBFOX2. Primary motifs in gold, secondary motifs in teal, controls in grey. Error bars indicate +/- SD for three replicates. 74 75 Figure S3. RBFOX2 iCLIP demonstrates broad agreement with nsRBNS. a. Some secondary motifs show sharp peaks near 0 in a metaplot centered at the motif in introns (black) and 3' UTRs (grey) in RBFOX2 iCLIP data(Jangi et al., 2014). 5' ends of iCLIP reads containing the motif of interest were aligned with position one of the pentamer at 0 and normalized to the minimum read count in an 80-nt window (50-nt window shown). Y-axis range was reduced for secondary motifs. See Methods for read counts. b. AU-rich nsRBNS motifs do not show characteristic read peaks near 0 in a metaplot centered at the motif in introns (black) and 3' UTRs (grey) in RBFOX2 iCLIP data(Jangi et al., 2014). iCLIP reads containing the motif of interest were aligned with position one of the pentamer at 0 and normalized to the minimum read count in an 80-nt window (50-nt window shown). Y-axis range was reduced for secondary motifs. See methods for read counts. c. Schematic showing the generation of a clip enrichment (CE) score from iCLIP data. After generation of a metaplot, the read count at the peak apex was divided by the read count at its lowest point to generate a CE score analogous to an enrichment. d. Correlation of iCLIP- and nsRBNS-enriched 5mers in 3' UTRs (n = 1024). CLIP enrichment (CE) scores were computed for iCLIP peaks. Secondary motifs indicated in teal, primary motifs indicated in gold. 76 Figure S4. Enrichment of 5mers in HiTS-CLIP. 5mer enrichment of top 200 5mers in two HiTS-CLIP datasets in both introns and 3' UTRs. 5mer enrichment was calculated by determining the frequencies of all 1,024 5mers in CLIP peaks in each region and dataset and subsequently normalizing to control peaks from that region. Peaks from (a) Mouse ventral spinal neuron 3' UTR HiTS-CLIP, (b) Mouse whole brain intronic HiTS-CLIP, and (c) Mouse whole brain 3' UTR HiTS-CLIP were analyzed. Gold indicates primary motifs, teal indicates secondary motifs. 77 Figure S5. Representative raw data from flow cytometry. Graphs were drawn with pseudocolor in FlowJo. a. Gating strategy to select for single, live, intact cells. Events were gated through three serial gates to obtain approximately 25000 events for downstream analysis. Total number of events in each graph, and the percentage of events within the gate in each graph are shown. (FSC: forward scatter; SSC: side scatter; A: area; H: height; W: width.) b,c. Compensated values of the three fluorophores used (dsRED, EGFP, Cerulean), in positive and control samples with (b) primary and (c) secondary motifs. 78 Figure S6. Secondary motifs promote inclusion in a splicing reporter in an RBFOX1- dependent manner at the protein level. a. Six secondary motifs approximate the exon inclusion of one primary motif in an Rbfox1-dependent manner at the protein level, replicate 2. RG6 plasmids containing one primary motif or six secondary motifs were co-transfected in HEK293T cells with fluorescently labelled Rbfox1 and monitored by flow cytometry for the inclusion isoform (GFP), exclusion isoform (dsRED), and Rbfox1 (Cerulean) expression at the single-cell level. Controls including a scrambled motif co-transfected with Rbfox1 (light grey) and scrambled and intact motifs without Rbfox1 (grey) are also shown. Bins detailed in Supplementary Table 5. b. The slope of linear fit of two flow cytometry replicates were null- subtracted and normalized to their permuted controls. Error bars represent standard error of the mean (SEM). 79 Figure S7. Secondary motifs become engaged at specific intervals of neuronal differentiation. Pearson correlation of secondary motif presence with exon inclusion at intervals of neuronal differentiation beginning with embryonic stem cells and progressing to mature 28- day glutamatergic neurons (ESC–NESC (n = 448), NESC–RG (n = 1478), RG–DS1 (n = 940), DS1–DS3 (n = 2189), DS3–MAT16 (n = 1600), MAT16–MAT21 (n = 378), MAT21–MAT28 (n = 373)). Size of point indicates correlation coefficient, intensity indicates p-value < 0.05. 80 Figure S8. Estimation of secondary motif-dependent Rbfox events across neuronal cell types. In a comparison of neuronal cell types with medium to highest Rbfox mRNA expression, exons likely to be regulated by Rbfox are significantly (P < .0084 Fisher’s exact test, ndown = 13; nup = 28) enriched in secondary motifs. Of 864 alternative exons with increased splicing, 11% are primary, 26.4% primary and secondary, and 3.2% are 4+ secondary motif-associated. Exons with one to three secondary motif instances are also significantly enriched (P < 0.0012, Fisher’s exact test, ndown = 263; nup = 354). 81 Figure S9. Affinity estimation of Rbfox secondary motifs. RBNS 7-mer enrichments (R- value) for 1.1 μΜ RBFOX2 (a) and 1.3 μΜ RBFOX3 (b) binding were first corrected for non- specific contributions (R’ see Methods) and then linearly correlated with known dissociation constants (Kd) for RBFOX1 binding(Auweter et al., 2006a; Stoltz, 2015b). Correlation coefficients between log(R’) and log(Kd) were r=-0.955, P-value=8.379 x 10-9 (a) and r=-0.915, P-value=6.7 x 10-7 (b). Scatter plots show estimated Kd as a function of the original, uncorrected R-value. Resulting 7-mer Kd estimates were highly correlated between RBFOX2 and RBFOX3 (c) with r=0.763, P-value ≈ 0. Data for all 7-mers are shown on a logarithmic scale. Primary motif containing 7-mers are highlighted in gold (GCAUG), yellow (GCACG), and teal (secondary motifs GCUUG, GAAUG, GUUUG, GUGUG, GUAUG, GCCUG). Grouping 7- mers by their 5-mer content allows to estimate average Kds for each 5-mer (see Methods). A histogram of these 5-mer dissociation constants is shown in (d), with primary and secondary motifs highlighted as in (c). Motifs GCUUG, GAAUG and GUUUG were considered strong motifs. 136 non-primary or secondary 5-mers with partial overlap to primary motifs GCAUG, GCACG were excluded. 82 83 Figure S10. A model for Rbfox secondary motifs. a. A high nuclear mRNA expression weighted histogram of potential intronic Rbfox binding sites (1,000,000 mRNAs/cell with average half-life time of 3 hours). Motif 5mers in gold (GCAUG), yellow (GCACG), and teal (GCUUG, GAAUG, GUUUG, GUAUG). b. A low nuclear mRNA expression weighted histogram of potential intronic Rbfox binding sites (10x lower mRNA copies/cell and a half-life time of 4 hours). c-d. Predicted average Rbfox occupancies on 5mer motifs as a function of the nuclear Rbfox concentration in low (c) and high (d) mRNA scenarios. The low mRNA scenario predicts that the fraction of Rbfox bound to secondary motifs surpasses primary motifs at Rbfox levels > 1 μΜ. This is lower than estimates from the high mRNA scenario in main Figure 6 (~14 μΜ). Non-specific binding depicted in grey. e. Filter binding with radiolabeled oligonucleotide containing three copies of a primary (GCAUG) or secondary (GCUUG, GAAUG, GUUUG) were incubated to equilibrium in the presence of unlabeled, single copy GCAUG oligonucleotide at six concentrations of RBFOX2. As protein concentration increased, so did the fraction bound of labeled RNA for both primary and secondary motifs. Error bars indicate +/- SD of three replicates. 84 85 Chapter 3. Future Directions The biological impact of Rbfox secondary motifs Rbfox secondary motifs (suboptimal motifs that enable a second “wave” of splicing) characterized in this study represent the first observation of spatiotemporal regulation through suboptimal motifs in RNA-binding proteins––establishing the complete biological significance of these motifs will helpful in understanding the regulatory potential of suboptimal motifs for other RBPs. In our study, dozens of exons regulated exclusively by secondary motifs were identified computationally in both in vitro glutamatergic neuron differentiation and across neuronal subtypes. These exons can be validated as bona fide Rbfox targets by performing Rbfox CLIP at important stages of the in vitro differentiation and identifying reproducible peaks in the introns flanking candidate exons. Functional validation of these targets is also possible, as exemplified in the case of the Cd47 tri-chromatic splicing reporter assay in our study. Because most splicing regulation occurs within 400-nt of the alternative exon (Jaganathan et al., 2019) and secondary motifs were evaluated in the first 250-nt of the downstream intron in our analyses, a medium-throughput validation of the putative Rbfox secondary motif targets identified in both differentiation and diversification is feasible. Candidate proximal intron sequences can be synthesized and directly incorporated into a splicing reporter and co-expressed with exogenous Rbfox; secondary motifs can be mutationally ablated or blocked with antisense oligonucleotides specific to the candidate proximal intron. Furthermore, it may be instructive to survey other available in vitro differentiation datasets of neuronal subtypes where Rbfox expression changes substantially to find secondary motif-regulated splicing events and similarly validate those exons (Jacko et al., 2018; Mazin et al., 2021). 86 Direct experimental validation of secondary motif-regulated exons will provide a clearer picture of the global impact of secondary motifs on neuronal differentiation, and validated targets can be evaluated for their developmental behavior, biological function, and conservation. Exons demonstrably regulated by secondary motifs in the splicing reporter may be evaluated for their individual contributions to differentiation by abrogating their regulation during a differentiation time course. Initially, total abrogation of Rbfox binding can be achieved by using antisense oligonucleotides targeted to secondary motif clusters in the relevant intron. For particularly strong candidates (with clear secondary motif clusters and prior evidence for being temporally important in neuronal development), CRISPR/Cas9 genome editing could be employed at their endogenous loci to add primary GCAUG motif(s) to the intron, inducing exon inclusion earlier in the differentiation process, to assess the impact of mis-timed expression of Rbfox secondary motif-regulated exons. Of particular interest may be exons regulated by the lowest-affinity secondary motif that was activated at the peak of Rbfox expression (GUUUG), which were also responsive to the small decrease in Rbfox expression after peak expression; these transiently- activated exons may have interesting biological roles. However, any single target exon studied may not have dramatic effects on differentiation or diversification, as RBPs promote general regulatory programs and their impact is often best evaluated globally. It will be challenging to disentangle the overall effect of Rbfox secondary motifs on neuronal development from the many essential Rbfox primary motif-regulated exons. Titrating the overall expression level of Rbfox at the peak of Rbfox expression should result in a greater reduction in binding to secondary motifs compared to primary motifs, preferentially dysregulating secondary motif-regulated exons. As such, one method to probe the overall contribution of secondary motif-regulated exon inclusion in late neuronal development would be 87 incremental knockdown of Rbfox at peak expression by titrating shRNA with the use of an inducible promoter; better control could possibly be achieved with the degradation tag system with titration of the dTAG small molecule. Secondary motif-regulated exons should be preferentially dysregulated upon knockdown of Rbfox, enabling observation of their global effect on late neuronal differentiation. Another approach would be alteration of the binding properties of Rbfox itself by weakening its affinity for secondary motifs. While it may be possible to engineer the Rbfox RRM to more strongly require an adenine at position three in its GCAUG primary motif, boosting its affinity for primary motif GCAUG (and secondary motif GAAUG), replacing the RRM with Pumilio- homology domain built to recognize GCAUG may better approximate a global “loss” of secondary motif recognition (Adamala et al. 2016). Because Rbfox mediates exon inclusion with its C-terminal domain (with no contribution from its RRM) and exon exclusion with its RRM (with no contribution from its C-terminal domain), this system would be best applied in contexts like glutamatergic neuron differentiation, where Rbfox predominantly activates exon inclusion (Sun et al., 2012; Ying et al., 2017). The engineered Rbfox’s specificity should be biochemically characterized with RBNS; following characterization, expression of this engineered protein under control of one or more endogenous Rbfox promoters during neuronal differentiation may show the global effect of a loss of secondary motif binding on this biological process. A robust set of validated Rbfox secondary motif-regulated exons also invites further bioinformatic analyses in comparison to exons regulated by Rbfox primary motifs. It would be interesting to evaluate the overall expression and tissue-specificity of these exons—are they neuronally-specific? What is their prevalence across neuronal subtypes, and are there similarities in the genes they regulate? Another open question is the extent to which these more-abundant 88 secondary motifs provide raw material for the rewiring of the splicing regulatory network over evolutionary timescales. Comparing the syntenic primary and secondary Rbfox motifs in proximal downstream introns in mouse and human genomes could reveal whether or not Rbfox motifs have tended to strengthen or weaken over evolutionary time, and the overall conservation of the upstream exon relative to its motifs may also shed light on the evolution of splicing. Suboptimal motifs in neuronal development Regulation of splicing through binding to suboptimal motifs might be a common way to transduce spatiotemporal information, from extracellular morphogens to post-transcription gene regulation in neuronal development; they are likely to be particularly relevant to the alternative splicing of the mammalian brain. The work presented herein suggests that in order to functionally bind secondary motifs, an RBP must have two key properties: robust affinity for a specific RNA motif and sufficient RBP expression in some cells or tissues such that secondary motifs would be substantially occupied, presumably coincident with saturation of primary-motif sites. Notably, several RBP families with defined specificity exhibit large changes in expression over neuronal development (Weyn-Vanhentenryck et al., 2018). Furthermore, alternative splicing is an inherently dynamic process; alternative exons are variably included across tissues, developmental time, and cellular state, and it follows that the nuanced regulatory control afforded by suboptimal motifs would be advantageous in complex biological contexts. It is no surprise, then, that alternative splicing is most widespread in the brain, contributing to the vast complexity of post-transcriptional gene regulation in its tissues (Barbosa-Morais et al., 2012; Merkin et al., 2012; Zhang et al., 2014) by remodeling the levels of and presumed activity of both coding and non-coding RNAs (Ellis et al., 2012; Mockenhaupt & Makeyev, 2015; Yang et 89 al., 2016). In addition to the variation in splicing cells of the adult brain, alternative splicing is also most pronounced in the brain during the course of neuronal development: a survey of alternative splicing changes across development in human and animal tissues found that the brain harbored the greatest number of developmentally-regulated exons and that the human brain was enriched in both conserved and novel splicing events (Mazin et al., 2021). Further, individual neurons are highly specialized cells that often rely on unique mechanisms of post-transcriptional gene regulation through splicing. The morphological complexity of neurons, which includes multi-polar structures such as axons and dendrites with sites of localized, non-somatic translation of mRNAs, requires precise regulation of the spatial transcriptome (Schieweck et al., 2021). Overall, there are myriad neuronal mechanisms for post- transcriptional control across time, space, and in response to stimuli. In neurons, intron retention not only abrogates gene expression, but can induce rapid gene expression: polyadenylated, intron-retained transcripts accumulate in the nucleus in the mouse neocortex, and activity- dependent splicing can then rapidly release transcripts for translation in response to a specific cellular signal (Mauger et al. 2016). Alternative splicing coupled with nonsense-mediated decay has been shown to mediate the spatiotemporal expression of multiple proteins in neurons, demonstrating the use of splicing for spatiotemporal regulation in the cell (Giorgi et al. 2007, Chen et al. 2008). Of course, precise and unique modes of post-transcriptional regulation often do not proceed in isolation, and are especially entangled in the brain, where mechanisms of regulation by ncRNAs as diverse as lncRNAs, circRNAs, and miRNAs can be intertwined (Kleaveland et al., 2018). The use of suboptimal motifs to tune splicing may be particularly robust and widespread in neuronal development, where RBP expression is dynamically regulated. An extensive 90 characterization of the temporal regulation of alternative splicing in mouse cortical neurons confirmed the variable expression of multiple tissue-specific RBP families—Rbfox, Nova, MBNL, and PTBP—in neuronal development, and analysis of thousands of significantly changing cassette exons showed that these four RBP families predominate in their contributions to splicing over developmental time (Weyn-Vanhentenryck et al., 2018). Thus, in the investigation of the functional scope of suboptimal motifs usage by RBPs, it would be instructive to begin by biochemically interrogating the affinity landscapes of these proteins to incorporate suboptimal motifs into similar integrative modeling approaches to better understand the contribution(s) of secondary motif–binding to the essential, alternative-splicing functions of these proteins in the course of neuronal development. 91 References Adamala, K. P., Martin-Alarcon, D. A., Boyden, E. S. (2016). Programmable RNA-binding protein composed of repeats of a single modular unit. Proceedings of the National Academy of Sciences of the United States of America, 113(19). https://doi.org/10.1073/pnas.1519368113 Anantharaman, V., Koonin, E. v, & Aravind, L. (2002). Comparative genomics and evolution of proteins involved in RNA metabolism. Nucleic Acids Research, 30(7). https://doi.org/10.1093/nar/30.7.1427 Augenlicht, L. H., & Lipkin, M. (1976). Appearance of rapidly labeled, high molecular weight RNA in nuclear ribonucleoprotein. Release from chromatin and association with protein. The Journal of Biological Chemistry, 251(9). Auweter, S. D., Fasan, R., Reymond, L., Underwood, J. G., Black, D. L., Pitsch, S., & Allain, F. H. T. (2006). Molecular basis of RNA recognition by the human alternative splicing factor Fox-1. EMBO Journal. https://doi.org/10.1038/sj.emboj.7600918 Baltz, A. G., Munschauer, M., Schwanhäusser, B., Vasile, A., Murakawa, Y., Schueler, M., Youngs, N., Penfold-Brown, D., Drew, K., Milek, M., Wyler, E., Bonneau, R., Selbach, M., Dieterich, C., & Landthaler, M. (2012). The mRNA-bound proteome and its global occupancy profile on protein-coding transcripts. Molecular Cell, 46(5). https://doi.org/10.1016/j.molcel.2012.05.021 Barash, Y., Calarco, J. A., Gao, W., Pan, Q., Wang, X., Shai, O., Blencowe, B. J., & Frey, B. J. (2010). Deciphering the splicing code. Nature, 465(7294), 53–59. https://doi.org/10.1038/nature09000 Barbosa-Morais, N. L., Irimia, M., Pan, Q., Xiong, H. Y., Gueroussov, S., Lee, L. J., Slobodeniuc, V., Kutter, C., Watt, S., Çolak, R., Kim, T., Misquitta-Ali, C. M., Wilson, M. D., Kim, P. M., Odom, D. T., Frey, B. J., & Blencowe, B. J. (2012). The Evolutionary Landscape of Alternative Splicing in Vertebrate Species. Science, 338(6114), 1587–1593. https://doi.org/10.1126/science.1230612 Barnby, G., Abbott, A., Sykes, N., Morris, A., Weeks, D. E., Mott, R., Lamb, J., Bailey, A. J., & Monaco, A. P. (2005). Candidate-Gene Screening and Association Analysis at the Autism- Susceptibility Locus on Chromosome 16p: Evidence of Association at GRIN2A and ABAT. The American Journal of Human Genetics. https://doi.org/10.1086/430454 Beckmann, B. M., Horos, R., Fischer, B., Castello, A., Eichelbaum, K., Alleaume, A.-M., Schwarzl, T., Curk, T., Foehr, S., Huber, W., Krijgsveld, J., & Hentze, M. W. (2015). The RNA-binding proteomes from yeast to man harbour conserved enigmRBPs. Nature Communications, 6. https://doi.org/10.1038/ncomms10127 Benavides-Piccione, R., Regalado-Reyes, M., Fernaud-Espinosa, I., Kastanauskaite, A., Tapia- González, S., León-Espinosa, G., Rojo, C., Insausti, R., Segev, I., & DeFelipe, J. (2019). Differential Structure of Hippocampal CA1 Pyramidal Neurons in the Human and Mouse. Cerebral Cortex. https://doi.org/10.1093/cercor/bhz122 Bhalla, K., Phillips, H. A., Crawford, J., McKenzie, O. L. D., Mulley, J. C., Eyre, H., Gardner, A. E., Kremmidiotis, G., & Callen, D. F. (2004). The de novo chromosome 16 translocations of two patients with abnormal phenotypes (mental retardation and epilepsy) disrupt the A2BP1 gene. Journal of Human Genetics. https://doi.org/10.1007/s10038-004- 0145-4 92 Biswas, J., Patel, V. L., Bhaskar, V., Chao, J. A., Singer, R. H., & Eliscovich, C. (2019). The structural basis for RNA selectivity by the IMP family of RNA-binding proteins. Nature Communications, 10(1). https://doi.org/10.1038/s41467-019-12193-7 Bray, N. L., Pimentel, H., Melsted, P., & Pachter, L. (2016). Near-optimal probabilistic RNA-seq quantification. Nature Biotechnology. https://doi.org/10.1038/nbt.3519 Brudno, M. (2001). Computational analysis of candidate intron regulatory elements for tissue- specific alternative pre-mRNA splicing. Nucleic Acids Research. https://doi.org/10.1093/nar/29.11.2338 Buckanovich, R. J., & Darnell, R. B. (1997). The neuronal RNA binding protein Nova-1 recognizes specific RNA targets in vitro and in vivo. Molecular and Cellular Biology, 17(6). https://doi.org/10.1128/MCB.17.6.3194 Buckanovich, R. J., Posner, J. B., & Darnell, R. B. (1993). Nova, the paraneoplastic Ri antigen, is homologous to an RNA-binding protein and is specifically expressed in the developing motor system. Neuron, 11(4). https://doi.org/10.1016/0896-6273(93)90077-5 Buckanovich, R. J., Yang, Y. Y., & Darnell, R. B. (1996). The onconeural antigen Nova-1 is a neuron-specific RNA-binding protein, the activity of which is inhibited by paraneoplastic antibodies. The Journal of Neuroscience: The Official Journal of the Society for Neuroscience, 16(3). Busch, A., & Hertel, K. J. (2012). Evolution of SR protein and hnRNP splicing regulatory factors. Wiley Interdisciplinary Reviews. RNA, 3(1). https://doi.org/10.1002/wrna.100 Bush, S. J., Chen, L., Tovar-Corona, J. M., & Urrutia, A. O. (2017). Alternative splicing and the evolution of phenotypic novelty. Philosophical Transactions of the Royal Society B: Biological Sciences, 372(1713), 20150474. https://doi.org/10.1098/rstb.2015.0474 Butter, F., Scheibe, M., Mörl, M., & Mann, M. (2009). Unbiased RNA-protein interaction screen by quantitative proteomics. Proceedings of the National Academy of Sciences of the United States of America, 106(26). https://doi.org/10.1073/pnas.0812099106 Castello, A., Fischer, B., Eichelbaum, K., Horos, R., Beckmann, B. M., Strein, C., Davey, N. E., Humphreys, D. T., Preiss, T., Steinmetz, L. M., Krijgsveld, J., & Hentze, M. W. (2012). Insights into RNA biology from an atlas of mammalian mRNA-binding proteins. Cell, 149(6). https://doi.org/10.1016/j.cell.2012.04.031 Christofk, H. R., vander Heiden, M. G., Harris, M. H., Ramanathan, A., Gerszten, R. E., Wei, R., Fleming, M. D., Schreiber, S. L., & Cantley, L. C. (2008). The M2 splice isoform of pyruvate kinase is important for cancer metabolism and tumour growth. Nature, 452(7184). https://doi.org/10.1038/nature06734 Cléry, A., Blatter, M., & Allain, F. H.-T. (2008). RNA recognition motifs: boring? Not quite. Current Opinion in Structural Biology, 18(3). https://doi.org/10.1016/j.sbi.2008.04.002 Conboy, J. G. (2017). Developmental regulation of RNA processing by Rbfox proteins. In Wiley Interdisciplinary Reviews: RNA. https://doi.org/10.1002/wrna.1398 Cox, R. A. (1981). Biophysical Chemistry Part III: The Behavior of Biological Macromolecules: by C. R. Cantor and P. R. Schimmel Freeman; San Francisco, 1980 xxx + 524 pages. £22.40 (board); £11.90 (paper). FEBS Letters. https://doi.org/10.1016/0014-5793(81)80070- 1 Crocker, J., Abe, N., Rinaldi, L., McGregor, A. P., Frankel, N., Wang, S., Alsawadi, A., Valenti, P., Plaza, S., Payre, F., Mann, R. S., & Stern, D. L. (2015). Low Affinity Binding Site Clusters Confer Hox Specificity and Regulatory Robustness. Cell, 160(1–2), 191–203. https://doi.org/10.1016/j.cell.2014.11.041 93 Damianov, A., Ying, Y., Lin, C. H., Lee, J. A., Tran, D., Vashisht, A. A., Bahrami-Samani, E., Xing, Y., Martin, K. C., Wohlschlegel, J. A., & Black, D. L. (2016). Rbfox Proteins Regulate Splicing as Part of a Large Multiprotein Complex LASR. Cell. https://doi.org/10.1016/j.cell.2016.03.040 Dominguez, D., Freese, P., Alexis, M. S., Su, A., Hochman, M., Palden, T., Bazile, C., Lambert, N. J., Van Nostrand, E. L., Pratt, G. A., Yeo, G. W., Graveley, B. R., & Burge, C. B. (2018). Sequence, Structure, and Context Preferences of Human RNA Binding Proteins. Molecular Cell. https://doi.org/10.1016/j.molcel.2018.05.001 Driever, W., & Nüsslein-Volhard, C. (1988). A gradient of bicoid protein in Drosophila embryos. Cell, 54(1), 83–93. https://doi.org/10.1016/0092-8674(88)90182-1 Eden, E., Lipson, D., Yogev, S., & Yakhini, Z. (2007). Discovering motifs in ranked lists of DNA sequences. PLoS Computational Biology. https://doi.org/10.1371/journal.pcbi.0030039 Eden, E., Navon, R., Steinfeld, I., Lipson, D., & Yakhini, Z. (2009). GOrilla: A tool for discovery and visualization of enriched GO terms in ranked gene lists. BMC Bioinformatics. https://doi.org/10.1186/1471-2105-10-48 Ellis, J. D., Barrios-Rodiles, M., Colak, R., Irimia, M., Kim, T., Calarco, J. A., Wang, X., Pan, Q., O’Hanlon, D., Kim, P. M., Wrana, J. L., & Blencowe, B. J. (2012). Tissue-specific alternative splicing remodels protein-protein interaction networks. Molecular Cell, 46(6). https://doi.org/10.1016/j.molcel.2012.05.037 Fairbrother, W. G., & Chasin, L. A. (2000). Human Genomic Sequences That Inhibit Splicing. Molecular and Cellular Biology, 20(18), 6816–6825. https://doi.org/10.1128/MCB.20.18.6816-6825.2000 Farley, E. K., Olson, K. M., Zhang, W., Brandt, A. J., Rokhsar, D. S., & Levine, M. S. (2015). Suboptimization of developmental enhancers. Science. https://doi.org/10.1126/science.aac6948 Feingold, E. A., Good, P. J., Guyer, M. S., Kamholz, S., Liefer, L., Wetterstrand, K., Collins, F. S., Gingeras, T. R., Kampa, D., Sekinger, E. A., Cheng, J., Hirsch, H., Ghosh, S., Zhu, Z., Patel, S., Piccolboni, A., Yang, A., Tammana, H., Bekiranov, S., … Harvey, S. C. (2004). The ENCODE (ENCyclopedia of DNA Elements) Project. In Science. https://doi.org/10.1126/science.1105136 Fu, X.-D., & Ares, M. (2014). Context-dependent control of alternative splicing by RNA-binding proteins. Nature Reviews. Genetics, 15(10). https://doi.org/10.1038/nrg3778 Gallagher, T. L., Arribere, J. A., Geurts, P. A., Exner, C. R. T., McDonald, K. L., Dill, K. K., Marr, H. L., Adkar, S. S., Garnett, A. T., Amacher, S. L., & Conboy, J. G. (2011). Rbfox- regulated alternative splicing is critical for zebrafish cardiac and skeletal muscle functions. Developmental Biology. https://doi.org/10.1016/j.ydbio.2011.08.025 Gao, G., Xie, A., Huang, S.-C., Zhou, A., Zhang, J., Herman, A. M., Ghassemzadeh, S., Jeong, E.-M., Kasturirangan, S., Raicu, M., Sobieski, M. A., Bhat, G., Tatooles, A., Benz, E. J., Kamp, T. J., & Dudley, S. C. (2011). Role of RBM25/LUC7L3 in Abnormal Cardiac Sodium Channel Splicing Regulation in Human Heart Failure. Circulation, 124(10). https://doi.org/10.1161/CIRCULATIONAHA.111.044495 Gaudet, J., & Mango, S. E. (2002). Regulation of organogenesis by the Caenorhabditis elegans FoxA protein PHA-4. Science. https://doi.org/10.1126/science.1065175 Gehman, L. T., Meera, P., Stoilov, P., Shiue, L., O’Brien, J. E., Meisler, M. H., Ares, M., Otis, T. S., & Black, D. L. (2012). The splicing regulator Rbfox2 is required for both cerebellar 94 development and mature motor function. Genes and Development. https://doi.org/10.1101/gad.182477.111 Gehman, L. T., Stoilov, P., Maguire, J., Damianov, A., Lin, C. H., Shiue, L., Ares, M., Mody, I., & Black, D. L. (2011). The splicing regulator Rbfox1 (A2BP1) controls neuronal excitation in the mammalian brain. Nature Genetics. https://doi.org/10.1038/ng.841 Gerstberger, S., Hafner, M., & Tuschl, T. (2014). A census of human RNA-binding proteins. Nature Reviews. Genetics, 15(12). https://doi.org/10.1038/nrg3813 Gomes, E., & Shorter, J. (2019). The molecular language of membraneless organelles. Journal of Biological Chemistry. https://doi.org/10.1074/jbc.TM118.001192 Gomez, A. M., Traunmüller, L., & Scheiffele, P. (2021). Neurexins: molecular codes for shaping neuronal synapses. Nature Reviews Neuroscience, 22(3), 137–151. https://doi.org/10.1038/s41583-020-00415-7 Guenther, U.-P., Yandek, L. E., Niland, C. N., Campbell, F. E., Anderson, D., Anderson, V. E., Harris, M. E., & Jankowsky, E. (2013). Hidden specificity in an apparently nonspecific RNA-binding protein. Nature, 502(7471), 385–388. https://doi.org/10.1038/nature12543 Hamada, N., Ito, H., Nishijo, T., Iwamoto, I., Morishita, R., Tabata, H., Momiyama, T., & Nagata, K. I. (2016). Essential role of the nuclear isoform of RBFOX1, a candidate gene for autism spectrum disorders, in the brain development. Scientific Reports. https://doi.org/10.1038/srep30805 Harrow, J., Frankish, A., Gonzalez, J. M., Tapanari, E., Diekhans, M., Kokocinski, F., Aken, B. L., Barrell, D., Zadissa, A., Searle, S., Barnes, I., Bignell, A., Boychenko, V., Hunt, T., Kay, M., Mukherjee, G., Rajan, J., Despacio-Reyes, G., Saunders, G., … Hubbard, T. J. (2012). GENCODE: The reference human genome annotation for the ENCODE project. Genome Research. https://doi.org/10.1101/gr.135350.111 Helder, S., Blythe, A. J., Bond, C. S., & Mackay, J. P. (2016). Determinants of affinity and specificity in RNA-binding proteins. In Current Opinion in Structural Biology. https://doi.org/10.1016/j.sbi.2016.05.005 Hentze, M. W., Castello, A., Schwarzl, T., & Preiss, T. (2018). A brave new world of RNA- binding proteins. Nature Reviews. Molecular Cell Biology, 19(5). https://doi.org/10.1038/nrm.2017.130 Hiller, M., Zhang, Z., Backofen, R., & Stamm, S. (2007). Pre-mRNA secondary structures influence exon recognition. PLoS Genetics, 3(11). https://doi.org/10.1371/journal.pgen.0030204 Hodgkin, J., Zellan, J. D., & Albertson, D. G. (1994). Identification of a candidate primary sex determination locus, fox-1, on the X chromosome of Caenorhabditis elegans. Development. Hogan, D. J., Riordan, D. P., Gerber, A. P., Herschlag, D., & Brown, P. O. (2008). Diverse RNA-Binding Proteins Interact with Functionally Related Sets of RNAs, Suggesting an Extensive Regulatory System. PLoS Biology, 6(10), e255. https://doi.org/10.1371/journal.pbio.0060255 Hudson, B. P., Martinez-Yamout, M. A., Dyson, H. J., & Wright, P. E. (2004). Recognition of the mRNA AU-rich element by the zinc finger domain of TIS11d. Nature Structural & Molecular Biology, 11(3). https://doi.org/10.1038/nsmb738 Hunter, J. D. (2007). Matplotlib: A 2D graphics environment. Computing in Science and Engineering. https://doi.org/10.1109/MCSE.2007.55 Jacko, M., Weyn-Vanhentenryck, S. M., Smerdon, J. W., Yan, R., Feng, H., Williams, D. J., Pai, J., Xu, K., Wichterle, H., & Zhang, C. (2018). Rbfox Splicing Factors Promote Neuronal 95 Maturation and Axon Initial Segment Assembly. Neuron. https://doi.org/10.1016/j.neuron.2018.01.020 Jaeger, S. A., Chan, E. T., Berger, M. F., Stottmann, R., Hughes, T. R., & Bulyk, M. L. (2010). Conservation and regulatory associations of a wide affinity range of mouse transcription factor binding sites. Genomics, 95(4), 185–195. https://doi.org/10.1016/j.ygeno.2010.01.002 Jaganathan, K., Kyriazopoulou Panagiotopoulou, S., McRae, J. F., Darbandi, S. F., Knowles, D., Li, Y. I., Kosmicki, J. A., Arbelaez, J., Cui, W., Schwartz, G. B., Chow, E. D., Kanterakis, E., Gao, H., Kia, A., Batzoglou, S., Sanders, S. J., & Farh, K. K.-H. (2019). Predicting Splicing from Primary Sequence with Deep Learning. Cell, 176(3). https://doi.org/10.1016/j.cell.2018.12.015 Jangi, M., Boutz, P. L., Paul, P., & Sharp, P. A. (2014). Rbfox2 controls autoregulation in RNA- binding protein networks. Genes and Development. https://doi.org/10.1101/gad.235770.113 Jankowsky, E., & Harris, M. E. (2015). Specificity and nonspecificity in RNA-protein interactions. Nature Reviews Molecular Cell Biology. https://doi.org/10.1038/nrm4032 Jarmoskaite, I., Denny, S. K., Vaidyanathan, P. P., Becker, W. R., Andreasson, J. O. L., Layton, C. J., Kappel, K., Shivashankar, V., Sreenivasan, R., Das, R., Greenleaf, W. J., & Herschlag, D. (2019). A Quantitative and Predictive Model for RNA Binding by Human Pumilio Proteins. Molecular Cell, 74(5), 966-981.e18. https://doi.org/10.1016/j.molcel.2019.04.012 Jens, M., & Rajewsky, N. (2015). Competition between target sites of regulators shapes post- transcriptional gene regulation. Nature Reviews Genetics. https://doi.org/10.1038/nrg3853 Jensen, K. B., Musunuru, K., Lewis, H. A., Burley, S. K., & Darnell, R. B. (2000). The tetranucleotide UCAY directs the specific recognition of RNA by the Nova K-homology 3 domain. Proceedings of the National Academy of Sciences of the United States of America, 97(11). https://doi.org/10.1073/pnas.090553997 Jiang, J., & Levine, M. (1993). Binding affinities and cooperative interactions with bHLH activators delimit threshold responses to the dorsal gradient morphogen. Cell, 72(5), 741– 752. https://doi.org/10.1016/0092-8674(93)90402-C Jin, Y., Suzuki, H., Maegawa, S., Endo, H., Sugano, S., Hashimoto, K., Yasuda, K., & Inoue, K. (2003). A vertebrate RNA-binding protein Fox-1 regulates tissue-specific splicing via the pentanucleotide GCAUG. EMBO Journal. https://doi.org/10.1093/emboj/cdg089 Jones, E., Oliphant, T., Peterson, P., & Others. (2001). SciPy.org. SciPy: Open Source Scientific Tools for Python2. Jourdain, A. A., Begg, B. E., Mick, E., Shah, H., Calvo, S. E., Skinner, O. S., Sharma, R., Blue, S. M., Yeo, G. W., Burge, C. B., & Mootha, V. K. (2021). Loss of LUC7L2 and U1 snRNP subunits shifts energy metabolism from glycolysis to OXPHOS. Molecular Cell, 81(9). https://doi.org/10.1016/j.molcel.2021.02.033 Kim, E., Magen, A., & Ast, G. (2007). Different levels of alternative splicing among eukaryotes. Nucleic Acids Research, 35(1), 125–131. https://doi.org/10.1093/nar/gkl924 Kim, K. K., Adelstein, R. S., & Kawamoto, S. (2009). Identification of neuronal nuclei (NeuN) as Fox-3, a new member of the Fox-1 gene family of splicing factors. Journal of Biological Chemistry. https://doi.org/10.1074/jbc.M109.052969 Kishore, S., Jaskiewicz, L., Burger, L., Hausser, J., Khorshid, M., & Zavolan, M. (2011). A quantitative analysis of CLIP methods for identifying binding sites of RNA-binding proteins. Nature Methods. https://doi.org/10.1038/nmeth.1608 96 Kleaveland, B., Shi, C. Y., Stefano, J., & Bartel, D. P. (2018). A Network of Noncoding Regulatory RNAs Acts in the Mammalian Brain. Cell, 174(2), 350-362.e17. https://doi.org/10.1016/j.cell.2018.05.022 König, J., Zarnack, K., Rot, G., Curk, T., Kayikci, M., Zupan, B., Turner, D. J., Luscombe, N. M., & Ule, J. (2010). iCLIP reveals the function of hnRNP particles in splicing at individual nucleotide resolution. Nature Structural & Molecular Biology, 17(7). https://doi.org/10.1038/nsmb.1838 Kosik, K. S. (2016). Life at Low Copy Number: How Dendrites Manage with So Few mRNAs. In Neuron. https://doi.org/10.1016/j.neuron.2016.11.002 Kuroyanagi, H. (2009). Fox-1 family of RNA-binding proteins. In Cellular and Molecular Life Sciences. https://doi.org/10.1007/s00018-009-0120-5 Kuroyanagi, H., Ohno, G., Mitani, S., & Hagiwara, M. (2007). The Fox-1 Family and SUP-12 Coordinately Regulate Tissue-Specific Alternative Splicing In Vivo. Molecular and Cellular Biology. https://doi.org/10.1128/mcb.01508-07 Kuwasako, K., Takahashi, M., Unzai, S., Tsuda, K., Yoshikawa, S., He, F., Kobayashi, N., Güntert, P., Shirouzu, M., Ito, T., Tanaka, A., Yokoyama, S., Hagiwara, M., Kuroyanagi, H., & Muto, Y. (2014). RBFOX and SUP-12 sandwich a G base to cooperatively regulate tissue-specific splicing. Nature Structural and Molecular Biology. https://doi.org/10.1038/nsmb.2870 Lai, Z. C., Maxson, R., & Childs, G. (1988). Both basal and ontogenic promoter elements affect the timing and level of expression of a sea urchin H1 gene during early embryogenesis. Genes & Development, 2(2), 173–183. https://doi.org/10.1101/gad.2.2.173 Lambert, N., Robertson, A., Jangi, M., McGeary, S., Sharp, P. A., & Burge, C. B. (2014). RNA Bind-n-Seq: Quantitative Assessment of the Sequence and Structural Binding Specificity of RNA Binding Proteins. Molecular Cell. https://doi.org/10.1016/j.molcel.2014.04.016 Lee, E. H. Y., Hsieh, Y. P., Yang, C. L., Tsai, K. J., & Liu, C. H. (2000). Induction of integrin- associated protein (IAP) mRNA expression during memory consolidation in rat hippocampus. European Journal of Neuroscience. https://doi.org/10.1046/j.1460- 9568.2000.00985.x Lee, J. A., Damianov, A., Lin, C. H., Fontes, M., Parikshak, N. N., Anderson, E. S., Geschwind, D. H., Black, D. L., & Martin, K. C. (2016). Cytoplasmic Rbfox1 Regulates the Expression of Synaptic and Autism-Related Genes. Neuron. https://doi.org/10.1016/j.neuron.2015.11.025 Lewis, H. A., Chen, H., Edo, C., Buckanovich, R. J., Yang, Y. Y., Musunuru, K., Zhong, R., Darnell, R. B., & Burley, S. K. (1999). Crystal structures of Nova-1 and Nova-2 K- homology RNA-binding domains. Structure (London, England : 1993), 7(2). https://doi.org/10.1016/S0969-2126(99)80025-2 Li, J. J., Bickel, P. J., & Biggin, M. D. (2014). System wide analyses have underestimated protein abundances and the importance of transcription in mammals. PeerJ. https://doi.org/10.7717/peerj.270 Lim, L. P., & Burge, C. B. (2001). A computational analysis of sequence features involved in recognition of short introns. Proceedings of the National Academy of Sciences, 98(20). https://doi.org/10.1073/pnas.201407298 Lovci, M. T., Ghanem, D., Marr, H., Arnold, J., Gee, S., Parra, M., Liang, T. Y., Stark, T. J., Gehman, L. T., Hoon, S., Massirer, K. B., Pratt, G. A., Black, D. L., Gray, J. W., Conboy, J. G., & Yeo, G. W. (2013). Rbfox proteins regulate alternative mRNA splicing through 97 evolutionarily conserved RNA bridges. Nature Structural and Molecular Biology. https://doi.org/10.1038/nsmb.2699 Lu, G., Dolgner, S. J., & Hall, T. M. T. (2009). Understanding and engineering RNA sequence specificity of PUF proteins. Current Opinion in Structural Biology, 19(1). https://doi.org/10.1016/j.sbi.2008.12.009 Lunde, B. M., Moore, C., & Varani, G. (2007). RNA-binding proteins: modular design for efficient function. Nature Reviews. Molecular Cell Biology, 8(6). https://doi.org/10.1038/nrm2178 Martin, C. L., Duvall, J. A., Ilkin, Y., Simon, J. S., Arreaza, M. G., Wilkes, K., Alvarez- Retuerto, A., Whichello, A., Powell, C. M., Rao, K., Cook, E., & Geschwind, D. H. (2007). Cytogenetic and molecular characterization of A2BP1/FOX1 as a candidate gene for autism. American Journal of Medical Genetics, Part B: Neuropsychiatric Genetics. https://doi.org/10.1002/ajmg.b.30530 Mayrand, S. H., Pedersen, N., & Pederson, T. (1986). Identification of proteins that bind tightly to pre-mRNA during in vitro splicing. Proceedings of the National Academy of Sciences of the United States of America, 83(11). https://doi.org/10.1073/pnas.83.11.3718 Mazin, P. v., Khaitovich, P., Cardoso-Moreira, M., & Kaessmann, H. (2021). Alternative splicing during mammalian organ development. Nature Genetics, 53(6), 925–934. https://doi.org/10.1038/s41588-021-00851-w McGeary, S. E., Lin, K. S., Shi, C. Y., Pham, T. M., Bisaria, N., Kelley, G. M., & Bartel, D. P. (2019). The biochemical basis of microRNA targeting efficacy. Science, 366(6472). https://doi.org/10.1126/science.aav1741 McNutt, P. M., Hubbard, K. S., Gut, I. M., & Lyman, M. E. (2013). Longitudinal RNA sequencing of the deep transcriptome during neurogenesis of cortical glutamatergic neurons from murine ESCs. F1000Research. https://doi.org/10.12688/f1000research.2-35.v1 Merkin, J., Russell, C., Chen, P., & Burge, C. B. (2012). Evolutionary Dynamics of Gene and Isoform Regulation in Mammalian Tissues. Science, 338(6114), 1593–1599. https://doi.org/10.1126/science.1228186 Minovitsky, S., Gee, S. L., Schokrpur, S., Dubchak, I., & Conboy, J. G. (2005). The splicing regulatory element, UGCAUG, is phylogenetically and spatially conserved in introns that flank tissue-specific alternative exons. Nucleic Acids Research. https://doi.org/10.1093/nar/gki210 Mockenhaupt, S., & Makeyev, E. v. (2015). Non-coding functions of alternative pre-mRNA splicing in development. Seminars in Cell & Developmental Biology, 47–48, 32–39. https://doi.org/10.1016/j.semcdb.2015.10.018 Mordue, K. E., Hawley, B. R., Satchwell, T. J., & Toye, A. M. (2017). CD47 surface stability is sensitive to actin disruption prior to inclusion within the band 3 macrocomplex. Scientific Reports. https://doi.org/10.1038/s41598-017-02356-1 Moulé, Y., & Chauveau, J. (1968). 40 S ribonucleoprotein particles in rat liver nuclei. Journal of Molecular Biology, 33(2). https://doi.org/10.1016/0022-2836(68)90203-9 Mukherjee, N., Wessels, H.-H., Lebedeva, S., Sajek, M., Ghanbari, M., Garzia, A., Munteanu, A., Yusuf, D., Farazi, T., Hoell, J. I., Akat, K. M., Akalin, A., Tuschl, T., & Ohler, U. (2019). Deciphering human ribonucleoprotein regulatory networks. Nucleic Acids Research, 47(2), 570–581. https://doi.org/10.1093/nar/gky1185 Murata, T., Ohnishi, H., Okazawa, H., Murata, Y., Kusakari, S., Hayashi, Y., Miyashita, M., Itoh, H., Oldenborg, P. A., Furuya, N., & Matozaki, T. (2006). CD47 promotes neuronal 98 development through Src- and FRG/Vav2-mediated activation of Rac and Cdc42. Journal of Neuroscience. https://doi.org/10.1523/JNEUROSCI.3981-06.2006 Nicastro, G., Candel, A. M., Uhl, M., Oregioni, A., Hollingworth, D., Backofen, R., Martin, S. R., & Ramos, A. (2017). Mechanism of β-actin mRNA Recognition by ZBP1. Cell Reports, 18(5). https://doi.org/10.1016/j.celrep.2016.12.091 Nicastro, G., Taylor, I. A., & Ramos, A. (2015). KH-RNA interactions: back in the groove. Current Opinion in Structural Biology, 30. https://doi.org/10.1016/j.sbi.2015.01.002 Nilsen, T. W., & Graveley, B. R. (2010). Expansion of the eukaryotic proteome by alternative splicing. Nature, 463(7280), 457–463. https://doi.org/10.1038/nature08909 Nutiu, R., Friedman, R.C., Luo, S., Khrebtukova, I., Silva, D., Li, R., Zhang, L., Schroth, G. P., Burge, C. B. Direct visualization of DNA affinity landscapes using a high-throughput sequencing instrument. Nature Biotechnology, 29(7). https://dx.doi.org/10.1038%2Fnbt.1882 Orengo, J. P., Bundman, D., & Cooper, T. A. (2006). A bichromatic fluorescent reporter for cell- based screens of alternative splicing. Nucleic Acids Research. https://doi.org/10.1093/nar/gkl967 Ottoz, D. S. M., & Berchowitz, L. E. (2020). The role of disorder in RNA binding affinity and specificity. Open Biology, 10(12). https://doi.org/10.1098/rsob.200328 Pandit, S., Zhou, Y., Shiue, L., Coutinho-Mansfield, G., Li, H., Qiu, J., Huang, J., Yeo, G. W., Ares, M., & Fu, X.-D. (2013). Genome-wide analysis reveals SR protein cooperation and competition in regulated splicing. Molecular Cell, 50(2). https://doi.org/10.1016/j.molcel.2013.03.001 Papasaikas, P., Tejedor, J. R., Vigevani, L., & Valcárcel, J. (2015). Functional Splicing Network Reveals Extensive Regulatory Potential of the Core Spliceosomal Machinery. Molecular Cell, 57(1), 7–22. https://doi.org/10.1016/j.molcel.2014.10.030 Park, S., Myszka, D. G., Yu, M., Littler, S. J., & Laird-Offringa, I. A. (2000). HuD RNA recognition motifs play distinct roles in the formation of a stable complex with AU-rich RNA. Molecular and Cellular Biology, 20(13). https://doi.org/10.1128/MCB.20.13.4765- 4772.2000 Park, S., Phukan, P. D., Zeeb, M., Martinez-Yamout, M. A., Dyson, H. J., & Wright, P. E. (2017). Structural Basis for Interaction of the Tandem Zinc Finger Domains of Human Muscleblind with Cognate RNA from Human Cardiac Troponin T. Biochemistry, 56(32). https://doi.org/10.1021/acs.biochem.7b00484 Perez-Perri, J. I., Rogell, B., Schwarzl, T., Stein, F., Zhou, Y., Rettel, M., Brosig, A., & Hentze, M. W. (2018). Discovery of RNA-binding proteins and characterization of their dynamic responses by enhanced RNA interactome capture. Nature Communications, 9(1). https://doi.org/10.1038/s41467-018-06557-8 Pollard, K. S., Hubisz, M. J., Rosenbloom, K. R., & Siepel, A. (2010). Detection of nonneutral substitution rates on mammalian phylogenies. Genome Research. https://doi.org/10.1101/gr.097857.109 Quattrone, A., & Dassi, E. (2019). The Architecture of the Human RNA-Binding Protein Regulatory Network. IScience, 21, 706–719. https://doi.org/10.1016/j.isci.2019.10.058 Ray, D., Kazan, H., Cook, K. B., Weirauch, M. T., Najafabadi, H. S., Li, X., Gueroussov, S., Albu, M., Zheng, H., Yang, A., Na, H., Irimia, M., Matzat, L. H., Dale, R. K., Smith, S. A., Yarosh, C. A., Kelly, S. M., Nabet, B., Mecenas, D., … Hughes, T. R. (2013). A 99 compendium of RNA-binding motifs for decoding gene regulation. Nature, 499(7457), 172–177. https://doi.org/10.1038/nature12311 Rengifo-Gonzalez, J. C., el Hage, K., Clément, M.-J., Steiner, E., Joshi, V., Craveur, P., Durand, D., Pastré, D., & Bouhss, A. (2021). The cooperative binding of TDP-43 to GU-rich RNA repeats antagonizes TDP-43 aggregation. ELife, 10. https://doi.org/10.7554/eLife.67605 Rowan, S., Siggers, T., Lachke, S. A., Yue, Y., Bulyk, M. L., & Maas, R. L. (2010). Precise temporal control of the eye regulatory gene Pax6 via enhancer-binding site affinity. Genes and Development. https://doi.org/10.1101/gad.1890410 Sanders, D. W., Kedersha, N., Lee, D. S. W., Strom, A. R., Drake, V., Riback, J. A., Bracha, D., Eeftens, J. M., Iwanicki, A., Wang, A., Wei, M. T., Whitney, G., Lyons, S. M., Anderson, P., Jacobs, W. M., Ivanov, P., & Brangwynne, C. P. (2020). Competing Protein-RNA Interaction Networks Control Multiphase Intracellular Organization. Cell. https://doi.org/10.1016/j.cell.2020.03.050 Scardigli, R., Bäumer, N., Gruss, P., Guillemot, F., & le Roux, I. (2003). Direct and concentration-dependent regulation of the proneural gene Neurogenin2 by Pax6. Development, 130(14), 3269–3281. https://doi.org/10.1242/dev.00539 Scheiba, R. M., de Opakua, A. I., Díaz-Quintana, A., Cruz-Gallardo, I., Martínez-Cruz, L. A., Martínez-Chantar, M. L., Blanco, F. J., & Díaz-Moreno, I. (2014). The C-terminal RNA binding motif of HuR is a multi-functional domain leading to HuR oligomerization and binding to U-rich RNA targets. RNA Biology, 11(10). https://doi.org/10.1080/15476286.2014.996069 Scheiffele, P., Fan, J., Choih, J., Fetter, R., & Serafini, T. (2000). Neuroligin Expressed in Nonneuronal Cells Triggers Presynaptic Development in Contacting Axons. Cell, 101(6), 657–669. https://doi.org/10.1016/S0092-8674(00)80877-6 Schieweck, R., Ninkovic, J., & Kiebler, M. A. (2021). RNA-binding proteins balance brain function in health and disease. Physiological Reviews, 101(3), 1309–1370. https://doi.org/10.1152/physrev.00047.2019 Schwanhüusser, B., Busse, D., Li, N., Dittmar, G., Schuchhardt, J., Wolf, J., Chen, W., & Selbach, M. (2011). Global quantification of mammalian gene expression control. Nature. https://doi.org/10.1038/nature10098 Sebat, J., Lakshmi, B., Malhotra, D., Troge, J., Lese-Martin, C., Walsh, T., Yamrom, B., Yoon, S., Krasnitz, A., Kendall, J., Leotta, A., Pai, D., Zhang, R., Lee, Y. H., Hicks, J., Spence, S. J., Lee, A. T., Puura, K., Lehtimäki, T., … Wigler, M. (2007). Strong association of de novo copy number mutations with autism. Science. https://doi.org/10.1126/science.1138659 Segal, E., Raveh-Sadka, T., Schroeder, M., Unnerstall, U., & Gaul, U. (2008). Predicting expression patterns from regulatory sequence in Drosophila segmentation. Nature, 451(7178), 535–540. https://doi.org/10.1038/nature06496 Sellier, C., Cerro-Herreros, E., Blatter, M., Freyermuth, F., Gaucherot, A., Ruffenach, F., Sarkar, P., Puymirat, J., Udd, B., Day, J. W., Meola, G., Bassez, G., Fujimura, H., Takahashi, M. P., Schoser, B., Furling, D., Artero, R., Allain, F. H. T., Llamusi, B., & Charlet-Berguerand, N. (2018). RbFOX1/MBNL1 competition for CCUG RNA repeats binding contributes to myotonic dystrophy type 1/type 2 differences. Nature Communications. https://doi.org/10.1038/s41467-018-04370-x Sharma, D., Zagore, L. L., Brister, M. M., Ye, X., Crespo-Hernández, C. E., Licatalosi, D. D., & Jankowsky, E. (2021). The kinetic landscape of an RNA-binding protein in cells. Nature, 591(7848), 152–156. https://doi.org/10.1038/s41586-021-03222-x 100 Sharova, L. V., Sharov, A. A., Nedorezov, T., Piao, Y., Shaik, N., & Ko, M. S. H. (2009). Database for mRNA half-life of 19 977 genes obtained by DNA microarray analysis of pluripotent and differentiating mouse embryonic stem cells. DNA Research. https://doi.org/10.1093/dnares/dsn030 Shen, M., & Mattox, W. (2012). Activation and repression functions of an SR splicing regulator depend on exonic versus intronic-binding position. Nucleic Acids Research, 40(1). https://doi.org/10.1093/nar/gkr713 Shen, S., Park, J. W., Lu, Z. X., Lin, L., Henry, M. D., Wu, Y. N., Zhou, Q., & Xing, Y. (2014). rMATS: Robust and flexible detection of differential alternative splicing from replicate RNA-Seq data. Proceedings of the National Academy of Sciences of the United States of America. https://doi.org/10.1073/pnas.1419161111 Skipper, M., Milne, C. A., & Hodgkin, J. (1999). Genetic and molecular analysis of fox-1, a numerator element involved in Caenorhabditis elegans primary sex determination. Genetics. Sommerville, J. (1973). Ribonucleoprotein particles derived from the lampbrush chromosomes of newt oocytes. Journal of Molecular Biology, 78(3). https://doi.org/10.1016/0022- 2836(73)90470-1 Stitzinger, S. H., Sohrabi-Jahromi, S., & Söding, J. (2021). Cooperativity boosts affinity and specificity of proteins with multiple RNA-binding domains. BioRxiv. Stoltz, M. (2015). Interactions of the alternative splicing factor RBFOX with non-coding RNAs. ETH Zurich. Sun, S., Zhang, Z., Fregoso, O., & Krainer, A. R. (2012). Mechanisms of activation and repression by the alternative splicing factors RBFOX1/2. RNA. https://doi.org/10.1261/rna.030486.111 Sundararaman, B., Zhan, L., Blue, S. M., Stanton, R., Elkins, K., Olson, S., Wei, X., Van Nostrand, E. L., Pratt, G. A., Huelga, S. C., Smalec, B. M., Wang, X., Hong, E. L., Davidson, J. M., Lécuyer, E., Graveley, B. R., & Yeo, G. W. (2016). Resources for the Comprehensive Discovery of Functional RNA Elements. Molecular Cell, 61(6), 903–913. https://doi.org/10.1016/j.molcel.2016.02.012 Taliaferro, J. M., Lambert, N. J., Sudmant, P. H., Dominguez, D., Merkin, J. J., Alexis, M. S., Bazile, C. A., & Burge, C. B. (2016). RNA Sequence Context Effects Measured In Vitro Predict In Vivo Protein Binding and Regulation. Molecular Cell. https://doi.org/10.1016/j.molcel.2016.08.035 Tanaka Hall, T. (2005). Multiple modes of RNA recognition by zinc finger proteins. Current Opinion in Structural Biology, 15(3), 367–373. https://doi.org/10.1016/j.sbi.2005.04.004 Tanay, A. (2006). Extensive low-affinity transcriptional interactions in the yeast genome. Genome Research, 16(8), 962–972. https://doi.org/10.1101/gr.5113606 Teplova, M., Malinina, L., Darnell, J. C., Song, J., Lu, M., Abagyan, R., Musunuru, K., Teplov, A., Burley, S. K., Darnell, R. B., & Patel, D. J. (2011). Protein-RNA and protein-protein recognition by dual KH1/2 domains of the neuronal splicing factor Nova-1. Structure (London, England : 1993), 19(7). https://doi.org/10.1016/j.str.2011.05.002 Tsai, A., Muthusamy, A. K., Alves, M. R., Lavis, L. D., Singer, R. H., Stern, D. L., & Crocker, J. (2017). Nuclear microenvironments modulate transcription from low-affinity enhancers. ELife, 6. https://doi.org/10.7554/eLife.28975 Ule, J., & Blencowe, B. J. (2019). Alternative Splicing Regulatory Networks: Functions, Mechanisms, and Evolution. Molecular Cell, 76(2), 329–345. https://doi.org/10.1016/j.molcel.2019.09.017 101 Ule, J., Jensen, K. B., Ruggiu, M., Mele, A., Ule, A., & Darnell, R. B. (2003). CLIP identifies Nova-regulated RNA networks in the brain. Science (New York, N.Y.), 302(5648). https://doi.org/10.1126/science.1090095 Ule, J., Stefani, G., Mele, A., Ruggiu, M., Wang, X., Taneri, B., Gaasterland, T., Blencowe, B. J., & Darnell, R. B. (2006). An RNA map predicting Nova-dependent splicing regulation. Nature, 444(7119). https://doi.org/10.1038/nature05304 Ule, J., Ule, A., Spencer, J., Williams, A., Hu, J.-S., Cline, M., Wang, H., Clark, T., Fraser, C., Ruggiu, M., Zeeberg, B. R., Kane, D., Weinstein, J. N., Blume, J., & Darnell, R. B. (2005). Nova regulates brain-specific splicing to shape the synapse. Nature Genetics, 37(8). https://doi.org/10.1038/ng1610 Uren, P. J., Bahrami-Samani, E., Burns, S. C., Qiao, M., Karginov, V. F., Hodges, E., Hannon, G. J., Sanford, J. R., Penalva, L. O. F., & Smith, A. D. (2012). Site identification in high- throughput RNA-protein interaction data. Bioinformatics. https://doi.org/10.1093/bioinformatics/bts569 Van Nostrand, E. L., Freese, P., Pratt, G. A., Wang, X., Wei, X., Xiao, R., Blue, S. M., Chen, J.- Y., Cody, N. A. L., Dominguez, D., Olson, S., Sundararaman, B., Zhan, L., Bazile, C., Bouvrette, L. P. B., Bergalet, J., Duff, M. O., Garcia, K. E., Gelboin-Burkhart, C., … Yeo, G. W. (2020). A large-scale binding and functional map of human RNA-binding proteins. Nature, 583(7818), 711–719. https://doi.org/10.1038/s41586-020-2077-3 Van Nostrand, E. L., Pratt, G. A., Shishkin, A. A., Gelboin-Burkhart, C., Fang, M. Y., Sundararaman, B., Blue, S. M., Nguyen, T. B., Surka, C., Elkins, K., Stanton, R., Rigo, F., Guttman, M., & Yeo, G. W. (2016). Robust transcriptome-wide discovery of RNA-binding protein binding sites with enhanced CLIP (eCLIP). Nature Methods. https://doi.org/10.1038/nmeth.3810 Vuong, C. K., Wei, W., Lee, J. A., Lin, C. H., Damianov, A., de la Torre-Ubieta, L., Halabi, R., Otis, K. O., Martin, K. C., O’Dell, T. J., & Black, D. L. (2018). Rbfox1 Regulates Synaptic Transmission through the Inhibitory Neuron-Specific vSNARE Vamp1. Neuron. https://doi.org/10.1016/j.neuron.2018.03.008 Wagner, S. D., Struck, A. J., Gupta, R., Farnsworth, D. R., Mahady, A. E., Eichinger, K., Thornton, C. A., Wang, E. T., & Berglund, J. A. (2016). Dose-Dependent Regulation of Alternative Splicing by MBNL Proteins Reveals Biomarkers for Myotonic Dystrophy. PLoS Genetics. https://doi.org/10.1371/journal.pgen.1006316 Wang, E. T., Sandberg, R., Luo, S., Khrebtukova, I., Zhang, L., Mayr, C., Kingsmore, S. F., Schroth, G. P., & Burge, C. B. (2008). Alternative isoform regulation in human tissue transcriptomes. Nature, 456(7221), 470–476. https://doi.org/10.1038/nature07509 Wang, J., Malecka, A., Trøen, G., & Delabie, J. (2015). Comprehensive genome-wide transcription factor analysis reveals that a combination of high affinity and low affinity DNA binding is needed for human gene regulation. BMC Genomics. https://doi.org/10.1186/1471-2164-16-S7-S12 Wang, X., & Tanaka Hall, T. M. (2001). Structural basis for recognition of AU-rich element RNA by the HuD. Nature Structural Biology, 8(2). https://doi.org/10.1038/84131 Wang, Y., Ma, M., Xiao, X., & Wang, Z. (2012). Intronic splicing enhancers, cognate splicing factors and context-dependent regulation rules. Nature Structural & Molecular Biology, 19(10). https://doi.org/10.1038/nsmb.2377 102 Wang, Z., Hoffmann, H. M., & Grabowski, P. J. (1995). Intrinsic U2AF binding is modulated by exon enhancer signals in parallel with changes in splicing activity. RNA (New York, N.Y.), 1(1). Wang, Z., Rolish, M. E., Yeo, G., Tung, V., Mawson, M., & Burge, C. B. (2004). Systematic identification and analysis of exonic splicing silencers. Cell, 119(6). https://doi.org/10.1016/j.cell.2004.11.010 Wang, Z., Xiao, X., van Nostrand, E., & Burge, C. B. (2006). General and Specific Functions of Exonic Splicing Silencers in Splicing Control. Molecular Cell, 23(1), 61–70. https://doi.org/10.1016/j.molcel.2006.05.018 Weyn-Vanhentenryck, S. M., Feng, H., Ustianenko, D., Duffié, R., Yan, Q., Jacko, M., Martinez, J. C., Goodwin, M., Zhang, X., Hengst, U., Lomvardas, S., Swanson, M. S., & Zhang, C. (2018). Precise temporal regulation of alternative splicing during neural development. Nature Communications. https://doi.org/10.1038/s41467-018-04559-0 Weyn-Vanhentenryck, S. M., Mele, A., Yan, Q., Sun, S., Farny, N., Zhang, Z., Xue, C., Herre, M., Silver, P. A., Zhang, M. Q., Krainer, A. R., Darnell, R. B., & Zhang, C. (2014). HITS- CLIP and Integrative Modeling Define the Rbfox Splicing-Regulatory Network Linked to Brain Development and Autism. Cell Reports. https://doi.org/10.1016/j.celrep.2014.02.005 Wiśniewski, J. R., Hein, M. Y., Cox, J., & Mann, M. (2014). A “proteomic ruler” for protein copy number and concentration estimation without spike-in standards. Molecular and Cellular Proteomics. https://doi.org/10.1074/mcp.M113.037309 Wright, J. E., Gaidatzis, D., Senften, M., Farley, B. M., Westhof, E., Ryder, S. P., & Ciosk, R. (2011). A quantitative RNA code for mRNA target selection by the germline fate determinant GLD-1. The EMBO Journal, 30(3), 533–545. https://doi.org/10.1038/emboj.2010.334 Xiao, X., Wang, Z., Jang, M., Nutiu, R., Wang, E. T., & Burge, C. B. (2009). Splice site strength-dependent activity and genetic buffering by poly-G runs. Nature Structural and Molecular Biology. https://doi.org/10.1038/nsmb.1661 Xiong, H. Y., Alipanahi, B., Lee, L. J., Bretschneider, H., Merico, D., Yuen, R. K. C., Hua, Y., Gueroussov, S., Najafabadi, H. S., Hughes, T. R., Morris, Q., Barash, Y., Krainer, A. R., Jojic, N., Scherer, S. W., Blencowe, B. J., & Frey, B. J. (2015). RNA splicing. The human splicing code reveals new insights into the genetic determinants of disease. Science (New York, N.Y.), 347(6218). https://doi.org/10.1126/science.1254806 Yang, P., Wang, D., & Kang, L. (2021). Alternative splicing level related to intron size and organism complexity. BMC Genomics, 22(1), 853. https://doi.org/10.1186/s12864-021- 08172-2 Yang, X., Coulombe-Huntington, J., Kang, S., Sheynkman, G. M., Hao, T., Richardson, A., Sun, S., Yang, F., Shen, Y. A., Murray, R. R., Spirohn, K., Begg, B. E., Duran-Frigola, M., MacWilliams, A., Pevzner, S. J., Zhong, Q., Trigg, S. A., Tam, S., Ghamsari, L., … Vidal, M. (2016). Widespread Expansion of Protein Interaction Capabilities by Alternative Splicing. Cell, 164(4). https://doi.org/10.1016/j.cell.2016.01.029 Yang, Y. Y., Yin, G. L., & Darnell, R. B. (1998). The neuronal RNA-binding protein Nova-2 is implicated as the autoantigen targeted in POMA patients with dementia. Proceedings of the National Academy of Sciences of the United States of America, 95(22). https://doi.org/10.1073/pnas.95.22.13254 Ye, X., & Jankowsky, E. (2020). High throughput approaches to study RNA-protein interactions in vitro. Methods, 178, 3–10. https://doi.org/10.1016/j.ymeth.2019.09.006 103 Yeo, G., & Burge, C. B. (2004). Maximum Entropy Modeling of Short Sequence Motifs with Applications to RNA Splicing Signals. Journal of Computational Biology, 11(2–3), 377– 394. https://doi.org/10.1089/1066527041410418 Yeo, G. W., Coufal, N. G., Liang, T. Y., Peng, G. E., Fu, X. D., & Gage, F. H. (2009). An RNA code for the FOX2 splicing regulator revealed by mapping RNA-protein interactions in stem cells. Nature Structural and Molecular Biology. https://doi.org/10.1038/nsmb.1545 Yeo, G. W., van Nostrand, E. L., Nostrand, E. L. van, & Liang, T. Y. (2007). Discovery and analysis of evolutionarily conserved intronic splicing regulatory elements. PLoS Genetics, 3(5). https://doi.org/10.1371/journal.pgen.0030085 Ying, Y., Wang, X. J., Vuong, C. K., Lin, C. H., Damianov, A., & Black, D. L. (2017). Splicing Activation by Rbfox Requires Self-Aggregation through Its Tyrosine-Rich Domain. Cell. https://doi.org/10.1016/j.cell.2017.06.022 Zandvakili, A., Campbell, I., Gutzwiller, L. M., Weirauch, M. T., & Gebelein, B. (2018). Degenerate Pax2 and Senseless binding motifs improve detection of low-affinity sites required for enhancer specificity. PLOS Genetics, 14(4), e1007289. https://doi.org/10.1371/journal.pgen.1007289 Zhang, C., Frias, M. A., Mele, A., Ruggiu, M., Eom, T., Marney, C. B., Wang, H., Licatalosi, D. D., Fak, J. J., & Darnell, R. B. (2010). Integrative modeling defines the Nova splicing- regulatory network and its combinatorial controls. Science (New York, N.Y.), 329(5990). https://doi.org/10.1126/science.1191150 Zhang, C., Zhang, Z., Castle, J., Sun, S., Johnson, J., Krainer, A. R., & Zhang, M. Q. (2008). Defining the regulatory network of the tissue-specific splicing factors Fox-1 and Fox-2. Genes and Development. https://doi.org/10.1101/gad.1703108 Zhang, L., Vielle, A., Espinosa, S., & Zhao, R. (2019). RNAs in the spliceosome: Insight from cryoEM structures. WIREs RNA, 10(3). https://doi.org/10.1002/wrna.1523 Zhang, Y., Chen, K., Sloan, S. A., Bennett, M. L., Scholze, A. R., O’Keeffe, S., Phatnani, H. P., Guarnieri, P., Caneda, C., Ruderisch, N., Deng, S., Liddelow, S. A., Zhang, C., Daneman, R., Maniatis, T., Barres, B. A., & Wu, J. Q. (2014). An RNA-Sequencing Transcriptome and Splicing Database of Glia, Neurons, and Vascular Cells of the Cerebral Cortex. Journal of Neuroscience, 34(36), 11929–11947. https://doi.org/10.1523/JNEUROSCI.1860-14.2014 104 105 Article Exon-Mediated Activation of Transcription Starts Graphical Abstract Authors Ana Fiszbein, Keegan S. Krick, Bridget E. Begg, Christopher B. Burge Correspondence cburge@mit.edu RNAPII In Brief Fiszbein et al. show that splicing of internal exons activates transcription from nearby upstream promoters, especially those that are weak or even cryptic, suggesting that regulation of splicing can be used to alter the transcriptional output of RNAPII mammalian genes. Highlights d New promoters arise near evolutionarily new internal exons d Splicing of internal exons activates proximal upstream weak promoters d Splicing may recruit transcription machinery locally to influence promoter selection d These impacts of splicing on transcription are widespread Fiszbein et al., 2019, Cell 179, 1551–1565 December 12, 2019 ª 2019 Elsevier Inc. https://doi.org/10.1016/j.cell.2019.11.002 EMATS activation EMATS inhibition Article Exon-Mediated Activation of Transcription Starts Ana Fiszbein,1 Keegan S. Krick,1 Bridget E. Begg,1 and Christopher B. Burge1,2,* 1Department of Biology, Massachusetts Institute of Technology, Cambridge, MA 02138, USA 2Lead Contact *Correspondence: cburge@mit.edu https://doi.org/10.1016/j.cell.2019.11.002 SUMMARY can reduce levels of histone 3 lysine 4 trimethyl (H3K4me3), a chromatin mark associated with active transcription (Bieberstein The processing of RNA transcripts from mammalian et al., 2012). genes occurs in proximity to their transcription. Several components of the splicing machinery associate with Here, we describe a phenomenon affecting thou- RNA polymerase II (RNAPII) and other transcription machinery sands of genes that we call exon-mediated activation (Das et al., 2007; Emili et al., 2002; Kameoka et al., 2004; Morris of transcription starts (EMATS), in which the splicing and Greenleaf, 2000; Mortillaro et al., 1996; Neugebauer and of internal exons impacts promoter choice and the Roth, 1997; Vincent et al., 1996). The U1 and U2 small nuclear ribonucleoprotein particles (snRNPs) associate with general expression level of the gene. We observed that transcription factors (GTFs) TFIIH/GTF2H1 (Kwek et al., 2002), evolutionary gain of internal exons is associated TFIIF/GTF2F2 (Kameoka et al., 2004), and the carboxy-terminal with gain of new transcription start sites (TSSs) domain (CTD) of RNAPII (Emili et al., 2002; Morris and Greenleaf, nearby and increased gene expression. Inhibiting 2000). In addition to its role in splicing, U1 snRNP acts as a exon splicing reduced transcription from nearby pro- general repressor of proximal downstream premature cleavage moters, and creation of new spliced exons activated and polyadenylation (PCPA) sites (Gunderson et al., 1998; Kaida transcription from cryptic promoters. The strongest et al., 2010). The relative abundance of U1 snRNP binding effects occurred for weak promoters located prox- sites upstream in the antisense orientation from promoters imal and upstream of efficiently spliced exons. contributes to frequent termination of antisense transcripts at Together, our findings support a model in which PCPA sites, yielding short and unstable transcripts (Almada splicing recruits transcriptionmachinery locally to in- et al., 2013). Alternative transcription initiation and termination sites drive a fluence TSS choice and identify exon gain, loss, and substantial portion of transcript isoform differences between hu- regulatory change asmajor contributors to the evolu- man tissues (Reyes and Huber, 2018). Recent analyses of full- tion of alternative promoters and gene expression in length mRNAs suggest that transcription starts and splicing mammals. may be coordinated (Anvar et al., 2018). However, whether exon splicing commonly impacts transcription start site (TSS) INTRODUCTION location and activity remains unknown. Here, we describe a phe- nomenon we call exon-mediated activation of transcription RNA transcripts from mammalian genes are processed within starts (EMATS) in which the splicing of internal exons, especially seconds or minutes after their synthesis, providing opportunities those near gene 50 ends, alters gene expression by activating for functional connections between transcription and splicing nearby TSSs, contributing to expression regulation of thousands (Custódio and Carmo-Fonseca, 2016). Several links between of genes. splicing and transcription are known, and both transcription rate and chromatin structure can influence splicing outcomes RESULTS (Bentley, 2014; Kornblihtt et al., 2013; Schor et al., 2013). How- ever, more recent evidence suggests that splicing also feeds Increased Exon Splicing Is Associated with Increased back on transcription (Braunschweig et al., 2013). Adding an Gene Expression and Alternative TSS Usage intron to an intron-less gene often boosts gene expression in We used a comparative approach to explore potential connec- plants, animals, and fungi; although themechanisms are not fully tions between splicing and TSS usage, examining transcript pat- understood, impacts on transcription, nuclear export, mRNA terns in orthologous genes of mouse and rat that differed by the stability, and/or translation have been reported (Furger et al., presence or absence of an internal exon. Previously, we identi- 2002; Shaul, 2017). Splicing can impact the rate of transcription fied over 1,000 such exons that were not detected in RNA-seq elongation (Fong and Zhou, 2001), and in yeast, the presence of data from diverse tissues of other mammals including rat, ma- an intron can generate a transcriptional checkpoint that is asso- caque, and cow, and therefore likely arose recently in the murine ciated with pre-spliceosome formation (Chathoth et al., 2014). lineage. We also identified a similar number of exons that are Furthermore, recruitment of the spliceosome complex can stim- unique to the rat, as well as several hundred exons uniquely ulate transcription initiation by enhancing preinitiation complex lost in mouse or in rat (Figure 1A). Most of these evolutionarily assembly (Damgaard et al., 2008), and inhibition of splicing new exons are located in 50 untranslated regions (UTRs) and Cell 179, 1551–1565, December 12, 2019 ª 2019 Elsevier Inc. 1551 A Mouse-specific Rat-specific B new exon lost exon **** 25 mouse rs f ye a 0.4 o 65 rats llio n mi 20 macaque 0.2 190 cow 0.0 ψ < 0.05 < ψ ψ > chicken 0.05 < 0.95 0.95 n=6482 n=579 n=2865 C D **** genes without 10 all genes 1.0 new exon genes with mouse-specific new exons 0.5 genes in tissues w/ new exon 0.0 ψ < 0.05 5 −0.5 genes in p value (KS) <<< 0.01 tissues w/ new exon −1.0 ψ > 0.05 no. of TSS > 1 (RNA-seq) 0 10 20 30 40 50 60 no. of TSS (Start-seq) E all genes F G genes with mouse-specific new exons **** 2 **** 2 genes w/ mouse TSS 1 > rat TSS 1 genes w/ 0 0 mouse TSS = rat TSS −1 −1 genes w/ −2 mouse TSS −2 < rat TSS genes genes in genes in all genes with without tissues w/ tissues w/ genes rat-specific new exon new exon new exon lost exons 0 10 20 30 40 50 ψ < 0.05 ψ > 0.05 proportion of genes (%) Figure 1. Splicing of New Exons Is Associated with Increased Gene Expression and Gain of TSSs (A) Phylogenetic tree representing the main species used for dating evolutionarily new exons and approximate branch lengths in millions of years. The patterns of inclusion and exclusion used to infer mouse-specific new exons (n = 1,089) and rat-specific lost exons (n = 515) are shown. (B) Fold change in gene expression in genes with mouse-specific exons (assessed by fragments per kilobase of exon per million mapped reads [FPKM]) between mouse and rat in 9 organs, binned by c value of the new exon in each tissue. Number of gene-tissue pairs in each category is indicated. Mean ± SEM of displayed distributions is shown. (C) Fold change in gene expression between corresponding tissues of mouse and rat in genes with multiple TSSs in mouse (no. of TSS > 1) for mouse control genes with no new exons (white), genes with mouse-specific new exons in tissues where inclusion of the new exon is not detected (PSI < 0.05 [gray]), and genes with new mouse-specific exons in tissues were the exon is included (PSI > 0.05 [pink]). (D) Distribution of the number of TSSs per gene using Start-seq data frommurine macrophages for all genes expressed in mouse and genes with mouse-specific new exons. Start-seq peaks located within 50 bp from each other were merged. Genes with mouse-specific new exons have increased numbers of TSSs (p < 2.2e–16 by Kolmogorov-Smirnov test). (legend continued on next page) 1552 Cell 179, 1551–1565, December 12, 2019 (FPKM mouse / FPKM rat), log2 *** *** (no. TSS mouse / no. TSS rat), log proportion of genes (%)2 (FPKM mouse / FPKM rat), log2 (no. TSS rat / no. TSS mouse), log2 NS are spliced in an alternative and tissue-specific fashion (Merkin (Figures 2A and S2A), confirming a tight connection between et al., 2015). Comparing closely related species, we have evolution of promoters, internal exons, and gene expression observed that genes with evolutionarily new internal exons levels. Only 10%of genes with new TSSs gainedmouse-specific tend to have increased gene expression but only in those tissues new exons (Figure S2B), not different from the fraction of where the new exons are included inmRNAs (Figure S1A and Ta- analyzed genes overall (Figure 2B). This directional bias sug- ble S1) (Merkin et al., 2015). This trend was stronger for exons gests that the gain of species-specific new exons favors the that were efficiently spliced—assessed by percent spliced in gain of new TSSs rather than vice versa. (PSI or c) values >0.95, indicating that more than 95% of mRNAs We observed a positional effect in which the increase in the from the gene include the exon (Figure 1B)—suggesting an asso- number of TSSs per gene was associated predominantly with ciation between the extent of exon splicing and level of gene new exons located in 50 UTRs (Figure 2C). We examined the dis- expression. tribution of the locations of all mouse TSSs relative to the loca- Grouping genes by the number of promoters used, we tions of mouse-specific new exons (Figure S2C and Table S3) observed a positive association between inclusion of new exons and compared it to the distribution of rat TSSs relative to sites and gene expression for genes with multiple TSSs, but this homologous to mouse-specific exons. This comparison showed association was not observed for genes with only one TSS (Fig- an enrichment of TSSs in mouse within a few kilobases (kb) up- ure 1C). Furthermore, our RNA-seq data (from Merkin et al., stream of new exons (Figure 2D). Thus, evolutionary gain of new 2012) showed that genes with mouse-specific new exons were internal exons was specifically associated with gain of proximal, far more likely to have multiple TSSs compared to all expressed upstream TSSs. genes in mouse (Figures S1B and S1C). We confirmed that We then asked about the relationship between splicing levels genes with new mouse-specific exons are more likely to have and usage of alternative TSSswithin the same gene. Considering multiple TSSs using other methods to define TSS locations, relative TSS usage (representing the fraction of transcripts from a including H3K4me3 ChIP-seq peaks (Yu et al., 2015) and data gene that derive from a given TSS), we found that use of themost from high-resolution sequencing of polymerase-associated proximal upstream TSS (designated TSS –1) was positively RNA (Start-seq) (Scruggs et al., 2015) (Figure 1D, S1D, Table correlated with new exon inclusion, especially for TSSs located S2). Genes with rat-specific new exons (n = 1517) were also far within about 1 kb upstream of the new exon (Figures 2E and more likely to have multiple TSSs than rat genes overall (Fig- S2D). Furthermore, absolute expression of transcripts from ure S1E). Furthermore, genes that gained new species-specific nearby TSSs increased specifically in tissues in which new exons exons weremore likely to have gained TSSs in the same species, were included at moderate or high levels (Figure 2F). These ob- suggesting that the evolutionary gain of an internal exon is servations suggest a positive influence of splicing on nearby connected to evolutionary gain of TSSs in a locus (Figures 1E transcription. and S1F). To investigate this connection further, we examined the new Manipulation of Exon Splicing Impacts Upstream exons and TSSs used by a gene across different tissues. We Transcription Initiation observed that genes containing mouse-specific exons used To directly testwhether splicing impacts nearby transcription, we more TSSs than their rat orthologs (Figure S1G) and that this as- chose two candidate mouse genes, Gper1 (G protein-coupled sociation was specific tomouse tissues where the new exon was estrogen receptor 1) and Tsku (Tsukushi, small leucine rich pro- included with PSI > 0.05 (Figures 1F and S1H), showing a teoglycan). These genes both have widespread, moderate connection between splicing and TSS usage across mammalian expression and contain a mouse-specific 50 UTR internal exon organs. We also observed higher PSI values for new exons in whose splicing is positively correlated with the expression of genes with multiple alternative TSSs relative to genes with a sin- the gene across mouse tissues (Spearman r = 0.64 and 0.57, gle TSS (Figure S1I). Conversely, loss of internal exons was respectively; Figures 3A and 3B, left panels). When cultured associated with TSS loss and decreased gene expression levels mouse fibroblasts were treated with morpholino antisense oligo- (Figures 1G and S1J). Together, these observations indicate that nucleotides (MO) targeting splice sites of the new exons in these the usage of new TSSs and the splicing of new internal exons genes, exon inclusion decreased by about 4-fold in both Gper1 tend to occur in the same genes, tissues, and species, suggest- (Figure 3A) and Tsku (Figure 3B). Moreover, gene expression ing an intimate connection between splicing, increased gene levels of these two geneswere depressed to a similar extent (Fig- expression, and new TSSs. ure S3A), consistent with a positive effect of exon inclusion on gene expression. We observed similar levels of repression TSSs Arise Proximal and Upstream of New Exons when assaying metabolically labeled nascent RNA (Figures 3A We observed that increased gene expression in mouse relative and 3B) and total mRNA (Figure S3A), indicating that the effect to rat was restricted to those genes that gained TSSs in mouse is primarily at the level of transcription rather thanmRNA stability. (E) Proportion of genes that gained TSSs in mouse (mouse TSS > rat TSS), genes that lost TSSs in mouse (mouse TSS < rat TSS), and genes with same number of TSSs in both species (mouse TSS = rat TSS) for all genes expressed in both species (gray) and genes with mouse-specific new exons (pink). (F) Fold change in the number of TSSs used per gene between mouse and rat for 9 tissues for mouse genes grouped as in (C). (G) Ratio of number of TSSs used in rat over number used in mouse for all genes expressed in both species (gray) and for genes with rat-specific lost exons (blue). Statistical significance of one-way ANOVA, Tukey post hoc test, is indicated by asterisks (***p < 0.001, ****p < 0.0001), unless otherwise indicated. See also Figure S1. Cell 179, 1551–1565, December 12, 2019 1553 A *** B C ** 0.2 20 5 0.0 NS 10 0 −0.2 0 all genes genes gained −5 −0.4 TSS genes w/ genes w/ genes w/ 5’UTR coding 3’UTR mouse TSS mouse TSS mouse TSS position of new exon < rat TSS = rat TSS > rat TSS D mouse mouse rat 1.5 4 1 0 rat –1000 0 1000 position of TSSs relative to new exon (bp) 0.5 0 –1e+5 –0.5e+5 0 0.5e+5 1e+5 –100000 0 100000 position of TSSs relative to new exon (bp) E 1.0 **** F 10 ** 0.5 5 0.0 0 −0.5 −5 −1.0 < −2 –5 –1 0. 1 5 >−2 5 1 to to 2 < − – – 0 1 5 >5 t to t o o to 5 2 5 −10 − 25 5 1 .1 to to – –1 – 5 2 t t 2 5 0 1 5 to o o– t – o 5 2 5 .1 – 1 1 55 0.1 position of TSSs relative to new exon (kb) position of TSSs relative to new exon (kb) Figure 2. TSSs Arise Proximal and Upstream of New Exons (A) Fold change in gene expression between mouse and rat for genes that gained TSSs in mouse (mouse TSS > rat TSS), genes that lost TSSs in mouse (mouse TSS < rat TSS), and genes with same number of TSSs in both species (mouse TSS = rat TSS). (B) Proportion of genes that gained mouse-specific new exons from all genes expressed in mouse and rat and genes that gained new TSSs. (C) Ratio between number of TSSs used in mouse and in rat for genes with mouse-specific new exons, binned by location of the exon within the gene. (D) Histogram of TSS locations in mouse (pink) and rat (gray) in all 9 tissues for genes with mouse-specific new exons, centered on start of mouse new exon or homologous genomic position in rat. Inset shows zoom-in of locations within 1 kb of new exon. Distributions were smoothed with kernel density estimation. (E) Spearman correlations between relative TSS usage and new exon PSI across mouse tissues for TSSs binned by position relative to mouse-specific exon. (F) Difference in expression (in units of FPKM) in mouse tissues for transcripts including TSSs in tissues where new exon is moderately or highly included (PSI > 0.2) versus tissues where new exon is excluded (PSI < 0.05), grouped by TSS location relative to new exon. Statistical significance of one-way ANOVA, Tukey post hoc test, is indicated by asterisks (**p < 0.01, ***p < 0.001, ****p < 0.0001). See also Figure S2. 1554 Cell 179, 1551–1565, December 12, 2019 ρSpearrman -6 -1 (new exon ψ versus TSS ψ) TSS density (x 10 bp ) (FPKM mouse / FPKM rat), log2 proportion of genes with new exons (%) TSS FPKM tissues w/ NE ψ > 0.2 – TSS FPKM tissues w/ NE ψ < 0.05 no. TSSs mouse - no. TSSs rat A TSS –2 TSS –1 E3.CE B TSS –2 TSS –1 E2.new E3.SE E4.CEexon E2.new exon Gper1 Tsku 1.00 1.00 12 1.5 1.0 1.00 ρ = 0.57 ρ = 0.64 8 1.0 0.75 0.75 0.75 4 ** 0.5 0.50 0.50 0.5 0.50 0 ** **** 0.0 0.25 *** 0.250.25 -4 −0.5 -8 0.00 0.000.0 0.00 0.25 0.50 0.75 1.00 0.00 0.05 0.10 0.15 MO MO MO MO ψ value of new exon MO MO MO MOcontrol 5’ss control 5’ss ψ value of new exon control 5’ss control 5’ss C D 0.75 all genes TSS –1 TSS +1 TSS +2 Stoml1 ** ** 1 4 0.50 * 0.75 4 * 3 genes with mouse-specific new exons 3 0.5 0.252 *** 2 0.25 *** 1 1 0.00 0 0 0 0 2.5 5.0 7.5 10 WT 5’ss 5’ss WT 5’ss 5’ss WT 5’ss 5’ss PRO-seq reads (mouse RPKM / rat RPKM) mut 1 mut 2 mut 1 mut 2 mut 1 mut 2 TSS –1 TSS +1 TSS +2 E genes with skipped exons genes with new exons * * 0.75 1 4 4 4 0.75 *3 * 33 * 0.50 2 0.5 ** 2 2 1 0.25 1 1 0.25 tissues w/ tissues w/ new exon new exon 0 0 0 ψ < 0.05 ψ > 0.05 WT 5’ss 5’ss WT 5’ss 5’ss WT 5’ss 5’ss mut 1 mut 2 mut 1 mut 2 mut 1 mut 2 TSS -1 antisense TSS +1 antisense TSS +2 antisense 0.00 0 1 2 3 4 5 6 7 no. of sePCPA or nePCPA sites Figure 3. Manipulation of Exon Splicing Impacts Upstream Transcription Initiation (A) (Left) Relationship between fold change in gene expression betweenmouse and rat and new exon PSI value across 8 tissues forGper1 gene. (Right) qRT-PCR analysis of fold change in new exon PSI value (middle) and gene expression (right) in nascent RNAmetabolically labeled for 10minwith 5-ethynyl uridine, following treatment of NIH 3T3 cells withMO targeting new exon 50 splice site relative to control treatment.Mean±SEMof displayed distributions, n = 3 biological replicates. (B) Same as in (A) for mouse Tsku gene. (C) Fold change in nascent sense (top) and antisense (bottom) RNA levels of Stoml1 in CAD cells measured by qRT-PCR of RNA metabolically labeled for 10 min with 5-ethynyl uridine and normalized using housekeeping genes Gapdh, Hprt, and Hspcb. Wild-type cells (white) and CRISPR-Cas cells with mutations in the 50 splice site of the new exon (blue). Mean ± SEM of displayed distributions, n = 3 independent experiments. A schematic diagram of Stoml1 exon-intron or- ganization is shown at top. (D) Fold change in nascent RNA expression between mouse and rat for all genes expressed in both species (gray) and genes with mouse-specific new exons (pink) in CD4+ T cells using PRO-seq data. Distributions were smoothed with kernel density estimation. (E) Distribution of the number of polyadenylation sites used per gene located 2 kb upstream or downstream of a control set of mouse genes with skipped exons (gray, sePCPA) and genes with mouse-specific new exons (pink, nePCPA). In the inset, distribution of the number of polyadenylation sites used 2 kb upstream/ downstream of new exons per gene in tissues where new exon is excluded (PSI < 0.05, gray) or included (PSI > 0.05, pink) for genes with new exons and at least one nePCPA. Distributions are not significantly different by Kolmogorov-Smirnov test. Statistical significance of one-way ANOVA, Tukey post hoc test, is indicated by asterisks (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001). See also Figure S3. Cell 179, 1551–1565, December 12, 2019 1555 nascent RNA level nascent RNA level (FPKM mouse / FPKM rat), log2 nascent RNA level nascent RNA level ψ value relative to control nascent RNA level nascent RNA level nascent RNA level proportion of genes (FPKM mouse / FPKM rat), log2 density ψ value relative to control no. of nePCPA sites nascent RNA level A B C D E F G H Figure 4. Creation of a New Splice Site Activates the Use of a Cryptic Promoter Nearby (A) Schematic of 50 RACE products showing TSS usage defined by the fraction of clones (4, ‘‘phi’’) corresponding to each TSS in control NIH 3T3 cells and cells transfected with MO targeting the 30 and 50 splice sites of the new exon in Tsku, with a minimum of 25 clones for each sample from 2 biological replicates. *p < 0.05, one-tailed Fisher exact test. (B and C) ChIP-PCR analysis of H3K4me3 (B) and TFIIF (subunit TFIIF-alpha) (C) in Tsku gene in NIH 3T3 cells for regions indicated in (A). Mean ± SD of two independent immunoprecipitations normalized to input and mean value for control IgG antibody are shown. Data shown for control cells (gray) and cells treated with MOs targeting the 30 and 50 splice sites of the new exon (blue). (D) Luciferase activity in HeLa cells transfected with the Tsku minigene reporters (right). Promoter activities of the corresponding constructs (corrected for transfection efficiency) are presented as fold increase ofRenilla luciferase activity relative to firefly luciferase activity (both encoded on the same plasmid). Mean ± SD for n = 3 independent experiments. (legend continued on next page) 1556 Cell 179, 1551–1565, December 12, 2019 We next sought to confirm the directionality of this effect and inset and S3F). We also saw no relationship between the number to ask how splicing of new exons impacts the usage of different of nePCPA sites and gene expression changes between mouse TSSs. We chose for analysis the mouse Stoml1 (Stomatin Like 1) and rat (Figure S3G). Thus, our results suggest that effects on gene, which has three active alternative TSSs as well as a new PCPA do not contribute significantly to EMATS and that EMATS exon. Using CRISPR/Cas9 mutagenesis to generate cell lines impacts transcription initiation rather than later steps. with mutations abolishing the inclusion of the new exon (Fig- ure S3B), we observed that the three alternative TSSs of the Creation of a New Splice Site Activates the Use of a gene responded differently to inhibition of splicing of the new Cryptic Promoter Nearby exon. The upstream TSS –1 was downregulated by 4-fold, and We next sought to explore how splicing might affect the use of downstream +1 and +2 TSSs were upregulated to a similar different upstream TSSs. In the Tsku gene, the mouse-specific extent in the mutant cell lines (Figure 3C). Effects on antisense TSS in position –1 is located within 1 kb upstream of the transcription in these mutant cell lines mirrored those observed mouse-specific exon, and the conserved TSS –2 is located for sense transcription (Figure 3C), suggesting that inclusion of further upstream. Analysis by 50 RACE showed that both TSSs the new exon enhances transcription from the upstream pro- are used at similar levels in mouse fibroblasts. However, inhibit- moter in both directions. This pattern is distinct from a report ing splicing of the new exon by MO preferentially suppressed of intron-mediated enhancement in which sense-oriented in- TSS –1 (Figures 4A, S4A, and S4B). This shift was accompanied trons specifically inhibited antisense transcription (Agarwal and by a 3-fold decrease in H3K4me3 levels near TSS –1 in MO- Ansari, 2016) but is consistent with reported impacts on tran- treated cells (Figure 4B). However, levels of H3K4me3 near scription initiation resulting from changes in the position of an TSS –2 were unchanged, confirming that transcription from intron in a reporter gene (Gallegos and Rose, 2017). Levels of TSS –2 is not affected (Figure 4B). In cells treated with MOs, H3K4me3 and RNAPII decreased in the upstream TSS and levels of TFIIF and RNAPII decreased by almost 3-fold near increased in the downstream TSSs in the mutant cell lines, TSS –1 but were unchanged near TSS –2 (Figures 4C and consistent with the observed effects on nascent transcript pro- S4C). These observations suggest that splicing of the new duction (Figure S3C). exonmay contribute to recruitment of core transcriptionmachin- To assess the relationship between splicing and nearby tran- ery to the proximal TSS –1. Moreover, the loss of signal for TFIIF scription initiation on a genome-wide scale, we analyzed preci- and RNAPII near the new exon following MO treatment suggests sion run-on sequencing (PRO-seq) data from mouse and rat that inclusion of the new exon is associated with recruitment of CD4+ T cells (Danko et al., 2018). Genes with mouse-specific transcription factors, consistent with functional interactions be- new internal exons had increased nascent RNA expression in tween GTFs and splicing machinery (Damgaard et al., 2008; mouse relative to rat, compared to all expressed genes (Fig- Das et al., 2007). These observations confirm that splicing of ure 3D). Furthermore, the relative increase in nascent RNA is new exons can regulate the usage of alternative TSSs, with pre- driven by transcripts initiating upstream of the position of the dominant effects on proximal upstream promoters. new exon, specifically from TSSs within 2 kb upstream of new To dissect the impacts of individual splice sites and splicing exons (Figures S3D and S3E). levels, we created an exon corresponding to the mouse-specific PCPA can produce truncated, unstable transcripts but can be new exon in the rat Tsku gene and assessed effects on transcrip- inhibited by binding U1 snRNP near a PCPA site (Gunderson tion. In the rat Tsku locus, transcripts are predominantly et al., 1998; Kaida et al., 2010). If the observations above re- transcribed from the distal TSS –2. However, the regions homol- flected effects of splicing machinery on PCPA rather than on ogous to TSS –1 and the mouse-specific new exon have high transcription, this would require the presence of new exon prox- sequence identity with the mouse genome: both 50 splice sites imal PCPA (nePCPA) sites in affected genes. Using available are present in rat, but no consensus 30 splice site (YAG, where polyA-seq data from five mouse tissues, we observed that only Y = C or T) is present in rat near the location of the mouse 30 8.6% of genes with new exons had evidence of a nePCPA site, splice site, likely preventing splicing (Figure S4D). To introduce no higher than in control genes (Figure 3E). For the subset of the desired mutations, we cloned the 50 end of the rat Tsku genes that contain nePCPA site(s), we did not observe differ- gene upstream of the coding sequence of Renilla luciferase ences in usage of the site between tissues where the new exon and recreated the 30 splice site that is present in the mouse was spliced in and those where it was spliced out (Figure 3E, genome (50ss rn + 30ss mm), as well as a stronger 30 splice site (E) 50 RACE analysis like in (A). NIH 3T3 mouse cells were transfected with plasmids expressing the corresponding rat Tskumutants. *p < 0.05, one-tailed Fisher exact test. (F) (Top) Schematic of HNRNPU constructs are shown and (bottom) qRT-PCR analysis of fold change in new exon PSI value following treatment with a control siRNA or an siRNA targeting HNRNPU (siHNRNPU) and rescues in NIH 3T3 cells. Mean ± SEM of displayed distributions, n = 3 biological replicates. (G) ChIP-PCR analysis of TFIIF (subunit TFIIF-alpha) in Tsku gene in NIH 3T3 cells for regions indicated in (A). Mean ± SD of two independent immunoprecipi- tations normalized to input andmean value for control IgG antibody are shown. Data shown for control cells (gray) and cells treatedwith siHNRNPU (blue) rescued with HNRNPU RNA binding domain (RGG, pink) or HNRNPU full-length (FL, green). (H) Model in which creation of a splice site during evolution triggers inclusion of a new internal exon, which activates use of an upstream cryptic TSS. In themodel, exon recognition by HNRNPU in transcripts from the distal promoter recruits TFIIF that activates a TSS located proximal and upstream of the exon. Transcripts initiating from the proximal promoter also include the exon, further boosting activity of this promoter. Statistical significance of one-way ANOVA, Tukey post hoc test, is indicated by asterisks (*p < 0.05, **p < 0.01, ***p < 0.001), unless otherwise indicated. See also Figure S4. Cell 179, 1551–1565, December 12, 2019 1557 (50ss rn + 30ss stronger), while either maintaining or mutating the new exon recruits the spliceosome and perhaps also splicing native rat 50 splice site sequence (50ss mutant + 30ss mm). Strik- factors such as HNRNPU that act to recruit GTFs nearby. Tran- ingly, the creation of a 30 splice site promoted the inclusion of an scripts from the new promoter will also include the exon, further exon analogous to that observed in mouse in constructs with an activating the new promoter in a sort of positive feedback loop. intact 50 splice site (Figure S4E), indicating that this mutation is The newly activated TSS will also produce novel transcript iso- sufficient to create a new exon in the rat gene. In the presence forms and generate higher gene expression in tissues where of both 30 and 50 splice sites, but not when either splice site the upstream promoter is active and the exon is included was absent, total gene expression levels increased, as (Figure 4H). measured by luciferase activity (Figure 4D). By 50 RACE analysis, TSS –1 is used at basal levels in the minigene. However, the mouse-specific exon in the rat context activates the usage of Efficiently Spliced Exons Activate Use of Weak TSS –1 by 3-fold in the presence of a 50 splice site, demon- Proximal TSSs strating that the effect on TSS usage depends on splicing of To investigate the genomic scope of the relationship between the mouse-specific exon rather than merely the presence of a splicing and alternative TSS usage observed above, we asked 30 splice site sequence (Figures 4E and S4F). whether the inclusion of alternative skipped exons (SE) in gen- Our findings above imply the existence of mechanisms that eral—not just those that evolved recently—can influence TSS coordinate splicing with TSS usage. To explore factors that selection. We identified 49,488 SEs in mouse RNA-seq data, might be involved in this coordination, we analyzed the enrich- distributed across 13,491 genes using conservative criteria (Ta- ment of binding motifs for splicing factors in mouse novel exons. ble S4). Analyzing unique SEs with TSS-exon distances match- We observed that the binding motifs of splicing factors RBM22, ing those of new exons, we observed no significant association HNRNPU, and ELAVL1 were at least 2-fold enriched in new between SE inclusion and use of proximal upstream TSSs overall exons whose inclusion was correlated with usage of nearby (Figure 5A). In addition, we observed a symmetrical distribution TSS (r > 0.3 compared to r < 0.3) (Figure S4G). This observation of TSSs around the locations of SEs, distinct from the up- raised the possibility that some splicing factorsmay contribute to stream-biased distribution seen relative to new exons (Fig- splicing-dependent regulation of TSSs, perhaps by recruitment ure 5B). These differences suggest that genes with new exons of GTFs near sites of RNA splicing as seen above (Figure 4C). have distinct properties that favor the linkage of splicing and To explore this possibility, we analyzed the recruitment of TFIIF transcription. to the Tsku locus following depletion of HNRNPU, which is Examining other features of gene loci with new exons, we known to interact with TFIIF via its N-terminal domain (Kim and observed that, although new exons tend to have lower PSI Nikodem, 1999). HNRNPU motifs were enriched downstream values than SEs overall (Figure S5A), those new exons with of the mouse-specific exon in Tsku (Figure S4H), and splicing proximal upstream TSSs tended to have higher PSI values and of this exon was reduced by about 3-fold following depletion of stronger 50 splice sites (Figure S5B). Furthermore, although the HNRNPU (Figure 4F). Levels of TFIIF near TSS –1 and near the distribution of relative TSS usage values was similar in genes new exon decreased following HNRNPU depletion, while levels with new exons and genes with SEs generally (Figure S5C), those near TSS –2 and the constitutive exon were not affected (Fig- TSSs located proximal and upstream of new exons had lower ure 4G). Consistently, we observed downregulation of tran- relative usage across tissues than TSSs in other locations (Fig- scripts from TSS –1 and no change in transcripts from TSS –2 ure 5C). Thus, the link between splicing and TSS usage is most following HNRNPU depletion by qRT-PCR (Figure S4I). The ef- pronounced when the promoter is intrinsically weak and splicing fects of HNRNPU depletion on exon splicing, transcription activity is high. Consistently, previous studies have observed from TSS –1, and TFIIF levels at the proximal promoter were stronger intron-mediated enhancement in the presence of rescued by overexpression of full length HNRNPU. However, a weaker promoters (Callis et al., 1987). To test this idea, we truncated version of HNRNPU that contained the C-terminal grouped SEs and their most proximal and upstream TSS into RNA binding and splicing regulatory domain, but lacked the four bins from weak to strong on the basis of the relative TSS us- N-terminal domain that interacts with TFIIF, partially rescued age value, and separately for the SE PSI value, and analyzed the splicing but failed to rescue TFIIF levels and TSS –1 expression correlation between relative TSS usage and SE PSI separately (Figures 4G and S4I). Together, these observations and those within each bin. Notably, we observed that TSS usage was above (Figure 4C) suggest that HNRNPU acts both to activate most highly correlated with exon inclusion for the lowest quartile splicing of the new exon and tomediate splicing-dependent acti- of relative TSS usage values (Figures 5D, S5D, and S5E) and for vation of transcription, perhaps by recruitment of TFIIF to the the highest quartile of SE PSI (Figures 5E and S5F). This obser- proximal promoter. vation provides evidence that the EMATS observed for new In some examples studied previously, species-specific alter- exons may occur for a subset of SEs generally. Robust effects native splicing alters protein function (Gracheva et al., 2011; were observed when a weak promoter is located upstream of Gueroussov et al., 2015). Our observations support the exis- a highly included SE—an arrangement we call ‘‘EMATS organi- tence of a distinct evolutionary pathway in which, following a zation’’—which occurred in 3,833 mouse genes. The strongest mutation that generates a new internal exon, splicing of the effects were observed when the weak promoter was within 2 new exon in transcripts from a distal upstream promoter acti- kb of the SE—a pattern we call ‘‘EMATS structure’’—which vates transcription from a cryptic proximal upstream promoter. occurred in 1,777 mouse genes (Figure 5F and Table S4). In hu- A possible model consistent with our data proposes that the mans, we identified 3,548 genes with EMATS organization and 1558 Cell 179, 1551–1565, December 12, 2019 A 1.0 B C ** 8 new exons skipped 6 0.5 exons 4 0.0 4 −0.5 2 −1.0 0 < –5 – – 5 –1 1– 0 1 0 0 10 0 0 0 00 50 10 > 5 -1000 0 1000 < – –50 ,, 00 ,0 0 00 t , – 1 1 > 0 t o 1 0 t 00 00 0,0 1 0, 10 10 0 0 o t o 0 0, 00 00 0 00 0 1 t 00 ,00 to 0 t – o 0 t 0 0 0 o t – o – 1 – 00 50 o 1 t 0 0 to o 01 5 00 0 100 00 0 o , 50 0 00 0, ,0 position of TSS relative to 0 to – 10 10 0–1 10 00 0 ,000 00 0 00 0 0 0 0 00 skipped/new exon (bp) position of TSS relative to skipped exon (bp) position of TSS relative to new exon (bp) D E F 1.0 **** **** 1.0 rel. TSS use < median 0.5 0.5 0.0 0.0 rel. TSS use > median −0.5 −0.5 SE SE ψ > median ψ < median −1.0 −1.0 quartile quartile 1st 2nd 3rd 4th of rel. 1st 2nd 3rd 4th of SE ρ TSS use ψ value Spearman G E4.CE E6.CE TSS –3 TSS –2 TSS –1 E2.SE E3.SE E5.SE Zfp672 1 5 kb 1.00 1.001.00 1.00 1.00 0.75 * 0.75 * * * ** ** ** 0.50 0.50 *** ** 0.50 0.50 0.50 0.25 0.25 0.00 0.00 0.00 0.00 0.00 MO MO MO MO MO MO MO MO MO MO MO MO control E2.SE control E2.SE E2.SE E2.SE MO MO MOcontrol E3.SE control E3.SE E3.SE E3.SE control E2.SE E3.SE E4.CE E4.CE E2.SE E4.CE TSS –3 TSS –2 TSS –1 E3.SE E4.CE TSS –3 TSS –2 TSS –1 E4.CE E6.CE 8.99 bits 8.6 7.5 9.1 5.53 bits 8.6 7.5 9.1 Figure 5. Efficiently Spliced Exons Activate Weak Proximal TSSs (A) Spearman correlations between relative TSS usage (n = 49,911) and skipped exon PSE (SE, n = 13,491) in the same gene across mouse tissues for all ex- pressed TSSs in genes with SEs, binned by genomic position relative to the SE. (B) Comparison of distributions of TSS positions in 9 tissues for genes with mouse-specific new exons (blue) and genes with SEs in mouse (gray). Position 0 is set to the start coordinate of the new exon or skipped exon. Distributions were smoothed with Kernel density estimation. (C) Expression of alternative first exons (AFE) for all TSSs in genes with mouse-specific new exons in tissues where the new exon is included (PSI > 0.05), binned by position relative to the new exon. (D) Spearman correlation between relative TSS use and SE. PSI in the same gene across mouse tissues for TSSs within 1 kb upstream of the SE, binned by quartiles of mean relative TSS use. (E) Same as (D) but binned by quartiles of mean SE. PSI. (F) Heatmap showing the median Spearman correlation between relative TSS use and SE. PSI in the same gene across mouse tissues for SEs with at least one TSS located upstream, in four groups according to whether the mean relative TSS use (across tissues) and the mean SE. PSI were greater than or less than the corresponding median values (across all TSSs and SEs analyzed). (G) Exon-intron organization of mouse Zfp672 gene. qRT-PCR analysis of expression of Zfp672 in NIH 3T3 cells normalized to expression of housekeeping genes Hprt and Hspcb. Data for control cells and cells treated with MO targeting the indicated splice sites (E4.CE and E6.CE). E5.SE is not included in NIH 3T3 cells. Inclusion levels of the skipped exons, as well as levels of exon-excluding transcripts from the alternative TSSs (TSS –3, TSS –2, and TSS –1) and total gene expression are shown. Scores of 50 splice sites of skipped exons and first exons are listed in bits. Mean ± SEM of displayed distributions for n = 3 independent experiments. Statistical significance of one-way ANOVA, Tukey post hoc test, is indicated by asterisks (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001). See also Figure S5. Cell 179, 1551–1565, December 12, 2019 1559 ρSpearrman (skipped exon ψ versus rel. TSS use) relative expression ρSpearrman (skipped exon ψ versus rel. TSS use) relative expression ρ TSS density (x 10 - 6 bp- 1 ) relative expression Spearrman (skipped exon ψ versus rel. TSS use) relative expression mean TSS expression (FPKM) relative expression 0.2 0.0 −0.2 1,413 genes with EMATS structure. Considering also constitu- factors may broadly impact promoter choice and identify exten- tive exons, the number of genes increases 3-fold. sive interactionsbetween these factors andcore transcriptionma- To further investigate the distance-dependence of splicing ef- chinery, consistent with a recent study (Xiao et al., 2019). fects on TSS use, we analyzed changes in TSS usage when in- To investigate potential biological roles of gene expression hibiting the inclusion of a SE in the mouse Tsku locus that is regulation via EMATS, we analyzed the functions of genes with located more than 6 kb downstream of the TSSs. Perturbations EMATS structure. In both human and mouse, these genes of the splicing of this exon yielded no detectable changes in TSS were enriched for functions in brain development, neuron projec- usage (Figure S5G), consistent with a requirement for proximity tion, synapse organization, and related functions (Figures 6D and of the spliced exon and TSS for EMATS activity. Considering S6D). This observation raised the possibility that regulation another mouse gene, Zfp672 (Zinc Finger Protein 672)—chosen via EMATS might contribute to neuronal differentiation. For because it containedmultiple TSSs and SEs expressed inmouse example, in the Ehmt2 (Euchromatic histone-lysine N-methyl- fibroblasts—we observed that inhibition of the stronger up- transferase 2) gene, inclusion of an SE contributes to neuronal stream SE in the locus affected the usage of TSSs more dramat- differentiation (Fiszbein et al., 2016). Consistent with EMATS ically than inhibition of the weaker downstream SE (Figure 5G). A regulation of this locus, we observed that upregulation of the weaker distal TSS (TSS –2) was impacted to a similar degree as a SE during differentiation of mouse neuro2A (N2a) cells was stronger proximal TSS (TSS –1) by perturbations of the splicing accompanied by increased usage of upstreamTSSs and that us- of these SEs (Figure 5G). Together, these observations provide age of these TSSs decreased following inhibition of exon splicing further evidence that splicing of SEs can impact TSS usage, by MO (Figure S6E). To investigate whether neuro-related particularly when the TSS is intrinsically weak, the SE is highly splicing factors regulate expression via EMATS, we analyzed included, and the TSS is located proximal and upstream of the transcriptome-wide changes following depletion of PTBP1, SE. The generalization of EMATS from new exons to the much which plays a central role in neurogenesis (Linares et al., 2015), larger class of highly included SEs implies that gene expression using available ENCODE data (Van Nostrand et al., 2017). may commonly be regulated through effects on the splicing of Following PTBP1 knockdown, 758 genes had significant promoter-proximal exons. changes in SE splicing, TSS usage, and gene expression, including 255 genes with EMATS organization: a 1.7-fold enrich- Splicing Factors Impact TSS Use ment over the background frequency of EMATS genes (Fig- A mechanistic link between splicing and nearby transcription ure 6E). For example, in the human BMF (Bcl2 modifying factor) initiation could potentially be mediated by core splicing machin- gene, we observed reduced exon inclusion accompanied by ery, splicing factors, or exon junction complex (EJC) compo- decreased use of upstream proximal TSSs and decreased nents deposited during splicing, particularly those factors that gene expression following PTBP1 knockdown (Figure S6F). interact with transcription machinery, transcription factors, or chromatin. To explore functional links between RNA-binding EMATS Impacts Transcription Initiation and Translation proteins (RBPs) and TSS use, we analyzed transcriptome-wide Efficiency Globally changes in alternative TSS usage following knockdown of To investigate whether splicing of SEs affects gene expression RBPs using data from a recent ENCODE project (Van Nostrand by regulation of transcription on a genome-wide scale, we et al., 2017). Consistent with previous observations inDrosophila analyzed transient transcriptome sequencing (TT-seq) data cells showing that up to 30% of alternative splicing events that that sensitively monitors rapid changes in transcription following were affected by knockdown of 56 RBPs involved changes in stimulation of human T cells (Michel et al., 2017; Schwalb et al., promoter use (Brooks et al., 2015), our analysis detected large 2016). We observed that genes with decreased inclusion of SEs numbers of TSS changes (Figure S6A). Depletion of EJC compo- after 15 min of activation tend to have decreased transcription nents and factors involved in RNA splicing impacted larger (lower TT-seq read density) (Figure 7A). This trend was stronger numbers of TSSs than did depletion of other RBPs (Figure 6A). in genes with EMATS structure (Figure 7A), suggesting that Previous studies have linked the EJC to gene expression regula- EMATS contributes to gene regulation by modulation of tran- tion, and to the control of RNAPII pausing (Akhtar et al., 2019; scription globally. Wiegand et al., 2003). Switching between alternative 50 UTR isoforms by altered TSS Basedonour results indicating that splicing factors can regulate choice has recently been identified as an important regulator of the recruitment of transcription machinery near alternative exons, translation efficiency in yeast (Cheng et al., 2018). To answer we focused on the 10 splicing factors associated with the largest whether the splicing-dependent regulation of TSS selection by numbers of changes in TSS usage (Figures 6B and S6B), which EMATS impacts translation, we analyzed transcript isoforms in included theHNRNPUfactor studiedabove.Usingprotein-protein polysomes sequencing (TrIP-seq) data from human HEK293T interaction (PPI) data fromtheSTRINGdatabase (Szklarczyket al., cells (Floor and Doudna, 2016) to assess the ribosome occu- 2015), we observed that these ten splicing factors interact with 65 pancy of EMATS isoforms. We found that first exons with other proteins, including subunits of RNAPII and GTFs such as EMATS-associated TSSs are significantly more ribosome-asso- TFIIF (Figure 6C). Compared with the PPI partners of the 10 ciated by a median fold change of about 1.3-fold than gene- splicing factors whose depletion affected the fewest TSSs, these matched controls that lack EMATS structure (Figures 7B and 65 proteins were enriched for functions in enhancer binding, tran- S7A). These observations indicate that isoforms activated tran- scription factor activity, and promoter proximal binding (Fig- scriptionally by EMATS tend to have enhanced translational ac- ureS6C). Together, theseobservations indicate that somesplicing tivity, amplifying the impact of EMATS on protein production. 1560 Cell 179, 1551–1565, December 12, 2019 A B exon junction complex regulation of RNA splicing mRNA polyadenylation spliceosomal complex 9 regulation of translation mitochondrial RNA regulation ribosome & translation factors 6 RNA localization RNA stability & decay RNA export from nucleus 3 microRNA processing & regulation snoRNA / snRNA / telomere processing 0 50 100 150 0 0 100 200 300 no. of genes with altered TSS use no. of genes with altered TSS use following KD C seed of 10 splicing factors transcription from RNAPII promoter enhancer-box binding RNAPII promoter binding transcription co-factor binding DNA sequence-specific binding gene expression regulators snRNA binding pre-mRNA binding splicing factors RNA binding D E synapse assembly * actomyosin structure organization * genes with genes with changes in SE changes in GE peptidyl−serine phosphorylation ** 1493235 1075 dendrite development *** striated muscle cel differentiation ** 758 503 255 3293 neuron projection development 222 2710** neuron projection morphogenesis **** 1031 genes with changes in brain development **** SE, GE and TSSs genes with 0 1 2 3 genes with EMATS organization changes in TSSs fold enrichment of GO category all overlaps p-value < 0.00001 Figure 6. A Subset of Splicing Factors Have Wide Impacts on TSS Usage and Interact with Transcription Machinery (A) Distribution of the number of genes with significant changes in promoter usage associated with depletion of 250 RBPs, binned by Gene Ontology Biological Process categories of RBPs. Mean ± SEM between all RBPs in each Gene Ontology (GO) category for two cell lines (HepG2 and K562) is plotted. (B) Histogram of number of geneswith significant changes in TSS usage following depletion of 67 splicing factors. Top ten splicing factors with greatest number of changes shown in red. (C) PPI network for the top 10 splicing factors from (B), colored by Gene Ontology category. Nodes represent proteins and edges represent PPIs. Node and label size are proportional to protein connectivity. The 10 selected splicing factors in red primarily interact with other 65 proteins, generating a network with 75 nodes and 424 edges. PPI data are from STRING database (Szklarczyk et al., 2015). Networks were built using Gephi (http://gephi.org/). (D) Gene Ontology analysis of 1,777 mouse genes with the strongest EMATS potential. Fold enrichments shown for the most significant categories with asterisk indicating adjusted p-values and color indicating relation to neuron development. (E) (Left) Venn diagram showing the overlap between genes with significant changes in gene expression (GE), alternative splicing of SEs, and relative usage of TSSs following knockdown of PTBP1 in human HepG2 cells. (Right) Venn diagram showing the overlap between genes with changes in GE, SE, and TSSs following knockdown of PTBP1 for human genes with EMATS organization. The overlap is 1.7-fold above background expectation (p < 1.6e-20, hypergeo- metric test). See also Figure S6. Cell 179, 1551–1565, December 12, 2019 1561 category of RBP knocked down no. of splicing factors A **** C TF 5000 0 enhancer direct tx −5000 all genes genes with EMATS GTF RNAPII decreased genes with inclusion decreased of SE inclusion SF of SE B EMATS **** EMATS TSS4 control TSS direct tx reg. **** AAA3 AAA TF SF A 2 AA AAA direct tx reg. EMATS 1 AAA stimulus 0 High polysomes Low polysomes Figure 7. EMATS Impacts Transcription and Translation Initiation Globally (A) Change in nascent RNA levels (TT-seq read counts) after 15 min of T cell stimulation for all genes expressed in humans, genes with increased inclusion of a skipped exon (SE), and genes with EMATS structure and increased inclusion of SE. (B) EMATS TSS isoforms (pink) have increased translation efficiency (TE) relative to matched control isoforms from the same gene (gray). For each AFE, the TE was calculated across high (left, polysomes 6–8) or low (right, monosomes and polysomes 2–4) polysomes. Boxplots show the distributions of TE values. (C) Model for the role of EMATS in dynamic gene expression programs. Growth factor or other stimuli activate transcription factors (TF) and splicing factors (SF). TFs influence gene expression by direct effects on transcription (tx) and indirectly by regulating levels of SFs. Effects of SFs on splicing contribute to gene expression programs by EMATS. In genes with EMATS structure, SFs recruit general TFs (GTF) or RNAPII to activate weak TSS(s) proximal and upstream of the exon. Statistical significance of one-way ANOVA, Tukey post hoc test, is indicated by asterisks (****p < 0.0001), unless otherwise indicated. See also Figure S7. DISCUSSION machinery can recruit GTFs or modulate transcription activity (Damgaard et al., 2008; Fong and Zhou, 2001; Kwek et al., Here, we have shown that inclusion of a new internal exon in a 2002), and depletion of RBPs can impact promoter selection gene can activate transcription from an upstream TSS and on a large scale (Brooks et al., 2015) (Figure 6B). The involve- thereby increase gene expression levels, a phenomenon which ment of splicing machinery or proteins deposited on the tran- we refer to as EMATS. Our study highlights several features of script in connection with splicing would explain feature (1) this relationship: (1) it requires exon splicing, not merely pres- above, while the more efficient recruitment of splicing machin- ence of a 50 or 30 splice site; (2) it is more potent when the ery to more efficiently spliced exons would explain feature (2). exon is highly included; (3) it is more potent when the promoter Recruitment of RNAPII or GTFs might be expected to activate is intrinsically weak; (4) it is sensitive to genomic distance, occur- transcription more effectively at weaker promoters where ring most robustly when exon and promoter are within 1-2 kb; RNAPII recruitment is limiting than at strong promoters with and (5) the above features occur in thousands of mammalian higher intrinsic RNAPII occupancy, explaining feature (3). A genes (Table S4). requirement for direct physical interaction between splicing ma- The most straightforward model to explain the above proper- chinery and RNAPII or GTFs might constrain the genomic dis- ties would involve direct positive effects of cotranscriptionally tances involved, explaining feature (4). However, the varied recruited splicing components on recruitment of transcription chromatin conformations of different gene loci—which in machinery to nearby upstream promoters (Figure 7D). Splicing some cases may involve chromatin loops between promoters 1562 Cell 179, 1551–1565, December 12, 2019 translation efficiency Change in expression 15 minutes TrIP-seq reads (polysome / total cytoplasmic) post-stimulation, TT-seq and alternative exons (Mercer et al., 2013; Ruiz-Velasco et al., d QUANTIFICATION AND STATISTICAL ANALYSIS 2017)—might alter distance requirements for different genes. B Definition of species-specific exons Frequent occurrence of the evolutionary path outlined above B Transcription start site annotation (Figure 4I) and of alternative 50 UTRs (Singer et al., 2008) may ac- B Software for data analysis, graphical plots and statisti- count for the widespread occurrence of EMATS organization in cal analyses mammalian genomes, explaining feature (5). d DATA AND CODE AVAILABILITY Recent studies have broadened the definition of transcrip- B Data availability tional enhancers, showing that some promoters also function d ADDITIONAL RESOURCES as enhancers (Engreitz et al., 2016; Scruggs et al., 2015); our B New exon inclusion, TSS usage, and species-specific findings support further broadening of this definition to include expression some exons as well. Here we have focused on the conse- B Definition of new exon-proximal cleavage and polya- quences of regulating exon splicing; it is possible that regulatory denylation sites change in the splicing of 50 UTR introns may also impact the use B Effects on nascent and steady state RNA levels of upstream promoters. We propose that emergence of new internal exons and of new SUPPLEMENTAL INFORMATION TSSs are linked (Figure 4I). Once activated, the new TSS pro- duces new transcript isoforms and higher overall expression of Supplemental Information can be found online at https://doi.org/10.1016/j. cell.2019.11.002. the gene in specific tissues, providing a substrate for the regula- tory evolution of the gene. The most obvious regulatory role for ACKNOWLEDGMENTS EMATS would be as a means for splicing factors to contribute to gene expression programs involved in differentiation or We thank Karla Neugebauer, Alberto R. Kornblihtt, Phillip A. Sharp, and mem- cellular responses to stimuli (Figure 7C). Specifically, we pro- bers of the Burge lab for helpful discussions and comments, as well as Maria pose that external stimuli such as growth factors trigger gene Alexis, Jason Merkin, Peter Freese, Charles Danko, and Brenton R. Graveley expression changes not only via direct effects on TF activity for assistance with access to genomic data and analysis pipelines. This (Malladi et al., 2016; Rajbhandari et al., 2018) but also by work was supported by grants from the NIH to C.B.B. (grant numbers HG002439 and GM085319) and by the Pew Latin American Fellows Program changes in splicing factor levels downstream of affected TFs in the Biomedical Sciences to A.F. or effects on splicing factor activity (Reinhardt et al., 2011; van der Houven van Oordt et al., 2000), triggering additional gene AUTHOR CONTRIBUTIONS expression changes via EMATS. Another implication of our find- ings is that targeted activation of the expression of a gene for A.F. and C.B.B. designed the study and wrote the manuscript. A.F. conducted research or therapeutic purposes may be achievable by use of all computational analyses and designed and performed experiments. K.S.K. compounds such as antisense oligonucleotides or small mole- performed experiments in Figure 4. B.E.B. performed the polysome profile cules (Havens and Hastings, 2016) that enhance the splicing of analysis in Figure 7B. C.B.B. supervised the project. appropriately located alternative or cryptic promoter-prox- DECLARATION OF INTERESTS imal exons. C.B.B. and A.F. have filed a patent application related to this work. The authors STAR+METHODS declare no other competing interests. Received: April 2, 2019 Detailed methods are provided in the online version of this paper Revised: August 20, 2019 and include the following: Accepted: October 30, 2019 Published: November 28, 2019 d KEY RESOURCES TABLE d LEAD CONTACT AND MATERIALS AVAILABILITY REFERENCES B Materials Availability Statement d EXPERIMENTAL MODEL AND SUBJECT DETAILS Agarwal, N., and Ansari, A. (2016). Enhancement of Transcription by a B Cell lines, cell culture and treatments Splicing-Competent Intron Is Dependent on Promoter Directionality. PLoS d METHOD DETAILS Genet. 12, e1006047. B RNA-seq analysis and genome builds Akhtar, J., Kreim, N., Marini, F., Mohana, G., Brüne, D., Binder, H., and Roig- B CRISPR sgRNA design, genetic deletions and nant, J.-Y. (2019). Promoter-proximal pausing mediated by the exon junction complex regulates splicing. Nat. Commun. 10, 521. genotyping RNA Extraction, RT-PCR and qPCR Almada, A.E., Wu, X., Kriz, A.J., Burge, C.B., and Sharp, P.A. (2013). PromoterB directionality is controlled by U1 snRNP and polyadenylation signals. Nature B ChIP and antibodies 499, 360–363. B 50 RACE Anvar, S.Y., Allard, G., Tseng, E., Sheynkman, G.M., de Klerk, E., Vermaat, M., B Plasmids, RNAi and luciferase activity assay Yin, R.H., Johansson, H.E., Ariyurek, Y., den Dunnen, J.T., et al. (2018). Full- B PRO-seq data analysis length mRNA sequencing uncovers a widespread coupling between transcrip- B Motif enrichment analysis tion initiation and mRNA processing. Genome Biol. 19, 46. B TT-seq data analysis Bentley, D.L. (2014). Coupling mRNA processing with transcription in time and B Polysome profile analysis space. Nat. Rev. Genet. 15, 163–175. Cell 179, 1551–1565, December 12, 2019 1563 Bieberstein, N.I., Carrillo Oesterreich, F., Straube, K., and Neugebauer, K.M. Gracheva, E.O., Cordero-Morales, J.F., González-Carcacı́a, J.A., Ingolia, N.T., (2012). First exon length controls active chromatin signatures and transcrip- Manno, C., Aranguren, C.I., Weissman, J.S., and Julius, D. (2011). Ganglion- tion. Cell Rep. 2, 62–68. specific splicing of TRPV1 underlies infrared sensation in vampire bats. Nature Braunschweig, U., Gueroussov, S., Plocik, A.M., Graveley, B.R., and Blen- 476, 88–91. cowe, B.J. (2013). Dynamic integration of splicing within gene regulatory path- Gueroussov, S., Gonatopoulos-Pournatzis, T., Irimia, M., Raj, B., Lin, Z.Y., ways. Cell 152, 1252–1269. Gingras, A.C., and Blencowe, B.J. (2015). An alternative splicing event am- Brooks, A.N., Duff, M.O., May, G., Yang, L., Bolisetty, M., Landolin, J., Wan, K., plifies evolutionary differences between vertebrates. Science 349, 868–873. Sandler, J., Booth, B.W., Celniker, S.E., et al. (2015). Regulation of alternative Gunderson, S.I., Polycarpou-Schwarz, M., and Mattaj, I.W. (1998). U1 snRNP splicing in Drosophila by 56 RNA binding proteins. Genome Res. 25, inhibits pre-mRNA polyadenylation through a direct interaction between U1 1771–1780. 70K and poly(A) polymerase. Mol. Cell 1, 255–264. Callis, J., Fromm, M., and Walbot, V. (1987). Introns increase gene expression Havens, M.A., and Hastings, M.L. (2016). Splice-switching antisense oligonu- in cultured maize cells. Genes Dev. 1, 1183–1200. cleotides as therapeutic drugs. Nucleic Acids Res. 44, 6549–6563. Chathoth, K.T., Barrass, J.D., Webb, S., and Beggs, J.D. (2014). A splicing- Kaida, D., Berg, M.G., Younis, I., Kasim, M., Singh, L.N., Wan, L., and Drey- dependent transcriptional checkpoint associated with prespliceosome forma- fuss, G. (2010). U1 snRNP protects pre-mRNAs from premature cleavage tion. Mol. Cell 53, 779–790. and polyadenylation. Nature 468, 664–668. Cheng, Z., Otto, G.M., Powers, E.N., Keskin, A., Mertins, P., Carr, S.A., Jova- Kameoka, S., Duque, P., and Konarska, M.M. (2004). p54(nrb) associates with novic, M., and Brar, G.A. (2018). Pervasive, Coordinated Protein-Level the 50 splice site within large transcription/splicing complexes. EMBO J. 23, Changes Driven by Transcript Isoform Switching during Meiosis. Cell 172, 1782–1791. 910–923.e16. Katz, Y., Wang, E.T., Airoldi, E.M., and Burge, C.B. (2010). Analysis and design Custódio, N., and Carmo-Fonseca, M. (2016). Co-transcriptional splicing and of RNA sequencing experiments for identifying isoform regulation. Nat. the CTD code. Crit. Rev. Biochem. Mol. Biol. 51, 395–411. Methods 7, 1009–1015. Damgaard, C.K., Kahns, S., Lykke-Andersen, S., Nielsen, A.L., Jensen, T.H., and Kjems, J. (2008). A 50 splice site enhances the recruitment of basal tran- Kim, M.K., and Nikodem, V.M. (1999). hnRNP U inhibits carboxy-terminal scription initiation factors in vivo. Mol. Cell 29, 271–278. domain phosphorylation by TFIIH and represses RNA polymerase II elonga- tion. Mol. Cell. Biol. 19, 6833–6844. Danko, C.G., Choate, L.A., Marks, B.A., Rice, E.J., Wang, Z., Chu, T., Martins, A.L., Dukler, N., Coonrod, S.A., Tait Wojno, E.D., et al. (2018). Dynamic evolu- Kornblihtt, A.R., Schor, I.E., Alló, M., Dujardin, G., Petrillo, E., and Muñoz, M.J. tion of regulatory element ensembles in primate CD4+ T cells. Nat. Ecol. Evol. (2013). Alternative splicing: a pivotal step between eukaryotic transcription 2, 537–548. and translation. Nat. Rev. Mol. Cell Biol. 14, 153–165. Das, R., Yu, J., Zhang, Z., Gygi, M.P., Krainer, A.R., Gygi, S.P., and Reed, R. Kwek, K.Y., Murphy, S., Furger, A., Thomas, B., O’Gorman, W., Kimura, H., (2007). SR proteins function in coupling RNAP II transcription to pre-mRNA Proudfoot, N.J., and Akoulitchev, A. (2002). U1 snRNA associates with TFIIH splicing. Mol. Cell 26, 867–881. and regulates transcriptional initiation. Nat. Struct. Biol. 9, 800–805. Derti, A., Garrett-Engele, P., Macisaac, K.D., Stevens, R.C., Sriram, S., Chen, Lawrence, M., Huber, W., Pagès, H., Aboyoun, P., Carlson, M., et al. (2013). R., Rohl, C.A., Johnson, J.M., and Babak, T. (2012). A quantitative atlas of pol- Software for Computing and Annotating Genomic Ranges. PLOS Computa- yadenylation in five mammals. Genome Res. 22, 1173–1183. tional Biology 9, e1003118. Dobin, A., Davis, C.A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., Batut, Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., Marth, G., P., Chaisson,M., andGingeras, T.R. (2013). STAR: ultrafast universal RNA-seq Abecasis, G., and Durbin, R. (2009). The Sequence Alignment/Map format and aligner, 29, pp. 15–21. SAMtools. Bioinformatics 25, 2078–2079. Dominguez, D., Freese, P., Alexis, M.S., Su, A., Hochman, M., Palden, T., Ba- Linares, A.J., Lin, C.-H., Damianov, A., Adams, K.L., Novitch, B.G., and Black, zile, C., Lambert, N.J., Van Nostrand, E.L., Pratt, G.A., et al. (2018). Sequence, D.L. (2015). The splicing regulator PTBP1 controls the activity of the transcrip- Structure, and Context Preferences of Human RNA Binding Proteins. Mol. Cell tion factor Pbx1 during neuronal differentiation. eLife 4, e09268. 70, 854–867.e9. Malladi, S., Macalinao, D.G., Jin, X., He, L., Basnet, H., Zou, Y., de Stanchina, Emili, A., Shales, M., McCracken, S., Xie, W., Tucker, P.W., Kobayashi, R., E., and Massagué, J. (2016). Metastatic Latency and Immune Evasion through Blencowe, B.J., and Ingles, C.J. (2002). Splicing and transcription-associated Autocrine Inhibition of WNT. Cell 165, 45–60. proteins PSF and p54nrb/nonO bind to the RNA polymerase II CTD, 8, pp. 1102–1111. Mercer, T.R., Edwards, S.L., Clark, M.B., Neph, S.J., Wang, H., Stergachis, A.B., John, S., Sandstrom, R., Li, G., Sandhu, K.S., et al. (2013). DNase I-hy- Engreitz, J.M., Haines, J.E., Perez, E.M., Munson, G., Chen, J., Kane, M., persensitive exons colocalize with promoters and distal regulatory elements. McDonel, P.E., Guttman, M., and Lander, E.S. (2016). Local regulation of Nat. Genet. 45, 852–859. gene expression by lncRNA promoters, transcription and splicing. Nature 539, 452–455. Merkin, J., Russell, C., Chen, P., and Burge, C.B. (2012). Evolutionary dy- namics of gene and isoform regulation in Mammalian tissues. Science 338, Fiszbein, A., Giono, L.E., Quaglino, A., Berardino, B.G., Sigaut, L., von Bilderl- 1593–1599. ing, C., Schor, I.E., Steinberg, J.H., Rossi, M., Pietrasanta, L.I., et al. (2016). Alternative Splicing of G9a Regulates Neuronal Differentiation. Cell Rep. 14, Merkin, J.J., Chen, P., Alexis, M.S., Hautaniemi, S.K., and Burge, C.B. (2015). 2797–2808. Origins and impacts of new mammalian exons. Cell Rep. 10, 1992–2005. Floor, S.N., and Doudna, J.A. (2016). Tunable protein synthesis by transcript Michel, M., Demel, C., Zacher, B., Schwalb, B., Krebs, S., Blum, H., Gagneur, isoforms in human cells. eLife 5, 1276. J., and Cramer, P. (2017). TT-seq captures enhancer landscapes immediately Fong, Y.W., and Zhou, Q. (2001). Stimulatory effect of splicing factors on tran- after T-cell stimulation. Mol. Syst. Biol. 13, 920. scriptional elongation. Nature 414, 929–933. Morris, D.P., and Greenleaf, A.L. (2000). The splicing factor, Prp40, binds the Furger, A., O’Sullivan, J.M., Binnie, A., Lee, B.A., and Proudfoot, N.J. (2002). phosphorylated carboxyl-terminal domain of RNA polymerase II. J. Biol. Promoter proximal splice sites enhance transcription. Genes Dev. 16, Chem. 275, 39935–39943. 2792–2799. Mortillaro, M.J., Blencowe, B.J., Wei, X., Nakayasu, H., Du, L., Warren, S.L., Gallegos, J.E., and Rose, A.B. (2017). Intron DNA Sequences Can Be More Sharp, P.A., and Berezney, R. (1996). A hyperphosphorylated form of the large Important Than the Proximal Promoter in Determining the Site of Transcript subunit of RNA polymerase II is associated with splicing complexes and the Initiation. Plant Cell 29, 843–853. nuclear matrix. Proc. Natl. Acad. Sci. USA 93, 8253–8257. 1564 Cell 179, 1551–1565, December 12, 2019 Neugebauer, K.M., andRoth,M.B. (1997). Transcription units as RNAprocess- Shen, S., Park, J.W., Lu, Z.-X., Lin, L., Henry, M.D., Wu, Y.N., Zhou, Q., and ing units. Genes Dev. 11, 3279–3285. Xing, Y. (2014). rMATS: robust and flexible detection of differential alternative Paten, B., Herrero, J., Beal, K., Fitzgerald, S., andBirney, E. (2008). Enredo and splicing from replicate RNA-Seq data. Proc. Natl. Acad. Sci. USA 111, Pecan: genome-wide mammalian consistency-based multiple alignment with E5593–E5601. paralogs. Genome Res. 18, 1814–1828. Singer, G.A.C., Wu, J., Yan, P., Plass, C., Huang, T.H.-M., and Davuluri, R.V. Paz, I., Kosti, I., Ares, M., Jr., Cline, M., and Mandel-Gutfreund, Y. (2014). (2008). Genome-wide analysis of alternative promoters of human genes using RBPmap: a web server for mapping binding sites of RNA-binding proteins. Nu- a custom promoter tiling array. BMC Genomics 9, 349. cleic Acids Res. 42, W361–W367. Szklarczyk, D., Franceschini, A., Wyder, S., Forslund, K., Heller, D., Huerta- Quinlan, A.R., and Hall, I.M. (2010). BEDTools: a flexible suite of utilities for Cepas, J., Simonovic, M., Roth, A., Santos, A., Tsafou, K.P., et al. (2015). comparing genomic features. Bioinformatics 26, 841–842. STRING v10: protein-protein interaction networks, integrated over the tree of Rajbhandari, P., Thomas, B.J., Feng, A.-C., Hong, C., Wang, J., Vergnes, L., life. Nucleic Acids Res. 43, D447–D452. Sallam, T., Wang, B., Sandhu, J., Seldin, M.M., et al. (2018). IL-10 Signaling Trapnell, C., Pachter, L., and Salzberg, S.L. (2009). TopHat: discovering splice Remodels Adipose Chromatin Architecture to Limit Thermogenesis and En- junctions with RNA-Seq. Bioinformatics 25, 1105–1111. ergy Expenditure. Cell 172, 218–233.e17. Trapnell, C., Roberts, A., Goff, L., Pertea, G., Kim, D., Kelley, D.R., Pimentel, Ran, F.A., Hsu, P.D., Wright, J., Agarwala, V., Scott, D.A., and Zhang, F. (2013). H., Salzberg, S.L., Rinn, J.L., and Lior, P. (2012). Differential gene and tran- Genome engineering using the CRISPR-Cas9 system. Nat. Protoc. 8, script expression analysis of RNA-seq experiments with TopHat and Cufflinks. 2281–2308. Nature protocols 7, 562–578. Reinhardt, H.C., Cannell, I.G., Morandell, S., and Yaffe, M.B. (2011). Is post- van der Houven van Oordt, W., Diaz-Meco, M.T., Lozano, J., Krainer, A.R., transcriptional stabilization, splicing and translation of selective mRNAs a Moscat, J., and Cáceres, J.F. (2000). The MKK(3/6)-p38-signaling cascade al- key to the DNA damage response? Cell Cycle 10, 23–27. ters the subcellular distribution of hnRNP A1 and modulates alternative Reyes, A., and Huber, W. (2018). Alternative start and termination sites of tran- splicing regulation. J. Cell Biol. 149, 307–316. scription drive most transcript isoform differences across human tissues. Nu- Van Nostrand, E.L., Freese, P., Pratt, G.A., Wang, X., Wei, X., Xiao, R., Blue, cleic Acids Res. 46, 582–592. S.M., Chen, J.-Y., Cody, N.A., Dominguez, D., et al. (2017). A Large-Scale Roberts, A., Pimentel, H., Trapnell, C., and Pachter, L. (2011). Identification of Binding and Functional Map of Human RNA Binding Proteins. bioRxiv. novel transcripts in annotated genomes using RNA-Seq. Bioinformatics 27, https://doi.org/10.1101/179648. 2325–2329. Vincent, M., Lauriault, P., Dubois, M.F., Lavoie, S., Bensaude, O., and Chabot, Ruiz-Velasco, M., Kumar,M., Lai, M.C., Bhat, P., Solis-Pinson, A.B., Reyes, A., B. (1996). The nuclear matrix protein p255 is a highly phosphorylated form of Kleinsorg, S., Noh, K.-M., Gibson, T.J., and Zaugg, J.B. (2017). CTCF-Medi- RNA polymerase II largest subunit which associates with spliceosomes. Nu- atedChromatin Loops between Promoter andGene Body Regulate Alternative cleic Acids Res. 24, 4649–4652. Splicing across Individuals. Cell Syst. 5, 628–637.e6. Wiegand, H.L., Lu, S., and Cullen, B.R. (2003). Exon junction complexes Schor, I.E., Fiszbein, A., Petrillo, E., and Kornblihtt, A.R. (2013). Intragenic mediate the enhancing effect of splicing on mRNA expression. Proc. Natl. epigenetic changes modulate NCAM alternative splicing in neuronal differen- Acad. Sci. USA 100, 11327–11332. tiation. EMBO J. 32, 2264–2274. Xiao, R., Chen, J.-Y., Liang, Z., Luo, D., Chen, G., Lu, Z.J., Chen, Y., Zhou, B., Schwalb, B., Michel, M., Zacher, B., Frühauf, K., Demel, C., Tresch, A., Ga- Li, H., Du, X., et al. (2019). Pervasive Chromatin-RNA Binding Protein Interac- gneur, J., and Cramer, P. (2016). TT-seq maps the human transient transcrip- tions Enable RNA-Based Regulation of Transcription. Cell 178, 107–121.e18. tome. Science 352, 1225–1228. Yeo, G., and Burge, C.B. (2004). Maximum entropy modeling of short Scruggs, B.S., Gilchrist, D.A., Nechaev, S., Muse, G.W., Burkholder, A., Fargo, sequence motifs with applications to RNA splicing signals. J. Comput. Biol. D.C., and Adelman, K. (2015). Bidirectional Transcription Arises from Two 11, 377–394. Distinct Hubs of Transcription Factor Binding and Active Chromatin. Mol. Yu, H.-B., Yurieva, M., Balachander, A., Foo, I., Leong, X., Zelante, T., Zolezzi, Cell 58, 1101–1112. F., Poidinger,M., andRicciardi-Castagnoli, P. (2015). NFATc2mediates epige- Shaul, O. (2017). How introns enhance gene expression. Int. J. Biochem. Cell netic modification of dendritic cell cytokine and chemokine responses to dec- Biol. 91 (Pt B), 145–155. tin-1 stimulation. Nucleic Acids Res. 43, 836–847. Cell 179, 1551–1565, December 12, 2019 1565 STAR+METHODS KEY RESOURCES TABLE REAGENT or RESOURCE SOURCE IDENTIFIER Antibodies H3K4me3 antibody Invitrogen PA5-17420; RRID: AB_10977680 RNA polymerase II (CTD repeat YSPTSPS) antibody Abcam Ab817; RRID: AB_306327 Transcription Factor IIF1 (TFIIF-alpha) antibody Invitrogen PA5-30050; RRID: AB_2547524 Rabbit IgG antibody Cell Signaling Technology 2729; RRID: AB_1031062 Chemicals, Peptides, and Recombinant Proteins M-MLV reverse transcriptase Invitrogen 28025013 Lipofectamine 2000 Life Technologies 11668027 Dual-Luciferase Reporter Assay System Promega E1910 DMEM GIBCO 11965118 Fetal bovine serum (FBS) GIBCO A31406-02 Endo-Porter (in PEG) Gene Tools OT-EP-PEG-1 Lipofectamine RNAiMAX Thermo 13778150 DMEM/F12 GIBCO 11320033 Critical Commercial Assays RNA-easy kit QIAGEN 74104 Click-iT Invitrogen C10365 50 RACE System for Rapid Amplification of cDNA Ends Invitrogen 18374058 Experimental Models: Cell Lines NIH 3T3 ATCC CRL-1658 HeLa ATCC CCL-2 CAD ECACC 08100805 N2a ATCC CCL-131 Primers for qPCR experiments Tsku_INC_F This study GTGTCCTGCCAAAGCAAGTG Tsku_INC_R This study CAGGAACAGAGAGCACAGCA Tsku_INC_F_juntCE This study CCTGGCTGAGCAGGTGT Tsku_INC_R_juntCE This study ATCCAAAGGGATGGGCACAG Tsku_INC_TSS1_F This study GACCTGCCAGGACGCTG Tsku_INC_TSS1_R This study TCAGCCAGGTCTGCTCCTAT Tsku_antisense_TSS1_F This study GCTCAGGGAGCGTCGTTAAA Tsku_antisense_TSS1_R This study GGGAACCGCGCACTTTTTAG Tsku_INC_TSS2_F This study TGGCCAGGCTCAGAGGAC Tsku_INC_TSS2_R This study TCAGCCAGGTCTGCTCCTAT Tsku_antisense_TSS2_F This study CCCCTTTTCATCACAGCCCA Tsku_antisense_TSS2_R This study GGGACGAACCTTCCAATCCA Stoml1_TSS1_F This study GGAGTAAAGCCGGAAGCAGT Stoml1_TSS1_R This study TCATGCTTGGAAGGTCTGGC Stoml1_TSS2_F This study ATATGGGACCTCCGTGTCCA Stoml1_TSS2_R This study AGCATGCCACACTCCTTACC Stoml1_TSS3_F This study CATGCAGAGCACTGACCTAGT Stoml1_TSS3_R This study CTGGAGGCTGTACTCAAGGC Stoml1_TSS1_antisense_F This study TCCTGACCACCTCCTACCTG Stoml1_TSS1_antisense_R This study TGGCCTCAAACCATTCCTCC (Continued on next page) e1 Cell 179, 1551–1565.e1–e6, December 12, 2019 Continued REAGENT or RESOURCE SOURCE IDENTIFIER Stoml1_TSS2_antisense_F This study CTGGGGAGAACTGAGGGTTC Stoml1_TSS2_antisense_R This study GAACCCCAGAGGGGAGTCTAT Stoml1_TSS3_antisense_F This study CTGCCTCTTGATTCCCAGCA Stoml1_TSS3_antisense_R This study CCCTTCCAAGACTGTGGCTT Primers for 50 RACE experiments Tsku_GSP1 This study TGAATGGTAGGTGCAGGCAG Tsku_GSP2 This study GGGAAGCAGGCGATGGATAA Tsku_nestedGSP This study GATGTCACTCAAGGGGGAGC hRLUC_GST1 This study GAACCAAGCGGTGAGGTACT hRLUC_GST2 This study CGATATGAGCCATTCCCGCT hRLUC_nested This study ATGATGCATCTAGCCACGGG Primers for ChIP experiments Tsku_Promoter_F This study ACTTTAACGACGCTCCCTGA Tsku_Promoter_R This study ATGGGCCGGCGCTTTT Tsku_TSS1_Intron1_F This study GAGGCGACAACTGCAGACC Tsku_TSS1_Intron1_R This study CGACTCTATGGCTCGGTGTC Tsku_Intron1_TSS2_F This study TTCCCAAGGGATGGCCAATG Tsku_Intron1_TSS2_R This study AGTGACCGAATCTCAACGGG Tsku_TSS2_Intron2_F This study GTGGCGAGCTTAGCTGAAAG Tsku_TSS2_Intron2_R This study ACCCAGGATCAAAAGCTCGG Tsku_Intron2_NEx_F This study ACAGACTCGGCAAGAGATGGA Tsku_Intron2_NEx_R This study CTTCAGGAAACTCCCAGGCTCA Tsku_NEx_F This study ACGCTGAGCCTGGGAGTTTC Tsku_NEx_R This study TAGCACTTGCTTTGGCAGGA Tsku_NEx_Intron3_F This study GTCCCATAGGAGCAGACCTGG Tsku_NEx_Intron3_R This study TCCCAGCCTTTGGGTAACTC Tsku_Intron3_AltEx_F This study GCTCAGTTCTCCCTTAGTGGG Tsku_Intron3_AltEx_R This study TGGGGGCTTCATTCACCTTT Tsku_AltEx_Intron4_F This study ACCGTCCGGTCTAACAGATTT Tsku_AltEx_Intron4_R This study ACGGTTAAGGGTTGGACCAG Tsku_LastEx_F This study AGGGCATCCTCCATCTACCA Tsku_LastEx_R This study GCAAACCCAGGCCTGAAAAC Morpholino sequences G9a_mouse_E10_5ss This study GTCCCGGCAGTTGGCAATTAATTAC G9a_mouse_E10_3ss This study CCATTCACTCCTGACACAGAGACAG zfp672_m_ex2_50ss This study CTGCATACATCTCACATTACCTTTG zfp672_m_ex3_50ss This study GGGTGTTTGTTCTGCCATACCAATA zfp672_m_ex5_30ss This study GATCCTATGGAAGGACAGTATGTAT mTsku_2ndSE_30ss This study CTTTGCTGAAATGAAACCACAGGTC Tsku_30ss This study TCTGAGAAAGGATAGGGAACCCAAT Tsku_50ss This study ACCCCTGAGTAGAGAGAGTCACCTG Gpr30_50ss This study ACCTGAAAATTTAAAAGTACTCACG Deposited Data RNA-seq data from 9 tissues from mouse and rat Merkin et al., 2012 GEO: GSE41637 Start-seq data from murine macrophages Scruggs et al., 2015 GEO: GSE62151 H3K4me3 ChIP-seq data from mouse Yu et al., 2015 GEO: GSE59896 PolyA-seq data from five mouse tissues Derti et al., 2012 GEO: GSE30198 (Continued on next page) Cell 179, 1551–1565.e1–e6, December 12, 2019 e2 Continued REAGENT or RESOURCE SOURCE IDENTIFIER PRO-seq data from mouse and rat CD4+ T cells Danko et al., 2018 GEO: GSE93229 TT-seq data from human T cells Michel et al., 2017 GEO: GSE85201 TrIP-seq data from human cells Floor and Doudna, 2016 GEO: GSE69352 Software and Algorithms TopHat Trapnell et al., 2009 http://ccb.jhu.edu/software/tophat/index.shtml STAR Dobin et al., 2013 https://github.com/alexdobin/STAR MATS Shen et al., 2014 http://rnaseq-mats.sourceforge.net/mats3.0.8/ Cufflinks Trapnell et al., 2012 http://cole-trapnell-lab.github.io/cufflinks/ MISO Katz et al., 2010 http://genes.mit.edu/burgelab/miso/ Cuffcompare Roberts et al., 2011 http://cole-trapnell-lab.github.io/cufflinks/cuffcompare/ Bioconductor N/A https://www.bioconductor.org/ BEDTools Quinlan and Hall, 2010 https://bedtools.readthedocs.io/en/latest/ SamTools Li et al., 2009 http://samtools.sourceforge.net/ GenomicRanges Lawrence et al., 2013 https://www.bioconductor.org/packages/release/bioc/ html/GenomicRanges.html Integrative Genomics Viewer N/A https://igv.org/ R (v.3.4.2) N/A https://www.r-project.org/ LEAD CONTACT AND MATERIALS AVAILABILITY Please direct any requests for further information and resources to the Lead Contact, Christopher B. Burge (cburge@mit.edu), Department of Biology, Massachusetts Institute of Technology, Cambridge MA 02138. Materials Availability Statement All reagents generated in this study are available from the Lead Contact upon request. EXPERIMENTAL MODEL AND SUBJECT DETAILS Cell lines, cell culture and treatments NIH 3T3 and HeLa cells were grown in DMEM, with high glucose and pyruvate (GIBCO), supplemented with 10% fetal bovine serum (FBS). Mouse CAD (Cath.-a-differentiated) cells were grown in DMEM/F12 (GIBCO) supplemented with 10% FBS. N2a cells were grown in DMEM,with high glucose and pyruvate (Invitrogen), supplemented with 10% fetal bovine serum (FBS). N2a cells were differ- entiated with retinoic acid as in Fiszbein et al., 2016. For morpholino oligonucleotide (MO) treatment (Gene Tools), 20 mMof morpho- lino targeting 50 or 30 splice site orMO control was addedwith Endo-Porter (Gene Tools) followingmanufacturer’s instructions to cells plated at low confluence and left for 24 h. METHOD DETAILS RNA-seq analysis and genome builds We used the strand-specific paired-end RNA-seq data from 9 tissues frommouse and rat (3 individuals each) associated withMerkin et al. (Merkin et al., 2012), available at NCBI Gene Expession Omnibus (GEO) (accession number GSE41637). Reads weremapped to themm9 and rn4 genome builds, respectively, and processed using TopHat (Trapnell et al., 2009) and Cufflinks (Trapnell et al., 2012). Cufflinks was used to estimate transcript abundance in each library (in standard FPKM units), and these values were used for splicing estimates or summed to obtain gene expression values. Alternative splicing patterns and PSI values were validated usingMISO (Katz et al., 2010). Exons were defined as in Merkin et al. (Merkin et al., 2012), requiring FPKM R 2 and meeting splice site junction read requirements implicit in the TopHat mapping. Exons with 0.05 < PSI < 0.97 in at least one tissue and two individuals were categorized as skipped exons (SE). Exons with PSI > 0.97 in all expressed tissues were defined as constitutive exons (CE), if the gene was ex- pressed in at least three tissues and two individuals. Genomic and splicing ages were defined as previously described (Merkin et al., 2015) by the pattern of species with genomic regions aligned to the exon or with an expressed exon in the orthologous gene over- lapping the aligned region, respectively, using the principle of evolutionary parsimony. As in Merkin et al. (Merkin et al., 2015), orthol- ogous exons were identified by finding annotated exons that overlapped with the query exonic region in Ensembl Pecan 19 amniota e3 Cell 179, 1551–1565.e1–e6, December 12, 2019 genome alignments (Paten et al., 2008). Exon groups with multiple overlapping exons in any species were excluded. Exons were considered ‘‘lost’’ in a species if there was no syntenic region in that species or if no exon overlapping the syntenic region was iden- tified and spliced into transcripts identified herein with a PSI R 0. Open reading frames (ORFs) were annotated as described previ- ously (Merkin et al., 2012) and used to classify exons as located in the 50 UTR, 30 UTR or coding region. CRISPR sgRNA design, genetic deletions and genotyping CRISPR-Cas cell lines with the 50 splice site of Stoml1 deleted were generated using the protocol described by Ran and coworkers (Ran et al., 2013). The single-guide RNA was designed in silico to target the 50 splice site using the CRISPR Design Tool (http://tools. genome-engineering.org) and cloned into a Cas9 expression plasmid (pSpCas9). After transfecting CAD cells with the plasmid ex- pressing Cas9 and the appropriate sgRNA, clonal cell lines were isolated and insertion/deletion mutations were detected by the Surveyor nuclease assay (IDT). Positive clones detected were amplified by PCR, subcloned into TOPO-TA plasmids, and individual colonies were sequenced to reveal the clonal genotype. RNA Extraction, RT-PCR and qPCR Total RNA was extracted using the RNA-easy kit (QIAGEN) according to the manufacturer’s protocol. Reverse transcription using M-MLV reverse transcriptase (Invitrogen) and random primers was performed according to the manufacturer’s instructions. For nascent RNA extraction, RNA was metabolically labeled with 5-Ethynil Uridine for 10 min using Click-iT (Invitrogen) and labeled RNA was extracted and amplified according to the manufacturer’s instructions. Quantitative PCR analyses were performed with SYBR green labeling using a LightCycler 480 II (Roche). ChIP and antibodies Chromatin immunoprecipitation was performed using the MAGnify Chromatin Immunoprecipitation System (Invitrogen) according to the manufacturer’s recommendations. For each immunoprecipitation, we used 10 mg of H3K4me3 antibody (PA5-17420 from Invi- trogen), 10 mg of RNA polymerase II (CTD repeat YSPTSPS) antibody (Ab817 from Abcam), 10 mg of Transcription Factor IIF1 (TFIIF- alpha) antibody (PA5-30050 from Invitrogen) and 10 mg of Rabbit IgG antibody (Invitrogen) as a negative control. DNA was purified and quantitative PCR analysis was performedwith SYBR green labeling using a LightCycler 480 II (Roche). Immunoprecipitated chro- matin was normalized to input chromatin and control IgG antibody. 50 RACE 50 RACE experiments were performed with 50 RACE System for Rapid Amplification of cDNA Ends (Invitrogen) using three gene-spe- cific primers (GSP) that anneal to the known region and an adaptor primer that targets the 50 end. Products generated by 50 RACE were subcloned into TOPO-TA vectors and individual colonies were sequenced. Plasmids, RNAi and luciferase activity assay Rat Tsku genomic region and mutants were cloned into the psiCHECK backbone. HNRNPU full-length and mutant were cloned intro the RG6 plasmid (pcDNA3.1 backbone). For transfection assays, 1 mg plasmid was transfected into each well of a 6-well culture plate using Lipofectamine 2000 (Life Technologies) according to themanufacturer’s recommendations and cells were harvested after 24 h. For knock-down experiments, a siRNA targeting the 30UTR of HNRNPU or a scrambled siRNA was transfected together with either the control (empty) RG6 plasmid or with one of the HNRNPU-expressing constructs. To measure luciferase activity, we used the Dual-Luciferase Reporter Assay System (Promega). PRO-seq data analysis PRO-seq reads in mouse and rat CD4+ T cells were mapped as in Danko et al. (Danko et al., 2018) counting reads in the interval between 500 bp downstream of the annotated TSS and whichever was shorter: either the end of the gene or 60,000 bp into the gene body. Following the analysis in Danko et al., reads were transferred to the hg19 coordinates to be compared between mouse and rat using liftOver. For each gene PRO-seq reads were defined by the sum of read counts within the gene in the interval described above. The number of readsmapping to a gene (r) were then divided by the number of reads in the library (L). RPM values were calcu- lated for each gene as r/L x 1,000,000 and divided by gene length in bp (and multiply by 1000) to get RPKM values for both species. Motif enrichment analysis The number of binding motifs for each splicing factor was calculated using RBPmap (Paz et al., 2014) by mapping each binding motif to the query sequence.We used the 94 RNA binding proteins present in the RBPmap database and added 30 additional RNA binding proteins whose binding motifs were identified by RNA Bind-n-Seq (Dominguez et al., 2018). The whole sequence of the novel exons and 20bp into the upstream and downstream introns was taken for the analysis. The enrichment of splicing factors binding motifs in mouse novel exons for each protein was calculated by dividing the mean number of binding motifs in new exons with a correlation above 0.3 with the nearby TSS by the mean number of binding moths in new exons with a correlation below 0.3 with the analog TSS. Cell 179, 1551–1565.e1–e6, December 12, 2019 e4 TT-seq data analysis TT-seq reads in human T cells were taken from Michel et al. (Michel et al., 2017) and fold change of nascent gene expression was calculated between 0 and 15 min of T cell activation. Genes expressed in both samples were then assigned to EMATS or control genes depending of their genomic structure. For splicing analyses, total RNA-seq samples from Michel et al. (Michel et al., 2017) were mapped to the hg19 genome build using STAR (Dobin et al., 2013) and splicing fold changes were processed using MATS (Shen et al., 2014). After filtering for an FDR < 0.1, we obtained 9,379 significantly changing internal exons between 0 and 15 min after T cell stimulation. Distribution of fold change in TT-seq reads after 15min of T cell stimulation were assessed for all genes, genes with significant decreased SE inclusion and EMATS genes with significant decreased SE inclusion. Polysome profile analysis Cytoplasmic, monosomal, and polysomal samples from Floor et al. (Floor and Doudna, 2016) were mapped using STAR aligner (Do- bin et al., 2013) with standard ENCODE specifications, and with the requirement that each read map uniquely to the genome. Perfectly mapping reads were then assigned to EMATS or control alternative first exons (AFEs) if they overlapped at least 25 bases with the AFE. Each instance of an overlapping read was then tallied, and AFEs with at least 5 reads assigned were considered for further analysis. AFEs were then filtered for gene representation in both EMATS and control sets to preserve a gene-matched anal- ysis, ultimately including 177 EMATS AFEs and 313 control AFEs in the analysis. Normalized read counts were used to create a ‘‘translational efficiency’’ (TE) score for each AFE by dividing the normalized read counts in a given sample by the reads counts in the cytoplasm for that AFE. Results were assessed by Wilcoxon Rank-Sum. QUANTIFICATION AND STATISTICAL ANALYSIS Definition of species-specific exons Evolutionarily new exons were identified as in Merkin et al. (Merkin et al., 2015). Genomic mappings of mouse and rat RNA-seq data were combined with whole-genome alignments to classify the species distribution of exons. Only internal exons were considered in this analysis, excluding first and last exons, and only unique exons were considered, excluding exons that arose from intra-genic duplications to avoid issues related to possibly inaccurate genome assemblies, annotations or read mappings. In all, 1,089 mouse exons were classified as mouse-specific exons and 1,571 rat exons were classified as rat-specific exons, as they were detected in RNA-seq data from mouse or rat, respectively, but not from any other species analyzed (Tables S1 and S2). Most genes that con- tained a new exon had only one, with 159 mouse genes and 276 rat genes containing more than one new exon. Transcription start site annotation TSSs in the same RNA-seq data used to classify new exons, were identified using data from Merkin et al. (Merkin et al., 2012) (GEO accession number GSE41637) mapped with TopHat combined with Ensembl annotations. As in Merkin et al. (Merkin et al., 2012), Cufflinks version 1.0.2 was used to identify novel transcripts. The set of TSSs from each library identified from transcripts as the start site of the first exon were combined with the existing Ensembl annotations and merged into a single set of annotations using Cuff- compare (Roberts et al., 2011). Cufflinks was then applied to each library to quantitate the same set of transcripts. TSS FPKM was calculated by summing the FPKM of transcripts that used the TSS. The TSS FPKM was then divided by the sum of FPKM of tran- scripts that used any other TSS to calculate the relative TSS usage. Thus, relative TSS usage was calculated dividing the FPKM of transcripts that used the TSS by that of transcripts that used a different TSS. Extensive data has accumulated that relative TSS usage derived from RNA-seq data correlate with methods that assess 50 ends of nascent RNA. Expression in FPKM from different TSSs estimated by Cufflinks from RNA-seq data strongly correlates with those derived from TT-seq analysis of nascent RNA (r = 0.76). In this manuscript, TSSs in mouse were also identified using Start-seq data from Scruggs and coworkers (Scruggs et al., 2015) downloaded from GEO (accession number GSE62151); Start-seq uses high-throughput sequencing of nascent capped RNA species from the 50 end, allowing for definition of TSSs at nucleotide resolution. TSSs were defined in 2,000 bp search windows centered on RefSeq-annotated TSSs, using the location towhich the largest number of Start-RNA reads aligned. Very closely spaced TSSs separated by less than 50 bp were considered as a single TSS in Figure 1D. The number of TSSs was also estimated by the number of H3K4me3 peaks assigned to each gene with ChIP data from Yu et al. (Yu et al., 2015) (GEO accession nos. GSE59896 and GSE59998). Software for data analysis, graphical plots and statistical analyses For data analysis we used R Bioconductor, BEDTools, SamTools, GenomicRanges, the Integrative Genomics Viewer, MISO, Cuf- flinks, STAR and MATS. All statistical analyses were performed in R (v.3.4.2) and graphical plots were made using the R package ggplot2. Lower and upper hinges of boxplots correspond to the 25th and 75th percentiles, respectively. The upper and lower whiskers extend from the hinge to the largest and lowest value no further than 1.5 3 IQR (interquartile range), respectively. Notches give approximate 95% confidence interval for comparing the medians. Statistical significance of one-way ANOVA, Tukey post hoc test, is indicated by asterisks (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, *****p < 0.00001), unless otherwise indicated. e5 Cell 179, 1551–1565.e1–e6, December 12, 2019 DATA AND CODE AVAILABILITY Data availability Data of evolutionarily new exons is available in Merkin et al. (Merkin et al., 2015) as well as here in Tables S1 and S2. The RNA-seq data from 9 tissues from mouse and rat associated with Merkin et al. (Merkin et al., 2012) is available at GEO (accession number GSE41637). The Start-seq data from Scruggs et al. (Scruggs et al., 2015) is available at GEO (accession number GSE62151), as well as the H3K4me3 data from Yu et al. (Yu et al., 2015) (accession number GSE59896 and GSE59998). PolyA-seq data from five mouse tissues is available in Derti et al. (Derti et al., 2012) (accession number GSE30198). PRO-seq data from mouse and rat CD4+ T cells from Danko et al. (Danko et al., 2018) is available at GEO (accession number GSE93229). TT-seq data in human T cells from Michel et al. (Michel et al., 2017) is available at GEO (accession number GSE85201). Polysomes sequencing (TrIP- seq) data in human cells from Floor et al. (Floor and Doudna, 2016) is available at GEO (accession number GSE69352). ADDITIONAL RESOURCES New exon inclusion, TSS usage, and species-specific expression We considered genes with new exons as all genes with a new exon with PSI > 0.05 in any of the 9 tissues sequenced. We grouped genes as control genes with no new exons and genes with new exons divided by whether the exon was included or excluded in a given tissue. We calculated the number of TSSs used in each gene in each tissue and considered genes that gained TSSs in mouse, genes that gained TSSs in rat, and genes with same number of TSSs in both species based on the numbers of TSSs for each species in each gene in each tissue, or when considering all tissues together. Gene expression was calculated by estimating transcript abun- dance with Cufflink and summing standard FPKM units per gene. The new exons were included in the length normalization for spe- cies with the exons. The FPKM normalization was done by Cufflinks and was species-specific and isoform-specific. Each tissue was run individually and transcript expression was length normalized before combining. Gene expression in mouse was compared to that in rat by taking the ratio of expression in mouse to expression of the orthologous gene in the analogous tissue in rat. Definition of new exon-proximal cleavage and polyadenylation sites Polyadenylation sites were identified using available polyA-seq data from fivemouse tissues (brain, liver, kidney, muscle, testis) (Derti et al., 2012). Only reads aligning to unique loci were retained and ends of reads within 25 nt of each other on the same strand were clustered. Polyadenylation sites were considered to be new exon-proximal cleavage and polyadenylation (nePCPA) sites if they were located within 2 kb upstream or downstream of a new exon, and as skipped exon-proximal cleavage and polyadenylation (sePCPA) sites if they were located within 2 kb upstream or downstream of skipped exons. Effects on nascent and steady state RNA levels Effects on transcription initiation should be reflected in nascent RNA, while effects on RNA stability would only be visible in steady state mRNA. In the Tsku gene, nascent RNA levels were reduced to a similar extent as steady state mRNA (Figure 2d, Extended data Figure 3b, Extended data Figure 5a-d), in both sense and antisense orientations. For other genes studied here,Stoml1 andGper1, we also observed similar effects on nascent RNA in sense and antisense directions (Figure 2c, Extended data Figure 3b, Extended data Figure 4a-c). Furthermore, the model invoking inhibition of PCPA involves U1 snRNP binding at a 50 splice site, but we observed increased gene expression from creation of a 30 splice site. Thus, our observations are consistent with splicing-dependent regulation of transcription initiation but not with models involving PCPA. Cell 179, 1551–1565.e1–e6, December 12, 2019 e6 Supplemental Figures A 1.00 B C genes with no 40 new exons all genes 0.3 genes in tissues 0.75 new exon ψ < 0.05 genes with mouse-specific new exons all genes 30 genes in tissues genes with mouse-specific new exon ψ > 0.05 new exons 0.2 0.50 20 0.1 0.25 10 p-value = 1.53e-06 0.00 0.00 5 10 −2 0 2 5 10 15 no. of TSSs in brain (FPKM mouse / FPKM rat), log no. of TSS (RNA-seq)2 after match it D E F 0.6 genes w/ mouse TSS all genes < rat TSS 30 all genes genes with rat-specific new exons 0.4 genes with genes w/ mouse-specific new exons mouse TSS 20 = rat TSS 0.2 genes w/ 10 mouse TSSs > rat TSSs 0.0 2 4 6 8 0 0 10 20 30 40 50 no. of H3K4me3 peaks 5 10 15 proportion of genes (%) no. of TSS (RNA-seq) all genes genes with rat-specific new exons G H I J **** 3 * **** **** 1.50 5.0 **** 0.75 1.25 2.5 2 0.50 0.0 1.00 1 −2.5 0.25 0.75 −5.0 0 0.00 all genes with 0.50 ψ < 0.05 < ψ ψ > no. of no. of genes rat-specific mouse rat 0.05 < 0.5 0.5 TSS = 1 TSS > 1 lost exons Figure S1. Genes with Evolutionarily New Exons Have Increased Numbers of TSSs across Species, Related to Figure 1 (A) Fold change in gene expression between mouse and rat for mouse control genes with no evolutionarily new exons (black, dotted line), genes with mouse- specific new exons in tissues where inclusion of the new exon is not detected, PSI < 0.05 (gray), and genes with new mouse-specific exons in tissues were the exon is included, PSI > 0.05 (pink). Statistical significance by Mann-Whitney U test is indicated between genes with mouse-specific new exons in tissues with PSI < 0.05 and tissues with PSI > 0.05. (B) Distribution of the number of TSSs per gene using RNA-seq data acrossmultiple species andmultiple tissues, for all genes expressed inmouse and geneswith mouse-specific new exons. Distributions are significantly different by Kolmogorov-Smirnov test. (C) Distribution of the number of TSSs per gene in the mouse brain using RNA-seq data, for all genes expressed in mouse and genes with mouse-specific new exons, after matching the distribution of gene expression levels between the two groups using the MatchIt package in R. Distributions remain significantly different by Kolmogorov-Smirnov test aftermatching the gene expression levels between the groups, demonstrating that, independent of gene expression, genes with mouse-specific new exons are enriched in multiple TSSs. (D) Distribution of the number of -see3 peaks per gene using H3K4me3ChIP-seq data for all genes expressed inmouse (gray) and geneswithmouse-specific new exons (dark red). Distributions are significantly different by Kolmogorov-Smirnov test. Genes with mouse-specific new exons are enriched in H3K4me3 peaks. (E) Distribution of the number of TSSs per gene for all genes expressed in rat (gray) and genes with rat-specific new exons (green). Genes with rat-specific new exons are enriched in multiple TSSs (by Kolmogorov-Smirnov test). (F) The proportion of genes with fewer TSSs in mouse (genes w/ mouse TSS < rat TSS), genes with the same number of TSSs in both species (genes we mouse TSSs = rat TSS), and genes that have more TSSs in mouse, for all genes expressed in both species (gray) and for genes with rat-specific new exons (green). Statistical significance indicated by asterisks corresponds to one-way ANOVA, Tukey post hoc test (NS = not significant). (G) Fold change in the number of TSSs used per gene between tissueswheremouse-specific exons are included (PSI > 0.05) and excluded (PSI < 0.05), formouse genes and for the same tissues in rat. Evolutionary gain of internal exons and of TSSs are associated, but only in those tissues where new exons are included. (legend continued on next page) no. TSSs tissues ψ > 0.05 / no. TSSs tissues ψ < 0.05 proportion of genes (%) cumulative fraction of genes no. TSSs mouse / no. TSSs rat proportion of genes (%) proportion of genes (%) ψ value of new exon (FPKM rat / FPKM mouse), log2 proportion of genes (%) *** *** NS (H) Distribution of PSI values of new exons binned by the number of TSSs used in the same gene, for 9 tissues pooled together in mouse. (I) Fold change in the number of TSSs used per gene betweenmouse and rat for 9 tissues, for genes with mouse-specific new exons, binned by c value of the new exon in each tissue. (J) Fold change in gene expression between rat and mouse, for all genes expressed in rat (gray) and genes with rat-specific lost exons (blue). A *** *** B genes that gained 488 0.1 0.1 new exons 531 0.0 0.0 −0.1 −0.1 4729 genes genes w/ genes w/ genes w/ genes w/ mouse TSS mouse TSS mouse TSS mouse TSS that gained ≤ rat TSS ≥ rat TSS ≤ rat TSS ≥ rat TSS new TSS before MatchIt after MatchIt C D -2 -1 1 2 3e-5 1.0 ** 0.5 2e-5 0.0 1e-5 −0.5 −1.0 -2e+5 -1e+5 0 1e+5 2e+5 <−3 −2 −1 1 2 >3 position of TSSs relative to new exon (bp) order of TSSs relative to new exon Figure S2. Splicing of New Exons Is Associated with TSS Gain and Relative TSS Usage, Related to Figure 2 (A) Distribution of the fold change in gene expression levels between mouse and rat for genes with fewer or same number of TSSs used in mouse than rat (white) and genes with more TSSs used in mouse than rat (brown), before (left panel) and after (right panel) balancing the distribution of gene expression levels in mouse between the groups byMatchIt. Evolutionarily change in gene expression remain significantly different when balancing gene expression levels in mouse between groups, demonstrating that, independently of gene expression levels in one species, genes gaining TSSs in mouse have increased gene expression levels compared to rat. (B) Venn diagram showing the overlap between genes that gained mouse-specific new exons and genes that gained new TSSs in mouse. (C) TSS position relative to the start coordinate of the new exon in genes with mouse-specific new exons, for all TSSs used in 9 tissues in mouse. (D) Spearman correlations between the usage of a particular TSS and the PSI value of the new exon acrossmultiple tissues for all TSSs used in genes withmouse- specific new exons, binned by their relative position to the new exon with negative numbers for TSSs located upstream of the new exon and positive numbers for TSSs located downstream of the new exon. TSS density (FPKM mouse / FPKM rat), log (FPKM mouse / FPKM rat), log ρSpearman (new exon ψ versus rel. TSS use) A B RNAPII 1.00 1.0 10 min into cells 0.75 * 5-ethynyl uridine (EU) 0.9labeled nascent RNA 0.50 ** 0.5 *** 0.6 ** 0.25 streptavidin 0.3 magnetic beads *** 0.0 0.00 0.0 *** MO MO MO MO MO MO 3’ss control 5’ss control 3’ss 5’ss and 5’ss WT 5’ss 5’ss mut 1 mut 2 C TSS –1 TSS +1 D TSS +2 20 A B C D E F G H I J K L M H3K4me3 15 10 10 WT 8 5’ss mutant 5 6 0 4 −5 2 −2e+05 0 2e+05 position relative to new exon (bp) 0 A B C D E F G H I J K L M E genomic position 10 RNAPII 4 WT 5’ss mutant 3 5 2 1 0 < − −1 −20 100 200 >10 0 0 0 10 0 , ,00 0 0 t to o 2 to ,00 0 00 to 0 – 0− 100 0 10, 00 M 2 0 0 0A B C D E F G H I J K L 00 genomic position position relative to new exon (bp) F G 50 0.75 tissues w/ ρSpearman = 0.094 new exon 40 ψ < 0.05 p-value = NS tissues w/ 30 0.50 new exon ψ > 0.05 20 0.25 10 0 0.00 −4 −2 0 2 4 6 0 1 2 3 4 5 6 7 log(gene expression in mouse / gene expression in rat) no. of nePCPA sites for the homologues gene in the same tissue (legend on next page) density Fold enrichment (IP - input rel to IgG) Fold enrichment (IP - input rel to IgG) Expression relative to control Expression relative to control no. of nePCPA sites PRO-seq reads PRO-seq RPM in mouse - rat in mouse genes with new exons in genes with mouse-specific exons biotin fold change (INC/EXC) Figure S3. Splicing of New Exons Impacts Transcription from Nearby TSSs, Related to Figure 3 (A) A diagram representing the technique used to label nascent RNAwith 5-ethynyl uridine and pull down the nascent RNAwith the click-it method. Fold change in RNA levels ofGpr30 (left) and Tsku (right) in NIH 3T3 cells measured by qPCR. NIH 3T3 cells were transfected with 20 mMmorpholino (MO) targeting the 50 splice site of the new exon in Gpr30 or 20 mMMO targeting the 30 and/or the 50 splice sites of the new exon in Tsku for 24 h. Mean ± SEM of displayed distributions, n = 3 biological replicates. Statistical significance indicated by asterisks corresponds to one-way ANOVA, Tukey post hoc test. (B) Relative exon inclusion / exon exclusion of the mouse-specific new exon in Stoml1 gene is shown, measured by qPCR of nascent RNA in wild-type CAD cells and cells with CRISPR/cas-mediated mutations in the 50 splice site of the new exon. Mean ± SEM is shown for n = 3 biological replicates. Statistical significance indicated by asterisks corresponds to one-way ANOVA, Tukey post hoc test. (C) H3K4me3 and RNAPII (CTD repeat YSPTSPS) profiles in Stoml1 gene in CAD cells determined by ChIP assay followed by qPCR with the regions indicated in the top panel. Values of two independent immunoprecipitations normalized to input and the mean value for control IgG antibody are shown for each region. Wild- type cells (gray) and cells with CRISPR/cas-mediated mutations in the 50 splice site of the new exon (blue) are shown. (D) Fold change in nascent RNA expression (PRMunits) betweenmouse and rat relative to the start coordinate of the new exon for geneswithmouse-specific new exons in CD4+ T cells using PRO-seq data. (E) Distribution of nascent RNA expression levels in mouse genes with new exons binned by the position relative to the start coordinate of the new exon in CD4+ T cells using PRO-seq data. (F) Distribution of the number of polyadenylation sites used 2 kb upstream/downstream of new exons per gene in tissues new exon is excluded (PSI < 0.05, gray) and tissues with inclusion of new exons (PSI > 0.05, pink) for all genes with new exons. Distributions are not significantly different by Kolmogorov-Smirnov test. (G) Scatterplot showing the relationship between the number of nePCPA sites and the fold change in gene expression levels between mouse and rat. These variables are not significantly associated by Spearman correlation test. Polyadenylation sites for 5 tissues in mouse were analyzed using polyA-seq data (Derti et al., 2012). A E F B G C H D I Figure S4. Creation of a New Exon Can Activate Nearby TSSs, Related to Figure 4 (A) Isoform expression for Tsku gene in NIH 3T3 cells, measured by RT-PCR with primers targeting TSS –2 (top), TSS –1 (bottom) and exon E4.CE. Cells were transfected for 24 h with 20 mM MO targeting the 30 and/or 50 splice sites of the new exon as indicated. (legend continued on next page) (B) Fold change in inclusion and exclusion levels of the mouse-specific new exon in Tsku gene and antisense transcription levels from both TSSs measured by isoform-specific qPCR of nascent RNA with primers illustrated by arrows. Exon exclusion levels are measured from both alternative first exons to the following skipped exon downstream the mouse-specific new exon. NIH 3T3 cells were transfected with control MO or MO targeting the 30 and/or 50 splice sites of the new exon. A decrease in the inclusion level of transcripts starting at TSS –2 is compensated by an increase in the exclusion levels, while total level of transcripts starting at TSS –1 is reduced by MO treatment. Mean ± SEM of displayed distributions, n = 3 biological replicates. Statistical significance indicated by asterisks corresponds to one-way ANOVA, Tukey post hoc test. (C) RNAPII profile in Tsku gene in NIH 3T3 cells determined by ChIP assay followed by qPCR with the regions indicated in panel A. Values of two independent immunoprecipitations normalized to input and the mean value for control IgG antibody are shown for each region. NIH 3T3 cells transfected for 24 hours with 20 mM control MO or MO targeting both 30 and 50 splice sites of the new exon. (D) Alignments and identity between mouse (mm) and rat (rn) of the DNA sequence of TSS –1 and mouse-specific new exon in the Tsku gene. (E) Splicing patterns of the Tsku gene in HeLa cells transfected with the hybrid constructs shown in Figure 4d. The creation of the mouse 30 splice site (rn + mm 30ss) or of a stronger 30 splice site (rn + strong 30ss) of the mouse-specific new exon in the rat sequence promotes the inclusion of the mouse-specific new exon in the rat context only when the wild-type 50 splice site is maintained (but not in the mm 30ss + mut 50ss construct). (F) Sequence of the 50 end of Tsku transcripts generated by 50 RACE in HeLa cells transfected with rat Tsku constructs with the 30 splice site of the mouse-specific new exon (50 RACE clone A, clone B) aligned to the mouse sequence (mm) and the rat sequence (rn). For 80% of the sequenced transcripts, the 50 end mapped 1 bp upstream of the position of mouse TSS –1 (clone A), while in the remainder the 50 end mapped 19 bp upstream of (clone B). (G) Ratio of mean density of RBPmap binding motifs (Paz et al., 2014) in mouse-specific new exons with mild/high correlation with usage of proximal upstream TSS (Spearman correlation > 0.3) to density of motifs in mouse-specific new exons with low correlation (spearman correlation < 0.3) for selected splicing proteins (top ten splicing proteins by ratio are shown). (H) Z-scores of splicing factors binding sites relative to the position of the new exon in mouse Tsku. (I) Isoform-specific qRT-PCR analysis of fold change in exclusion levels of the mouse-specific new exon in Tsku gene from both TSSs following treatment with a control siRNA or an siRNA targeting HNRNPU (siHNRNPU) and rescues in NIH 3T3 cells. Mean ± SEM of displayed distributions, n = 3 biological replicates. Statistical significance indicated by asterisks corresponds to one-way ANOVA, Tukey post hoc test. A B 10 * C6 9 6 skipped 4 8 exons 7 new 2 skipped new exonsexons 2 exons 6 5 < - -10 -10 1 0 0.5 1 0 10 > 0 0.5 1 10 00 0 0 0 t 0 t 00 10 00 o o to - -1 1 to 001 0 1 00 0 00 0 00 0 distribution of rel. TSS use values distribution of ψ values for 0 0 for TSSs in mouse exons in mouse position of TSSs relative to new exon (bp) D **** E **** F ** 1.0 10 0.4 0.5 5 0.0 0 0.2 −0.5 −5 −1.0 0.0 rel. TSS rel TSS rel. TSS rel. TSS rel. TSS rel. TSS rel. TSS rel. TSS use 1st use 2nd use 3rd use 4th use 1st use 2nd use 3rd use 4th quartile quartile quartile quartile quartile quartile quartile quartile < −2 –5 – 1 5−2 50 0 10 00 10 00 > 50 0 0 0t 0 t 0 0 20 5 0 0 o o t 0 0 t to – o 1 to o 0 G 1 – 00 5 25 0 400 nt 5000 nt – 05 00 10 0 0 00 0 00 0 0 00 position of TSSs relative to skipped exon (bp) TSS –2 TSS –1 E2.new E3.SE E4.CEexon 1.00 1.00 0.9 1.0 0.75 0.75 * 0.6 0.50 0.50 0.5 0.3 0.25 0.25 0.00 0.00 0.0 0.0 MO MO MO MO MO MO MO MO control 5’ss control 5’ss control 5’ss control 5’ss E3.SE E4.CE TSS –2 E4.CE TSS –1 E4.CE E4.CE Figure S5. Highly Included Alternative Exons Activate Transcription from Weak Promoters, Related to Figure 5 (A) Distribution of the PSI values of mouse-specific new exons (blue) and SE (gray) in mouse across 9 tissues. (B) Distribution of 50 splice site scores of mouse-specific new exons, binned by the relative position to the next upstream TSS used in the same gene. 50 splice site scores were calculated using MaxEntScan (Yeo and Burge, 2004). (C) Distribution of the PSI values of first exons associated with TSSs in genes with mouse-specific new exons (blue) and in genes with SEs in mouse (gray). (D and E) Difference in TSS usage based on relative TSS use value (D) and FPKM (E) in tissues with high versus low inclusion of skipped exons (SE), in the same gene across multiple tissues for proximal and upstream TSSs (within 1kb upstream the SE) used in genes with SEs in mouse, binned by quartiles of PSI values of the TSSs. F, Difference between relative TSS usage values in tissueswith high versus low inclusion of skipped exons (SE), for all weak TSSs (bottom quartile) used in geneswith skipped exons inmouse, binned by their position relative to the SE. G, Fold change in inclusion (far left), exclusion levels (left center, right center) and total levels (far right) of the skipped exon in Tsku gene (E3.SE) from both TSSs measured by isoform-specific qPCRs with primers shown by arrows. Exclusion levels were measured from both alternative first exons relative to the next constitutive exon downstream of the skipped exon. NIH 3T3 cells were transfected with control MO or MO targeting the 50 splice site of E3.SE. Mean ± SEM of displayed distributions, n = 3 biological replicates. Statistical significance indicated by asterisks corresponds to one-way ANOVA, Tukey post hoc test. relative expression rel. TSS use tissues w/ SE ψ > 0.05 – density rel. TSS use tissues w/ SE ψ < 0.05 relative expression relative expression TSS FPKM tissues w/ SE ψ > 0.2 – TSS FPKM tissues w/ SE ψ < 0.05 5’ss score of new exon (bits) relative expression rel. TSS use tissues w/ SE ψ > 0.05 – rel. TSS use tissues w/ SE ψ < 0.05 density A B SRSF5NONOSRSF1 PUF60 MATR3 HNRNPF TIAL1 PTBP1 HNRNPU NCBP2 TFIP11 6 TARDBPTIA1 HNRNPL U2AF2 SMN1 SF1 PRPF8 GEMIN5 PCBP1 SRSF9 PSIP1 HNRNPK HNRNPLL SFPQ 4 RBM22RBM15 U2AF1 TRA2A RBM17 AKAP8L EFTUD2 KHDRBS1 HNRNPC PPP1R8 GPKOW HNRNPA1 SRSF3 RBM25 QKI 2 MBNL1 HNRNPM RAVER1 FMR1 SND1 RBM39 PPIG SUGP2 EWSR1 DAZAP1 SRSF4 PCBP2 CELF1 RBM47 0 RBFOX2HNRNPA2B1 SRSF7 0 2000 4000 KHSRP FUS no. of TSS with expression change ADARHNRNPD ZRANB2 following RBP knock-down PRPF6STAU1 BUD13 DDX5 CCAR2 C 0 100 200 300 no. of genes with changes in promoter usage following splicing factor knock-down seed of 10 splicing factors ribonuclease activity double/single stranded RNA binding RNA binding misfolding protein binding pre-mRNA binding / splicing DNA binding protein binding involved in 3’ processing gene expression regulators translation regulation stress response RNA interference transcription from RNAPII promoter transcription co-regulation mithocondrial fission telomerase regulation D E F reterograde neuronal dense core vesicle transport * TSS –2 TSS –1 anterograde neuronal dense core vesicle transport * E9.SE 1.3 kb 1.5 kb calcium ion homeostasis * regulation of receptor localization to synapse * NLS NLS EHMT2 membrane raft organization * 2 kb BMF histone lysine demethylation * 6 * 1.0 regulation of dendrite spine development ** 1.0 leukocyte degranulation * androgen receptor signaling pathway * 4 * * neuron maturation * * * IRE1-mediated unfolded response * 0.50.5 regulation of interferon type I * 2 **** regulation of dendrite morphogenesis *** * **** neutrophil activation ** regulation of neuron apoptosis ** 0 0.0 axon guidance * Undif Dif Dif MO MO MO MO 0.0 control E9.SE E9.SE E9.SE TSS –2 TSS –1 SE GE 0.0 2.5 5.0 7.5 10.0 12.5 fold enrichment of GO category Control siPTBP1 TSS –2 TSS –1 TSS –2 TSS –1 total Figure S6. Splicing Regulators Impact TSS Usage, Related to Figure 6 (A) Histogram and smoothed density of number of TSSs with significant expression change following depletion of each of 250 RNA binding protein genes. Mean two cell lines (HepG2 and K562) is plotted for each RBP. (B) Number of genes with significant difference in promoter usage associated with depletion of 67 splicing factors (SF). The red line indicates the cutoff for the top ten splicing factors driving the largest changes in promoter usage, while the green line indicates the cutoff for bottom ten control splicing factors driving the fewest changes in promoter usage. (C) Protein interaction network for 10 control splicing factors driving the fewest changes in promoter usage. The control 10 splicing factors in red primarily interact with 88 other proteins, generating a network with 98 nodes and 410 edges, a diameter of 3, an average weighted degree of 4.29, an average clustering coefficient of 0.43 and an average path length of 1.32. Nodes represent proteins and links represent the interactions among them. Node size and label size is proportional to the protein connectivity (number of interactions a protein establishes with others). Protein interaction data were collected from STRING (Szklarczyk et al., 2015) (legend continued on next page) density (x 10–4) SF knocked down relative expression relative expression relative level to control (RNA-seq) and networks were built using Gephi (http://gephi.org/). Colors represent the Gene Ontology analysis of 88 proteins included in the interaction network, excluding the 10 seed control splicing factors. (D) Gene Ontology analysis of 1414 human genes with the strongest EMATS potential. Fold enrichments shown for the most significant categories with asterisk indicating adjusted p-values and color indicating relation to neuron development. (E) Exon-intron organization of EHMT2 gene. qRT-PCR analysis of expression of TSS-2 and TSS-1 in (left) N2a differentiated cells and (right) cells treated withMO targeting the SE (E9.SE) normalized to expression of housekeeping genes Hprt and GAPDH. Mean ± SEM of displayed distributions for n = 3 independent experiments. In Fiszbein et al., 2016 the authors only worked with Ehmt2 isoforms expressed from the gene-distal promoter that contain only the one constitutive NLS near the internal alternative exon. Here, we also analyze the expression of isoforms starting from the gene-proximal promoter that contains another NLS in the alternative first exon. (F) Exon-intron organization of human BMF gene. RNA-seq analysis of expression of BMF in HepG2 cells following PTBP1 knockdown normalized to expression of control cells. Inclusion levels of the skipped exon, as well as levels of relative usage of the alternative TSSs (TSS –2, TSS –1) and total gene expression are shown. Mean ± SEM of displayed distributions for n = 2 replicates. 6 4 EMATS TSS control TSS 2 0 S ly2 ly3 4 5 ly 6 ly7 y8 80 po po l poly poly po po po Figure S7. EMATS-Regulated Isoforms Are More Highly Translated, Related to Figure 7 EMATS TSSs (pink) have increased translation efficiency (TE) relative to gene-matched control TSSs (gray). For each TSS, the TE was calculated for each polysomal fraction. Boxplots show the distributions of TEs. translation efficiency TrIP-seq reads (polysome / total cytoplasmic)