5’ end sequence contributes to transcript expression

First, we asked whether the composition of the transcribed sequence affects its own transcription. To this end, the library included TSS-proximal regions from different transcript classes. As a negative control, mRNA termination regions were included, which exhibited low abundance (Fig. 1c,​,d),d), as expected. TSS-proximal regions from mRNAs supported the highest levels of RNA and protein, followed by long intergenic non-coding RNAs (lincRNAs), and then uaRNAs and enhancer RNAs (eRNAs) (Fig. 1c,​,d).d). This result is consistent with the endogenous elongation potential for these transcript classes and suggests that the composition of the TSS-proximal sequence contributes to transcription regulation. Nonetheless, despite significant differences between transcript classes, the distributions of RNA levels are overlapping, with some highly abundant eRNA sequences observed and conversely, some mRNA sequences detected at low levels. These findings are consistent with the many similarities observed for transcription of coding and non-coding RNAs²¹. Next, we tested whether functional sequence elements were limited to the most TSS-proximal regions or could be found equally in TSS-distal regions (160-333bp from the TSS) (Fig. 1e, Extended Data Fig. 1b). Interestingly, TSS-proximal regions supported higher abundance than their distal counterparts. Thus, we suggest that the mammalian genome has evolved to contain signals within the initially transcribed region that directly influence elongation potential.

We specifically evaluated sequences from typical enhancers (TE) and super enhancers (SE, ref. ²²), since it has been proposed that transcription termination is particularly important within SEs to prevent DNA damage caused by collisions between closely spaced, convergently transcribing RNAPIIs²³. In agreement with this idea, we observed that eRNA sequences from SEs are consistently more negative than those in TEs (Fig. 1f, Extended Data Fig. 1c), suggesting that eRNA sequences themselves contribute to efficient termination observed in SEs.

Effect of the transcribed sequence is independent of context

To test whether these conclusions were generalizable, we repeated the RNA-level screen after integrating the sequence library downstream of the 4930461G14Rik lincRNA TSS. We selected this locus, since this lincRNA is many kb long, highly transcribed and, in comparison to the uaRNA locus, displays higher levels of H3K27ac, a larger domain of H3K4me3 and less H3K4me1 (Fig. 1g). Screening the library at this lincRNA locus, without introduction of a fluorescent reporter, revealed a strong correlation with the data obtained at the uaRNA locus (Fig. 1h, Extended Data Fig. 1d, Supplementary Table 4). We conclude that the signals contained in transcribed sequences can be conveyed independent of chromatin context or promoter identity.

Transcribed sequence directly affects transcription levels

Since uaRNAs and eRNAs have been reported to be unstable due to rapid degradation by the exosome complex^2,24–26, we tested whether the differences in RNA abundance observed could be explained by RNA stability. First, we repeated the RNA screen with the library at the uaRNA locus after knockdown of the exosome subunit EXOSC3 (RRP40; Extended Data Fig. 2a). EXOSC3 depletion resulted in a general stabilization of transcripts from the uaRNA reporter locus (Extended Data Fig. 2b). However, using internal normalization of INSERT-seq data, we found that the relative levels of inserts were unchanged after EXOSC3 knockdown (Fig. 2a, ​,b),b), with mRNAs regions remaining the most positive and eRNAs the most negative (Fig. 2b, Supplementary Table 5). Second, we isolated nascent RNA from the library-containing pool of cells using a modified version of PRO-seq²⁷, and used this in an RNA-based screen. Consistently, we found the pattern of mRNA>lincRNA>uaRNA>eRNA in insert abundance in the nascent RNA (Fig. 2c), and a good correlation between nascent and steady state RNA levels (Extended Data Fig. 2c). Agreement was also observed when performing INSERT-seq using chromatin-associated RNA (Extended Data Fig. 2d,​,e),e), which is enriched in nascent and recently synthesized transcripts and thus less sensitive to effects of RNA stability. We conclude that the observed effects are largely independent of RNA stability and that the initially transcribed sequence directly affects transcription. Our data support that this regulation occurs predominantly at the level of transcription elongation, since INSERT-seq selectively measures productive transcript formation; however, contributions of other mechanisms cannot be excluded.

GC content inherently affects transcriptional output

mRNAs are more GC-rich than average genomic regions, especially at 5’ ends, which often overlap with CpG islands¹. The GC content of non-coding RNAs is closer to the genome average of 42% G/C in mice²⁸ (Extended Data Fig. 3a). To test whether GC content influences expression, the library contained ~1000 synthetic control sequences generated over a range of GC contents. Notably, these synthetic controls were devoid of the most common PAS hexamers to remove canonical termination sequences. We observed that inserts with higher GC content had a higher abundance in all assays (Fig. 3a, Extended Data Fig. 3b,​,c),c), indicating a striking impact of GC content on RNA production. We tested if this relation was due to enrichment of CpG dinucleotides specifically²⁹, but found CpG content was not a better predictor of RNA level than overall GC content (Pearson correlation coefficient of 0.47 in both cases, Fig. 3a and Extended Data Fig. 3d). We therefore conclude that early-transcribed regions enriched in G and C nucleotides support higher levels of processive transcription than sequences with a high AT content.

Open in a separate window

Figure 3.

GC content inherently affects transcriptional output.

a, Relation between GC content of synthetic control sequences and their steady-state RNA levels (n = 1,059). The red line is the best linear fit through the data. Pearson r = 0.47, P < 0.0001. b, Elongation indices from mESC PRO-seq data at endogenous genomic locations of uaRNAs (left) and eRNAs (right) sequences tested in INSERT-seq. Elongation Index was calculated as the ratio of PRO-seq density corresponding to elongating RNAPII (window from +50 to +250 downstream of the TSS) divided by the promoter-proximal RNAPII density (from the TSS to +50). Higher values indicate more efficient elongation. Loci were grouped by the GC content of their TSS-proximal regions (TSS+6 to +179), as included in the library, comparisons by Kruskal-Wallis test, **** indicates P < 0.0001, higher P values are indicated in the plot. c, Density plot of steady-state RNA abundance levels corrected for the predicted abundance based on the GC content (based on the fit line in panel a) by subtraction. Groups as in Fig. 1c. All three groups differ from the controls by Kruskal-Wallis test, P < 0.0001. c, Like b, but showing TSS-proximal regions from uaRNAs (n = 1,743) and eRNAs from typical enhancers (TE; n = 1,506) and super enhancers (SE; n = 600) alongside controls. Only SE eRNA regions were significantly different from the controls by Kruskal-Wallis test, P < 0.0001. d, Density plot of GC-corrected steady-state RNA levels of all TSS-proximal and TSS-distal uaRNA and eRNA regions, grouped by the presence of a canonical PAS hexamer (AWTAAA, where W=A/T). Controls n = 1,059, without PAS n = 3,895, with PAS n = 511. The groups with and without PAS hexamer are not significantly different by Mann-Whitney test (P = 0.69).

We corroborated these results using PRO-seq data from wild-type mESCs to look at transcriptionally engaged RNAPII over the sequences included in the INSERT-seq library. For uaRNA and eRNA loci included in our screen, we found that transcribed regions with higher GC content exhibited elevated levels of RNAPII (Extended Data Fig. 3e). Furthermore, with increasing GC content these loci display a higher Elongation index (downstream over upstream PRO-seq signal), indicative of more efficient transcription elongation (Fig. 3b). Thus, even in the genomic context, where numerous other factors are at play to regulate transcriptional output, we find that a high GC content in the initially transcribed region promotes processive transcription.

To probe to what extent GC content explained the behavior of genomic sequences in our library, the relation between GC content and abundance was used to “correct” the abundance of all inserts. After GC-correction, mRNA and lincRNA 5’ sequences remained positive and mRNA 3’ ends negative compared to controls (Fig. 3c). However, the RNA levels of uaRNAs and typical enhancer (TE) eRNAs were indistinguishable from controls after GC-content correction (Fig. 3d). Only eRNAs from SEs showed a small but significant shift towards lower abundance (Fig. 3d), supporting the idea that these sequences have evolved to induce more efficient transcription termination.

Notably, our data imply that early termination of eRNAs and uaRNAs does not rely on specific sequence motifs. In fact, the presence of PAS hexamers in uaRNAs and eRNAs caused no additional reduction in RNA level in our screen (Fig. 3e). We conclude that the high AT content of the genomic regions from which most uaRNAs and eRNAs arise inherently causes RNAPII to be termination prone, and they have not evolved to enrich for specific motifs. Indeed, high AT content increases RNAPII pausing and backtracking³⁰, since AT-rich RNAs fail to form stable secondary structures that prevent RNAPII backtracking^31,32 and AT-rich RNA-DNA hybrids destabilize the ternary elongation complex. Consequently, T-rich sequences promote termination by several different RNA polymerases^33–36. Thus, we propose that frequent RNAPII stalling and backtracking within AT-rich non-coding regions increases the chance of early termination, as it does downstream of a PAS at mRNA 3’ ends³⁷.

Co-transcriptionally spliced introns boost processive transcription

We then focused on mRNA TSS-proximal regions. Searches for sequence motifs that could convey positive signals in the most abundant TSS-proximal mRNA inserts identified the same top hit using MEME³⁸ and Homer³⁹: a 5’ splice site (5’SS) motif (Fig. 4a). Splicing, and the 5’SS in particular, have been suggested to stimulate RNA production^8–17, and the mechanisms underlying this activity are a subject of great interest¹⁸. Thus, we used INSERT-seq to systematically assess the effect of intron sequences on RNA production. We first measured the co-transcriptional splicing efficiency of inserts (using nascent RNA, see methods) and found that in mRNA TSS-proximal regions, efficiently spliced inserts were most abundant (Fig. 4b). Further, in inserts containing annotated mouse introns, splicing efficiency also correlated with abundance (Fig. 4c).

Open in a separate window

Figure 4.

Co-transcriptionally spliced introns boost transcription.

a, De novo motif identified by MEME and Homer algorithms to be most significantly enriched in the top-10% over the bottom-50% of TSS-proximal mRNA regions by GC-corrected steady-state RNA level. 5’ splice site (5’SS) motif from the JASPAR database (SD000.1) is shown for comparison. b, Density plot of GC-corrected steady-state RNA levels of TSS-proximal mRNA regions grouped by splicing efficiency (spliced/ total). Controls n = 1059, <3% spliced n = 3,687, 3-30% spliced n = 102, >30% spliced n = 78. To focus on co-transcriptional splicing, splicing efficiency for panels a and b was calculated using nascent RNA data. All groups differ significantly from each other by Kruskal-Wallis test, P < 0.0001. c, Steady-state RNA levels of inserts containing wild-type introns (unbarcoded) grouped by splicing efficiency measured using the nascent RNA screen data. <3% spliced n = 76, 3-30% spliced n = 107, >30% spliced n = 198, significance tested by Kruskal-Wallis test. d, Reporter level measured in clonal mESC lines with an intron from the mouse Atp5a1 gene inserted upstream of the fluorescent reporter at the uaRNA locus, compared to clonal control lines (CRISPR clones without intron integration). Effects in chromatin-associated (Chr) and steady-state (Total) RNA were determined by RT-qPCR and normalized to an internal control gene (TBP). Fluorescence was measured by flow cytometry. Comparisons were made using two-way ANOVA and corrected with the šidák method. n = 3 independently generated cell lines per genotype, though fluorescence was measured for only 2 of the control lines. e, Steady-state RNA levels of variations of a subset of introns. Introns without (open) and with (filled) a strong PAS terminator⁴⁰ replacing part of the intronic sequence were compared in otherwise wild-type introns and 5’SS mutant introns (each n = 19). For reference, strong PAS terminators were inserted in a background of scrambled (scr) introns (n = 18). Comparisons by two-sided Wilcoxon tests.

For all panels, **** indicates P < 0.0001, all higher P values are indicated in the plots.

This result was consistent across all screen read-outs (Extended Data Fig. 4, Supplementary Table 6). Moreover, clonal cell lines containing at intron at this locus displayed significantly elevated RNA and protein expression as compared to intron-less controls (Fig. 4d). Importantly, as described below, analysis of splice site mutants excludes the possibility that this stimulation is caused by intronic enhancers. We note that the locus where the sequences were integrated contains no PAS hexamers, thus the stimulation caused by spliced introns does not rely on the suppression of PAS-mediated transcript cleavage or termination.

Previous work nonetheless indicates that a 5’SS suppresses the recognition of a downstream cryptic PAS’^{6–8,15–17}, and prior study of a strong PAS-containing terminator in a plasmid context showed that PAS recognition can be inhibited inside an intron⁴⁰. To rigorously test this model, we inserted the strong 49-nt PAS-containing terminator⁴⁰ inside a variety of introns. If the terminator leads to RNA cleavage inside an insert, that will prevent amplification and lead to low abundance in the sequencing library. While addition of the terminator sequence into wild-type, splicing-competent introns had no effect on RNA levels (Fig. 4e), mutation of the 5’ splice site enabled the terminator to significantly reduce transcript abundance. In fact, RNA levels dropped to the same level as that of inserts with a terminator in a scrambled sequence background (Fig. 4e). Thus, while the PAS terminator can robustly induce termination at this locus, this property is fully masked when the terminator is within a spliced intron. To our knowledge, these data are the first to demonstrate in a genomic context that introns can broadly suppress the usage of a strong, canonical PAS terminator.

Distinguishing splicing-dependent and splicing-independent effects of 5’ splice sites

Because splicing efficiency correlated with insert abundance, we sought to establish a causal link using splice site mutants that abrogated splicing in efficiently spliced wild-type introns (see methods). Small (1 and 3 nt) mutations in either 5’SS or 3’SS greatly reduced RNA and protein levels (Fig. 5a, Extended Data Fig. 5a), indicating that it is the process of splicing and not the sequence composition inside introns that promotes transcription.

Open in a separate window

Figure 5.

Splicing-dependent and -independent role of the 5’SS.

a, Steady-state RNA levels of intron-containing inserts with wild-type (wt) or mutant (m) splice sites. Only introns are shown of which the wild-type version was >30% spliced and mutants <3% spliced in nascent RNA. 5’SS mutants n = 51, 3’SS mutants n = 23, WT n = 52, comparisons by Kruskal-Wallis test. b, Reporter RNA level (Tbp-normalized RT-qPCR) and fluorescence measured by flow in different clonal cell lines containing integrated versions of an intron from the Smc1 gene. Splice sites were either wild-type (+) or mutated (−), results were normalized to the double mutant. Comparisons made by one-way ANOVA. c, Density plot of GC-corrected steady-state RNA levels of unspliced TSS-proximal mRNA regions grouped by the presence and strength (MaxEnt score⁴¹) of a 5’SS motif (see methods). None n = 1,392, medium (MaxEnt 5–10) n = 1,792, strong (MaxEnt10+) n = 503. Both 5’SS-containing groups are significantly different from the inserts without a 5’SS motif (P < 0.0001), and from each other (P < 0.001) by Kruskal-Wallis test. d, Relative steady-state RNA levels of 10nt annotated 5’SSs with a MaxEnt score of >5, embedded into several background sequences. Only unspliced inserts (<3% spliced in nascent-RNA) were considered. Scrambled (Scr, n = 24) and antisense (AS, n = 24) versions of 5’SSs were compared to sense (S) 5’SSs (n = 50) by Kruskal-Wallis test. For all panels, **** indicates P < 0.0001, all higher P values are indicated in the plots.

Notably, 5’SS mutations had greater effects on RNA abundance than 3’SS mutations (Fig. 5a), consistent with the 5’SS having a splicing-independent role in stimulating transcription. The stronger effect of the 5’SS mutation was preserved when considering only inserts that did not contain PAS hexamers (Extended Data Fig. 5b), demonstrating that the splicing-independent 5’SS function does not rely on suppressing PAS-mediated transcript cleavage or termination.

We confirmed the effects of both 5’SS and 3’SS mutations in clonal cell lines that had a wild-type or mutant intron integrated at the reporter locus. Compared to cell lines where both SSs were mutated, restoring just the 3’SS had no effect on RNA or protein abundance (Fig. 5b). Restoring the 5’SS alone significantly increased the reporter level, while the effect of restoring both splice sites, and thereby splicing, was much larger (Fig. 5b, Extended Data Fig. 5c). This experiment supports the conclusion that the 5’SS has a splicing-independent role in stimulating transcription. Importantly, the double mutant behaves similarly to the 5’SS mutant, implying that the 3’SS only provides a positive signal in the context of a spliced intron and that the process of splicing stimulates transcription. We note that the process of splicing could strengthen the splicing-independent role of the 5’SS, or the splicing-dependent and splicing-independent stimulation could function through independent mechanisms.

Design of insert library

A library of 16,461 173nt sequences was designed, primarily consisting of sequences derived from specific regions in the mouse genome (mm10). The largest class (10,774 sequences) were sequences just downstream of transcription start sites as defined by Start-seq data (GSE43390)⁴⁷, using GENCODE M20 annotations (removing confidence level 4/5 transcripts). TSS locations were refined from annotation using TSScall with a Start-seq 5’ read count threshold of 16, search window of 250bp and joining distance of 500bp (https://github.com/lavenderca/TSScall). Using these TSSs, TSS-proximal (+6-179) and TSS-distal regions (+160-333) from different transcript classes were selected as described in Supplementary Table 1. In addition, 503 regions of annotated protein-coding transcript ends were included as described in Supplementary Table 1.

High-confidence introns were defined using GENCODE M20, selected to be between 50 and 158nt in size, and filtered based on evidence of splicing in mESCs (more details in Supplementary Table 1). The 173nt sequences included in the library start 12nt upstream of the 5’SS. For a subset of introns, multiple barcoded wild-type and mutant versions were added to the library. 1nt and 3nt mutations were made for both the 5’SS and 3’SS. A python script (https://github.com/audy/barcode-generator) was used to generate a set of 7nt barcodes with a minimum editing distance of 5 and a maximum stretch of 3. These barcodes replaced the first 7nt of each sequence.

Shorter sequences were placed in the context of the same 5 background sequences, see Supplementary Table 1. In these backgrounds, nucleotides 10 to 19 were replaced by 10nt 5’SSs (−3 to +7). A subset of 5’SSs was also inserted in the antisense direction, or was scrambled (no GGT remained after scrambling).

As controls, a set of 1100 completely randomly generated sequences were included, which cover a range of GC contents that was comparable to the regions taken from the genome. Sequences containing AWTAAA, AGGTR or an out-of-frame ATG were excluded.

Library cloning

The 173nt sequences described above were flanked by common regions to allow for amplification of the whole library and sequencing using common Illumina sequencing primers: TTTCCCTACACGACGCTCTTCCGATCTA[N₁₇₃]TAGATCGGAAGAGCACACGTCTGAACTCC. The reverse complement of these sequences was ordered as an oligo pool from Agilent Technologies, and was amplified for 10 cycles using NEBNext Ultra II Q5 (New England Biolabs). pHV141 (Oct4 uaRNA) and pHV152 (lincRNA) were linearized using AflII and the amplified library was inserted using Gibson assembly. The resulting DNA was concentrated and cleaned up using AMPure XP beads and electroporated into E. cloni 10G ELITE electrocompetent cells (Lucigen). A small fraction of electroporated cells was plated to estimate the number of transformants (~1.5M or more), the remainder was grown up overnight in liquid LB+Amp cultures and plasmids were isolated.

Library integration

For integration at the Oct4 uaRNA allele, first a gRNA was designed to target a site +157 nt downstream of the uaTSS, which overlays a SNP between 129 and CAST. Thus, this guide should only target the 129 allele in the hybrid mESCs. Using TIDE⁴⁸, the efficiency of this gRNA (pHV137) was measured to be ~75% on the 129 allele and negligible on the CAST allele. The gRNA plasmid (pHV137) and library-containing plasmid pool (pHV142) were co-transfected into the clonal cell line with the integrated TagBFP2 and mKate2 reporters at a large scale. 48 hours post-transfection, 6M GFP+ cells were sorted and propagated further. Genotyping suggested that integration had occurred in ~25% of this pool of cells. During propagation, care was taken to maintain library representation by always passaging >12M cells.

For integration at the 4930461G14Rik lincRNA allele, the same approach was used, but this locus did not include a fluorescent reporter, and the gRNA in pHV149 targeted TSS+102 specifically at the CAST allele. pHV149 and the library-containing plasmid pool (pHV154) were co-transfected into wild-type F121-9 cells. 48 hours post-transfection, 5.5M GFP+ cells were sorted and propagated further. Genotyping suggested that integration had occurred in ~33% of this pool of cells. During propagation, care was taken to maintain library representation by always passaging >12M cells.

RNA-based INSERT-seq screens

For steady-state RNA read-out, >5M library-containing cells were pelleted and resuspended in Trizol (Ambion). RNA was cleaned up using Direct-zol RNA miniprep columns (Zymo Research), treated with RQ1 DNase (Promega), then cleaned up again using a total RNA purification kit (Norgen Biotek Corp).

For nascent RNA, library-containing cells were permeabilized and run on largely as described⁴⁹, with the following changes: For permeabilization, buffers W and P were supplemented with 1 mM EGTA, and buffer F was supplemented with 1 mM EDTA. Cells were resuspended in 500 μL Buffer W, 10 mL Buffer P was added, and this was incubated on a Nutator for 1 minute, then put on ice for 6 minutes. Cells were washed once with 10mL buffer W, before being resuspended in Buffer F. Run-on was performed on 30M permeabilized cells, and in the 2X run-on reaction mix MgCl₂ was present at 10 mM, SUPERase-In at 0.4 U/μL and (biotin-)rNTPs were at 50 μM each. Run-on was performed at 37°C for 15 minutes, and Biotin-11-C/UTP (Perkin-Elmer) mixed with unlabeled rGTP and rATP (New England BioLabs). Immediately after RNA isolation using the Total RNA Purification Kit (Norgen Biotek Corp), biotinylated RNAs were bound to Streptavidin M-280 magnetic beads (ThermoFisher) and binding and washing was done as described.

Chromatin-associated RNA (Chr-RNA) was isolated as published previously⁵⁰, but 0.1% Triton-X-100 was added after pellets were resuspended in buffer B. Incubations were all performed on ice, 8 minutes in buffer A and two incubations in buffer B for 15 minutes. The washed chromatin pellets were resuspended in Trizol, RNA was extracted using Direct-zol RNA miniprep columns (Zymo Research), treated with RQ1 DNase (Promega), then cleaned up again using a total RNA purification kit (Norgen Biotek Corp).

Purified steady-state and chromatin-associated RNA and nascent RNA bound to beads was reverse transcribed with Superscript IV (Thermo Fisher) using a gene-specific primer. For experiments with the library integrated at the Oct4 uaRNA reporter locus, the RT primer annealed in mKate2, 221nt downstream of the variable region. For experiments with the library integrated at the 4930461G14Rik lincRNA locus, the RT primer annealed 287nt downstream of the variable region. cDNA from steady-state RNA was RNase-treated, cDNA from nascent RNA was removed from the beads with heat. cDNA was concentrated using 2x RNA XP beads (Beckman Coulter) to limit the number of parallel PCR reactions that had to be set up.

Sequencing libraries were generated through 17-20 PCR cycles using NEBNext Ultra II Q5 (New England Biolabs) and NEBNext Multiplex oligos, uniquely indexing each sample. Amplicons were cleaned up using 1.5x AMPure XP beads (Beckman Coulter) and concentrations were quantified using the NEBNext Library Quant Kit for Illumina, before mixing and checking the size distribution by TapeStation (Agilent). Paired-end sequencing was performed on MiSeq, NextSeq500, HiSeq4000 or Novaseq6000 (Illumina, 110+50 cycles or 150+150 cycles).

For normalization purposes, genomic DNA was extracted using QuickExtract (Lucigen) in parallel with RNA extraction, and sequencing libraries were generated from genomic DNA as described below for the sort-seq protocol. In these experiments, either the Oct4 uaRNA or the 4930461G14Rik lincRNA locus was amplified in the first round of PCR, depending on where the library was integrated.

Steady-state RNA screens after EXOSC3 depletion or non-targeting control siRNA treatment were performed similarly to other steady-state RNA screens, except that 48 hours prior to harvesting, cells were transfected with 40 nM siRNAs, using the RNAiMAX reagent (Invitrogen) in a suspension transfection. For EXOSC3 depletion, a mix of four siGENOME siRNAs was used (Catalog ID MQ-064537-01-0002, Horizon Discovery). As a control, the siGENOME Non-Targeting Control siRNA #2 was used.

Using these methods, four replicates of steady-state RNA and two replicates of nascent RNA and chromatin-associated RNA at the Oct4 uaRNA reporter locus were performed, as well as three replicates each of steady-state RNA after siRNA treatment (siEXOSC3 and control) and with the library integrated at the 4930461G14Rik lincRNA locus. We note that the steady-state RNA results inherently have a larger dynamic range than the results with nascent and chromatin-associated RNA: using RT-qPCR and comparing to a standard of in vitro transcribed RNA, we have estimated that there are ~15 copies of the reporter transcript in steady-state RNA per cell. In comparison, on the ~650bp between the annealing site of the RT primer and the PAS downstream of the reporter, we would not expect more than 1 RNAPII.

Sort-seq screen

Single-cell suspensions of library-containing cells were made in PBS with 1% FBS and 1mM EDTA. PI (0.μ/mL final concentration) was added and PI-negative cells were sorted into six bins based on the ratio between red (594 nm, 610/20) over blue (405 nm, 450/40) fluorescence using a FACSAria cell sorter (BD Biosciences). Sorting on the ratio between integrated fluorescent proteins reduces noise introduced by cell-to-cell variation in cell size and cell cycle stage. Gates were set so that the middle four bins were of equal width in the blue x red scatter plot, while the outer two bins were wider to include even cells with larger changes in protein expression. Compared to control cells in which the pHV137 Cas9+gRNA plasmid had been co-transfected with a repair template that only made a point mutation to prevent re-cleavage, the population of library-containing cells had the same median red/blue ratio, but a wider distribution.

Cells collected in each fluorescence bin, as well as unsorted cells, were pelleted and genomic DNA extracted using QuickExtract (Lucigen). Genomic DNA was also extracted from a pellet of unsorted cells. In a first round of PCR, the Oct4 uaRNA locus was amplified from these samples with NEBNext Ultra II Q5, using primers that would not amplify remaining plasmid or non-specific integrations at other loci. Amplicons were cleaned up using 0.9X AMPure XP beads and a small fraction used for a genotyping PCR. This confirmed an enrichment of insert-containing alleles in the outer bins, whereas alleles that did not integrate an insert were mostly found in the middle two bins.

The amplicons from the first round of PCR were used as template in a PCR with the NEBNext multiplex oligos to generate the sequencing library. This PCR and all subsequent steps were done as described for the RNA screen, with two changes: samples were amplified for ~22 cycles total between the first and second PCR step, and after pooling the sample was gel-extracted from a 1% agarose gel to size select. Two replicate sort-seq experiments were performed and had strong agreement (Spearman’s rho of 0.91).

Counting reads per insert

To map inserts amplified from genomic DNA, first Cutadapt v1.14⁵¹ was used to trim off adapters and low quality ends, then bowtie2 (version 2.3.4.3)⁵² was used to map against an index made of the library of 16,462 insert sequences. For bowtie2 mapping, the following parameters were set: --no-discordant --no-mixed --rdg 20,5 --rfg 20,5 --score-min L,0,−0.11. This basically disallowed gaps and allowed up to 3 mismatches on high-quality bases, with the editing distance between any two sequences in the library being >5.

For the cDNA-derived samples, splice junctions were called with STAR (version 2.7.0)⁵³, using as input a concatenated fastq file with the trimmed reads of replicates 2–4 (which had much higher coverage than replicate 1) of steady-state RNA-derived libraries from the untreated Oct4 uaRNA reporter locus screen. The following flags were used to run STAR: --outSJfilterOverhangMin 10 4 4 4 --peOverlapNbasesMin 10 --peOverlapMMp 0.08 --alignEndsType EndToEnd --outFilterMismatchNmax 5 --scoreDelOpen −10 --scorelnsOpen −10. The detected splice junctions were filtered to include only junctions on the + strand with >10 read counts and >25% of the maximum counts for other junctions detected in the same insert. A new reference genome was generated that contained both unspliced and spliced versions of the full library of sequences. Bowtie2 was run with the same settings as above to map the cDNA libraries to a bowtie2 index of this new reference.

Samtools v1.9⁵⁴ was used to output the number of reads mapped to each sequence in the ‘reference genome’.

Analysis of INSERT-seq RNA data

Data of steady-state RNA, nascent RNA and corresponding genomic DNA (gDNA) was first normalized by setting the total of mapped reads to 1. Abundances of all unspliced and spliced versions of an insert were added up for a total abundance per insert per sample. A cutoff was applied to filter out inserts with very low abundance in gDNA, for which no reliable RNA/gDNA ratio could be calculated. RNA data was first normalized to abundance in the gDNA data, and then normalized to controls, so that the median RNA/gDNA ratio of ~1100 synthetic (randomly generated) control sequences was 1. Per insert, a splicing efficiency was calculated by calculating the spliced/total ratio, where ‘spliced’ is the summed abundance of all spliced versions (0 if no spliced versions were detected). After averaging replicates and before plotting data on a log2-scale, inserts with an RNA/gDNA ratio of 0 were set to the minimum non-zero ratio for that dataset.

Sort-seq analysis

A cutoff was applied to filter out very lowly abundant inserts for which no reliable distribution could be calculated. To calculate a distribution over the bins for each insert, it has to be considered that the percentage of insert-containing cells can vary between bins and that sequencing yields no data on this percentage, necessitating a different method to calculate the scaling factor for each bin. After setting the total of mapped reads to 1 in each sequencing sample, a linear model was fit to predict the composition of the unsorted sample using the composition of each of the six bins. For the first replicate, two unsorted samples were sequenced, the replicate starting from more input material was weighed 75%, the lower input replicate was weighed 25%. The coefficients of this model (which fit the data well, with a Pearson correlation of >0.9 between predicted and observed) were used to scale bins and calculate the percentage of each insert found in each of the bins. From this the sort-seq score was calculated by assigning each bin a score of 1 through 6 and calculating a weighted average per insert.

Analysis of sequences embedded in constant background regions

Short sequences were all evaluated in the context of the same 5 background sequences. For each background, the median signal of 30 different scrambled 10-nt 5’SSs was used to normalize all sequences in that background. For RNA assays, normalization involved dividing by the median of scrambled 10-nt 5’SSs, for sort-seq this median was subtracted. The mean effect of an embedded sequence across all backgrounds was then calculated.

Analysis of barcoded wild-type and mutated introns

Original introns were included in the library with three barcodes, and 1bp and 3bp mutations of the 5’SS and 3’SS with two barcodes each. Barcodes were randomly matched up with mutant versions, allowing us to determine that none of the barcodes affected the results. Since the 1bp and 3bp mutations per splice site behaved very similarly, they were considered together as simply 5’SS mutant or 3’SS mutant. The median of all barcoded replicates across all replicate experiments was calculated as the average value for each intron version (original, 5’SS mutant, 3’SS mutant). To assess the effects of splice site mutations, introns were only considered if mutations reduced splicing efficiency to <3% in nascent RNA from >30% in original introns.

Motif searching algorithms

To find sequence motifs that may convey positive signals, the MEME³⁸ and Homer³⁹ algorithms were used to identify motifs enriched in the top-10% versus bottom-50% of TSS-proximal mRNA regions. Searches were limited to the given strand, corresponding to the sequence of the RNA. For MEME the ZOOPS Motif Site Distribution was expected, and motifs were 6-12 nt long. The Homer search was limited to motifs of 6, 8 or 10 nt.

Finding 5’SS motifs using MOODS and MaxEnt

MOODS v1.9.3⁵⁵ was used to find all 5’SS motif matches (motif SD0001.1 from the JASPAR database) in the insert sequences using a lenient cutoff (P < 0.05). Only motif matches on the given strand (RNA sequence) were considered, and the MaxEnt score⁴¹ for the sequence of each motif match was obtained. Each insert was assigned the maximum MaxEnt score of any of its 5’SS motif matches. Sequences with a MaxEnt score of >5 were considered a medium-strength 5’SS, sequences with a MaxEnt score of >10 were considered strong 5’SSs.

PRO-seq

Cell permeabilization and PRO-seq library construction was performed as described before⁴⁹, with the following modifications: 1 million permeabilized mESCs were spiked with 5% permeabilized Drosophila S2 cells. PRO-seq libraries were generated from 4 biological replicates by the Nascent Transcriptomics Core at Harvard Medical School, Boston, MA. 2x nuclear run-on buffer comprised of 10 mM Tris (pH 8), 10 mM MgCl2, 1 mM DTT, 300mM KCl, 20μM/ea biotin-11-NTPs (Perkin Elmer), 0.8ÛL SuperaseIN (Thermo), and 1% sarkosyl. The run-on reaction was performed at 37°C. The NEB 5’DNA adenylation kit was used to adenylate the 3’adapter. Adenylated 3’adapter was ligated overnight at 16°C by T4 RNA ligase 2, truncated KQ (NEB), per manufacturer’s instructions with 15% PEG-8000 final. Following the 3’adapter ligation and before the second bead binding, samples were incubated in 180μL of betaine blocking buffer (1.25M Betaine in binding buffer with 0.6μM blocking oligo (TCCGACGATCCCACGTTCCCGTGG/3InvdT/)) for 5 min at 65°C and then 2 min on ice. Following treatment with T4 polynucleotide kinase (NEB), beads were washed once each in high salt, low salt, and blocking oligo wash (0.25X T4 RNA ligase buffer (NEB), 0.3μM blocking oligo) buffers. Beads were then resuspended in 5’ adapter mix (10 pmol 5’ adapter, 30 pmol blocking oligo, water). For the 5’adapter ligation, 15% PEG-8000 was used. Eluted cDNA products were amplified for 9-11 cycles of PCR with NEBNext Ultra II Q5 master mix (NEB) and Illumina TruSeq PCR primers RP-1 and RPI-X per manufacturer instructions. Final libraries were pooled and sequenced on the Illumina NovaSeq.

PRO-seq data preprocessing was performed as described elsewhere⁴⁹. PRO-seq around uaRNA and eRNA TSSs included in the INSERT-seq library were investigated by generating count matrices from the PRO-seq bedgraph files using the make_heatmap script (https://github.com/AdelmanLab/NIH_scripts/tree/main/make_heatmap). Loci with <15 promoter read counts (PRO-seq reads of which the 5’ ends map in the TSS +/− 10nt window) were filtered out, and groups were divided up based on GC content. Metagene profiles show the mean counts per 25nt bin for each group. One outlier in group 2 of the eRNAs was removed. Elongation index was calculated as the ratio of PRO-seq density in the window from +50 to +250 downstream of the TSS (representing elongating Pol II), divided by the density from the TSS to +50 (representing promoter-proximal Pol II).

RT-qPCR

Steady-state RNA was isolated using either Trizol (Ambion) and Direct-zol columns (Zymo Research) or RNeasy columns (Qiagen). Chromatin-associated RNA was isolated as described above for INSERT-seq screens. All RNA was DNase treated and then purified using a total RNA purification kit (Norgen Biotek Corp). cDNA was generated using Superscript IV (Thermo Fisher) with random hexamer oligos. qPCRs were performed using a home-made SYBR mastermix on a Bio-Rad CFX384 qPCR instrument.

Immunoblotting

Single-cell suspensions were pelleted, and pellets were resuspended in 1x Laemmli sample buffer (BioRad) with 1:40 beta-mercaptoethanol. Sampled were boiled for 10 minutes and spun down at 13,000 g for 5 minutes. The lysates of 100k cells were run on 4–20% Mini-Protean TGX Precast Protein Gels (BioRad) according to manufacturer’s instructions. Proteins were transferred to nitrocellulose membranes for 70 minutes at 300 mA. Membranes were blocked using 5% dried milk in PBS, and incubated with anti-EXOSC3 (1:2000, A303-909A, Bethyl Labs) and anti-actin (1:2000, sc-1616, SantaCruz) in 5% dried milk in TBS-T, followed by HRP-conjugated secondary antibodies (1:10,000, 111-035-144 and 705-035-147, Jackson ImmunoResearch). After incubation with SuperSignal West Pico PLUS Chemiluminescent Substrate (BioRad) according to manufacturer’s instructions, blots were imaged using the BioRad ChemiDoc, using the ImageLab software.

Flow analysis on clonal lines

Single-cell suspensions of library-containing cells were made in PBS with 1% FBS and 1mM EDTA. Red fluorescence of Atp5a1 intron-containing lines was measured on a FACSAria cell sorter (BD Biosciences), with an excitation laser of 594 nm and 610/20 nm detector. Red fluorescence of Smc1 intron-containing lines was measured on a HS800S cell sorter (Sony), with an excitation laser of 561 nm and 600/60 nm detector.

ChIP-seq genome profile snapshots

Genomic profile snapshots were generated using the the fluff software package⁵⁶.

Published H3K4me1 data⁵⁷ was downloaded and converted to FASTQ from the SRA using the SRA toolkit (accessions SRR1202460 and SRR1202461). H3K4me3 and H3K27ac ChIP-seq data were generated as follows. mESCs were cross-linked for 5 (H3K4me3) or 10 minutes (H3K27ac) with 1% formaldehyde and quenched with 0.25μM glycine. Chromatin fragmentation was performed in sonication buffer (20mM Tris pH 8.0, 2mM EDTA, 0.5mM EGTA, 0.5% SDS, 0.5mM PMSF) with a Qsonica Q800R bath sonicator at 70% amplitude pulsing with 15 second “ON” and 45 seconds “OFF” cycles for 30 (H3K4me3) or 20 (H3K27ac) minutes total sonication time. Fragment sizes were checked by agarose gel. Fragmented chromatin from 7.5M (H3K4me3) or 2.5M (H3K27ac) mESCs was used as input for immunoprecipitation. For the H3K27ac ChIP-seq, mESCs were spiked with sonicated chromatin from 0.5M S2 cells. ChIP material was diluted with 1mL of IP buffer (20mM Tris pH 8.0, 2mM EDTA, 0.5% Triton X-100, 150mM NaCl, 10% Glycerol, 5% BSA) and pre-cleared with 10μL DynaBeads Protein A (ThermoFisher) for 1 hr at 4°C. The supernatant was then further diluted with 250μL IP buffer and incubated with 16 μL H3K4me3 antibody (EpiCypher Cat#13-0028, Lot #18303001) or 10μL H3K27ac antibody (Active Motif cat#39133, lot#22618011) overnight at 4°C. Bound chromatin was isolated with 150 μL (H3K4me3) or 50μL (H3K27ac) of pre-blocked DynaBeads Protein A for 2 hrs at 4°C. Beads were washed 1x in low salt (20mM Tris pH 8.0, 2mM EDTA, 1% Triton X-100, 150mM NaCl, 0.1% SDS), 3x in high salt (20mM Tris pH 8.0, 2mM EDTA, 1% Triton X-100, 500mM NaCl, 0.1% SDS), 1x in LiCl (20mM Tris pH 8.0, 2mM EDTA, 250mM LiCl, 1% IGEPAL, 1%[w/v] sodium deoxycholate), and 2x in TE (10mM Tris pH 8.0, 0.1mM EDTA) prior to elution with 2x 250uL of elution buffer (10% SDS, 100mM sodium bicarbonate). 0.2M NaCl was added and crosslinks were reversed by incubation at 65°C overnight. The IP mixture was treated with Proteinase K (Invitrogen) for 1 hr at 50°C. DNA was extracted with phenol:chloroform:isoamyl alcohol and precipitated with ethanol. Sequencing libraries were prepared using the NEBNext Ultra II DNA Kit (New England Biolabs) as per the manufacturer’s protocol. H3K4me3 and H3K27ac ChIP-seq libraries were sequenced on an Illumina NextSeq using a paired end 160 (read1:80, index:6, read2:80) and 80 cycle (read1:43, index:6, read2:42) kits respectively.

Data was mapped to mm10 (GRCm38) using Bowtie v1.2.2⁵⁸, using the parameters -v2 -X1000 –best and -k1 (for H3K4me1/me3) or -m1 (for H3K27ac). Using the extract_fragments script (https://github.com/AdelmanLab/NIH_scripts/tree/main/extract_fragments), reads were filtered for insert fragment sizes of 75-500bp, duplicate reads were removed using Samtools v1.3 and bedgraphs mapping the midpoints of deduplicated fragments were generated. For H3K4me1 data, which was sequenced as single-read, fragment lengths of 150bp were assumed. Outputted bedgraphs were at 1-bp resolution for H3K4me1/me3 data and in 25-bp bins for H3K27ac data. Bedgraphs from replicate experiments were merged using the bedgraph2stdbedgraph script (https://github.com/AdelmanLab/NIH_scripts/tree/main/bedgraphs2stdBedGraph), and 1-bp merged bedgraphs were rebinned to 25-bp bins using binBedGraph (https://github.com/benjaminmartin02/binBedGraph).

TT-seq

TT-seq was performed essentially as described⁵⁹, with some minor modifications. Briefly, mESCs, treated with 0.03 % DMSO for 4 hours, were labelled with 500 μM 4-thiouridine (Sigma) for 20 minutes. Cells were counted and collected in Trizol, to which was added 5% S2 cells, labelled with 4sU for 2 hours. RNA was chloroform extracted, DNase-treated, and chloroform extracted. To 60 μg total RNA was added 2X fragmentation buffer (final concentration: 75 mM Tris Cl, pH 8.3 112.5 mM KCl, and 4.5 mM MgCl₂) and heated to 95 °C for 5 minutes. Fragmentation was stopped by adding EDTA to 50 mM and placing samples on ice. RNA was ethanol precipitated and resuspended in H₂O. Fragment sizes were checked on the TapeStation (peak size of ~800 nt). Biotinylation reaction was performed with 0.025 mg/mL MTSEA-Biotin XX (Biotum) in reaction buffer (20% N,N-Dimethylformamide, 1mM EDTA, 10 mM HEPES pH 8) for 45 minutes in the dark. Labelled RNA was chloroform extracted and ethanol precipitated, and resuspended in 90 μL H₂O. 75 μL DynabeadsTM M-280 Streptavidin (ThermoFisher) were pre-washed with decon solution (0.1M NaOH + 50mM NaCl), 2X 100mM NaCl, 2X High Salt buffer (100mM Tris Cl pH 7.4 10mM EDTA pH 8, 1M NaCl, 0.05% (v/v) Tween 20), and resuspended in High salt buffer. Labelled RNA was denatured by heating to 65 °C for 5 minutes and placing on ice for 2 minutes. 10 μL High salt buffer was added to labelled RNA. Prewashed beads were placed on magnet, supernatant discarded, and then beads were resuspended in the labelled RNA. Beads and RNA were rotated for 30 minutes in the dark. Beads were then washed 4X for 1 minute with high salt buffer. Labelled RNA was eluted 2X with 100 mM DTT (1,4-Dithiothreitol). RNA was cleaned up with microelute columns (Norgen). Sequencing libraries were prepared from 250 ng of labelled RNA using the TruSeq Stranded Total RNA sequencing kit with RiboZero rRNA depletion, shortening the fragmentation step to 3 minutes and using 8 cycles of PCR amplification.

Reads were trimmed for a minimum quality score of 20 using a custom script (https://github.com/AdelmanLab/NIH_scripts/tree/main/trim_and_filter_PE), adapter sequences were removed using cutadapt version 1.14⁵¹. Reads were first mapped to dm6 using bowtie⁵⁸, before aligning to mm10 using STAR (version 2.7.3a)⁵³. Strand-specific bedgraphs were generated from deduplicated BAM files using STAR.

We thank Kanae Sasaki for her help in optimizing the run-on protocol for screening purposes and Emily Kaye for discussions on library design. We thank Stephen Buratowski for useful discussions on the project, and Daria Shlyueva and Torben Heick Jensen for feedback on the manuscript. We are also grateful to the Flow Cytometry Facility at HMS’ Department of Immunology for cell sorting help and advice, the HMS Nascent Transcriptomics Core for PRO-seq library construction and the HMS Biopolymers Facility and The Bauer Core Facility at Harvard University for next-generation sequencing. This research was supported by the European Molecular Biology Organization (ALTF 531-2017 to H.V.), Human Frontier Science Program (LT000651/2018-L to H.V.), the National Institutes of Health (NIH R01 GM139960 to K.A.), startup funding from Harvard Medical School to K.A, the National Science Foundation Graduate Research Fellowship (DGE1745303 to C.A.M) and the Canadian Institutes of Health Research (Banting fellowship to B.J.E.M.).