TT-seq maps the human transient transcriptome

See allHide authors and affiliations

Science  03 Jun 2016:
Vol. 352, Issue 6290, pp. 1225-1228
DOI: 10.1126/science.aad9841

TT-Seq maps a transient transcriptome

RNA expression is related to protein abundance and cellular function. However, the amounts of RNA generated at any one time-point have been difficult to determine. Schwalb et al. developed a method, transient transcriptome sequencing (TT-Seq), to collect and sequence all RNA segments synthesized over 5 minutes. Because 5 minutes is not long enough to fully degrade even the most transient RNA, this method can detect the synthesis of most RNA without bias. Applying this method to human K562 cells, TT-Seq detected thousands of noncoding transcripts, providing a snapshot of RNA synthesis rates and RNA half-lives, and full-length maps of short-lived RNAs such as enhancers and short intergenic noncoding RNAs.

Science, this issue p. 1225


Pervasive transcription of the genome produces both stable and transient RNAs. We developed transient transcriptome sequencing (TT-seq), a protocol that uniformly maps the entire range of RNA-producing units and estimates rates of RNA synthesis and degradation. Application of TT-seq to human K562 cells recovers stable messenger RNAs and long intergenic noncoding RNAs and additionally maps transient enhancer, antisense, and promoter-associated RNAs. TT-seq analysis shows that enhancer RNAs are short-lived and lack U1 motifs and secondary structure. TT-seq also maps transient RNA downstream of polyadenylation sites and uncovers sites of transcription termination; we found, on average, four transcription termination sites, distributed in a window with a median width of ~3300 base pairs. Termination sites coincide with a DNA motif associated with pausing of RNA polymerase before its release from the genome.

Transcription of eukaryotic genomes produces protein-coding mRNAs and diverse noncoding RNAs (ncRNAs), including enhancer RNAs (eRNAs) (1, 2). Most ncRNAs are rapidly degraded, difficult to detect, and thus far have not been mappable in their full range. Mapping of transient RNAs is required, however, for analysis of RNA sequence, function, and fate.

We developed transient transcriptome sequencing (TT-seq), a protocol that maps transcriptionally active regions and enables estimation of RNA synthesis and degradation rates. TT-seq is based on 4sU-seq, which involves a brief exposure of cells to the nucleoside analog 4-thiouridine (4sU) (Fig. 1A) (3). 4sU is incorporated into RNA during transcription, and the resulting 4sU-labeled RNAs are isolated and sequenced. 4sU-seq is more sensitive than RNA-seq in detecting transient RNAs. However, 4sU-seq fails to map human transcripts uniformly, because only a short 3′ region of nascent transcripts is labeled during a 5-min exposure to 4sU, and the long preexisting 5′ regions dominate the sequencing data. To remove this 5′ bias, TT-seq uses RNA fragmentation before isolation of labeled RNA fragments (Fig. 1A). Thus, TT-seq measures only newly transcribed RNA fragments and provides the number of polymerases transcribing a genomic position within 5 min.

Fig. 1 TT-seq enables uniform mapping of the human transient transcriptome.

(A) Workflow of 4sU-seq and TT-seq protocols. (B) Metagene coverage, comparing TT-seq with other transcriptomic methods. Shown is the average coverage for 2323 TUs lacking paused and active genes (10) around the first TSS (left), the 5′ splice site of intronic sequences > 10 kbp (first intron excluded), and the last pA site. Signals are relative to the maximum signal in the first kilobase pair from the first TSS.

When applied to human K562 cells, TT-seq samples newly transcribed regions uniformly, whereas 4sU-seq produces a 5′ bias (fig. S1A). The coverage of short-lived introns with respect to exons is estimated (4) to be 60% for TT-seq, whereas it is 23 and 8% for 4sU-seq and RNA-seq, respectively (figs. S1A and S2). TT-seq is highly reproducible (fig. S3) and enables complete mapping of transcribed regions, complementing the GRO-cap (5) and CAGE (6) protocols, which detect RNA 5′ ends (Fig. 1B). TT-seq monitors RNA synthesis, whereas protocols such as PRO-seq (7), NET-seq (8), and mNET-seq (9) detect RNAs attached to polymerase. Therefore, the latter protocols yield peak signals near the promoter where polymerase pauses (Fig. 1B), whereas TT-seq does not. For paused and active genes (10), TT-seq reveals higher rates of RNA synthesis near the promoter relative to other regions (fig. S1B).

Using TT-seq data and the segmentation algorithm GenoSTAN (4, 11), we identified 21,874 genomic intervals of apparently uninterrupted transcription (transcriptional units, TUs) (Fig. 2 and fig. S4A). TT-seq is highly sensitive, recovering 65% of transcription start sites (TSSs) obtained by GRO-cap (overlapping annotations within ± 400 bp) (5). A total of 8543 TUs overlapped GENCODE annotations (12) in the sense direction of transcription (50% reciprocal overlap of annotated regions; fig. S4B). This analysis detected 7810 mRNAs, 302 long intergenic noncodinzg RNAs (lincRNAs), and 431 antisense RNAs (asRNA). The 2916 TUs that shared less than 50% of their length with GENCODE annotations were not classified. The remaining 10,415 TUs (48%) represented newly detected ncRNAs that we characterized further.

Fig. 2 Annotation of RNAs mapped by TT-seq.

Example genome browser view showing RNAs from 5 of 7 transcript classes and 3 of 18 chromatin states [chromosome (chr.) 18, 9.00 to 9.54 million base pairs (Mbp)] (4). Arrows indicate the direction of transcription. Units on the y axes are read counts per 200-bp bin.

Transcripts arise from promoters but also from enhancers, which are regulatory elements with characteristic chromatin modifications (13, 14). To detect chromatin regions comprising putative enhancers and promoters (chromatin states), we applied GenoSTAN (11) to ENCODE ChIP-seq (chromatin immunoprecipitation–sequencing) data (15) for the coactivator p300 and a series of histone modifications (H3K27me3, H3K36me3, H4K20me1, H3K4me1, H3K4me3, H3K9ac, and H3K27ac) and to deoxyribonuclease I hypersensitivity data (fig. S5A). Of the resulting strong enhancer state regions, 81% overlapped at least one TSS from GRO-cap (5) and 68% overlapped a polymerase II (Pol II) peak (15), compared with 52 and 37%, respectively, for ENCODE enhancer states (fig. S5, B to D).

The 10,415 nonannotated TUs were classified based on GenoSTAN-derived chromatin states and their positions relative to known GENCODE annotations (fig. S6A). TUs within 1 kilo–base pair (kbp) of a GENCODE mRNA TSS included 685 upstream antisense RNAs (uaRNAs) (16) and 778 convergent RNAs (conRNAs) (8). The 3115 TUs on the strand opposite an mRNA were classified as asRNAs when they were more than 1 kbp away from the GENCODE TSS. Remaining TUs were grouped according to their GenoSTAN chromatin state at their TSS. The 2580 TUs that originated from promoter state regions (4) were classified as short intergenic ncRNAs (sincRNAs) (fig. S6B). Most sincRNAs (67%) were located within 10 kbp of a GENCODE mRNA TSS. The remaining 3257 TUs originated from enhancer state regions (4) and were classified as eRNAs (13, 14). The newly mapped ncRNAs are short (fig. S6C). On average, lincRNAs are five times as long as sincRNAs, and eRNAs have a median length of ~1000 nucleotides.

Kinetic modeling of TT-seq and RNA-seq data enabled us to estimate rates of RNA synthesis and degradation (Fig. 3 and fig. S7A) (4). We estimated rates of phosphodiester bond formation or breakage at each transcribed position and averaged these within TUs, thus obtaining estimates of relative transcription rates and RNA stabilities (4). We found that mRNAs and lincRNAs had the highest synthesis rates and longest half-lives. We determined a median mRNA half-life of ~50 min, compared with a previous estimate of ~139 min (17). Other transcript classes had low synthesis rates and short half-lives, explaining why short ncRNAs are difficult to detect. eRNAs had half-lives of a few minutes, consistent with prior data (17). Short RNA half-lives correlated with a lack of secondary structure (fig. S7B). The folding energy of eRNAs was comparable to the genomic background level (fig. S7C), and only 10% of their sequence was predicted to be structured, compared with 52% in mRNAs (fig. S7D).

Fig. 3 Estimated RNA half-lives for different transcript classes.

Black bars represent medians, boxes represent upper and lower quartiles, and whiskers represent 1.5 times the interquartile range (4).

We further found differences in transcription from promoters versus from enhancers (2). Enhancers showed lower occupancy of initiation factors TBP and TAF1 than mRNA promoters did (12- and 3.5-fold less, respectively; P < 10−16, Fisher’s exact test), whereas TFIIB and TFIIF had similar occupancies in enhancers and promoters. Occupancies were also similar for factors involved in polymerase pausing, such as NELF-E and the P-TEFb subunit cyclin T2 (fig. S8A). Synthesis of eRNAs terminated early (fig. S6C), probably because eRNAs are not enriched in U1 small nuclear ribonucleoprotein–binding sites (U1 signals; GGUAAG, GUGAGU, or GGUGAG) that can counteract early termination and lead to RNA stabilization (1820). eRNAs contained U1 signals at the genomic background level (47%), whereas mRNAs were enriched (69%; P < 10−16) (fig. S8B). In all transcript classes, longer RNAs were enriched with U1 signals in the first 1000 nucleotides (fig. S8C), suggesting that evolution of stable RNAs generally involves acquisition of U1 signals.

TT-seq also enabled us to uncover transcription termination sites (TTSs). TT-seq detected transient RNA downstream of the polyadenylation (pA) site (Figs. 1B and 4 and fig. S9). Such RNA is difficult to detect because RNA cleavage at the pA site leads to an unprotected 5′ end and RNA degradation by the XRN2 exonuclease (21, 22). For a total of 6977 mRNA genes, we derived, on average, four TTSs (figs. S10 and S11) (4). TTSs were located within a termination window that extended from the last pA site to an “ultimate TTS,” where RNA coverage dropped to background levels (Fig. 4, A to C). The termination window had a median width of ~3300 bp and could be up to 10 kbp wide (fig. S10D), consistent with Pol II occupancy data (Fig. 4D) (23). For the 5113 TTSs with the strongest drop in TT-seq signal (4), Pol II peaks were obtained by PRO-seq (7), NET-seq (8), and mNET-seq (9) (fig. S12), indicating that Pol II pauses at the TTS.

Fig. 4 Transcription termination sites.

(A) TT-seq coverage for two replicates (red and blue) downstream of the pA site in the ALDH1B1 gene locus. Arrows indicate TTSs obtained from segmentation (solid arrow, ultimate TTS). The last annotated pA site (per GENCODE) was aligned at zero. (B) Generic gene architecture. The first TSS was aligned at zero, and the last pA site was set at a rescaled distance of 5000 bp from the TSS (the real median distance is 24,079 bp for 6977 investigated genes) (4). The ultimate TTS is depicted at a median distance of 3359 bp from the last pA site (rescaled). (C) TT-seq coverage over quantiles 0.05 to 0.95 (pink area; black line, median), rescaled and aligned as in (B). (D) Pol II occupancy determined by the ChIP-exonuclease method (23), visualized as in (C) (black line, median; white line, upper quartile). (E) (C/G)(2–6)A and (T/A)(3–6) sequence count within ±100 bp of ultimate TTSs. (F) PWM (position weight matrix) logo representation of nucleotides at positions –9 to +2 around the ultimate TTS (position 0). (G) Predicted polymerase states at the T-rich stretch downstream of the TTS (top) and after backtracking to the TTS (bottom) (gray, DNA; red, RNA).

The derived TTSs are strongly enriched for the sequence (C/G)(2–6)A (window ± 5 bp; P <10−16, Fisher’s exact test; odds ratios, 2.98 and 1.63 for C3A and G3A, respectively; subscript numbers denote the number of nucleotide repeats) (Fig. 4E and figs. S13 and S14) (4). This sequence can contain up to six cytosines or guanines (Fig. 4F and fig. S10E). A G3A element has also been found in a known termination signal (24). The C(2–6)A and G(2–6)A sequences are generally followed by a T-rich [T(3–6)] or an A-rich [A(3–6)] stretch, respectively (P < 10−16; odds ratios, 2.39 and 1.24 for T4 and A4, respectively), that is located, on average, 15 bp downstream of the TTS (Fig. 4E). Such sequences were found at TTSs of all TUs (fig. S15) and can even be derived from published data (fig. S12B). In summary, the detected TTSs were highly enriched with the consensus motif (C/G)(2–6)ANx(T/A)(3–6), where Nx is a short stretch of nucleotides.

To test for the in vivo functionality of the derived TTS motif, we transfected expression plasmids into K562 cells that either lacked or contained four C3AN8T4 or C7AN8T4 motifs within 600 bp downstream of the pA site (fig. S16A and tables S1 to S4) (4). When the TTS motifs were present, significantly less RNA was detected downstream of the motifs, indicating termination of a fraction of polymerases (Wilcoxon test; fig. S16B). This experiment supports the functionality of the derived TTS motif in vivo. Termination depended on an upstream pA signal (fig. S16B), consistent with an occurrence of the motif in gene bodies, where they do not lead to transcription termination.

Transcription over the (C/G)(2–6)ANx(T/A)(3–6) sequence is predicted to destabilize the polymerase complex (24, 25) because the melting temperature (4) of the DNA-RNA hybrid is low in the T/A-rich region (fig. S17). This may trigger backtracking and trap polymerase at the TTS (Fig. 4G). At the TTS, the hybrid is C/G-rich and stable, and RNA may be cleaved from its 3′ end to yield a terminal A residue. Polymerase can then be released from DNA by XRN2 (21, 22). TT-seq has afforded insights into the determinants of human genome transcription and provides a complementary tool for transcriptome analysis.

Supplementary Materials

Materials and Methods

Figs. S1 to S17

Tables S1 to S4

References (2639)

References and Notes

  1. Materials and methods are available as supplementary materials on Science Online.
  2. Acknowledgments: The sequencing data and the annotation file have been deposited in the Gene Expression Omnibus database under accession code GSE75792. We thank H. Blum, S. Krebs, and A. Graf (Laboratory for Functional Genome Analysis, Gene Center, Ludwig-Maximilians-Universität München) for help with sequencing and J. A. Feuillet, J. Soeding, A. Sawicka, and C. Bernecky for help and discussions. K.F. was supported by the Center for Innovative Medicine (CIMED) at Karolinska Institutet and by the Science for Life Laboratory (SciLifeLab) in Stockholm. C.D. was supported by a Deutsche Forschungsgemeinschaft (DFG) fellowship at the Graduate School of Quantitative Biosciences Munich. A.T. was supported by a German Federal Ministry of Education and Research (BMBF) e:Bio grant and by DFG grant SFB 680. J.G. was supported by the Bavarian Research Center for Molecular Biosystems and the Bundesministerium für Bildung und Forschung, Juniorverbund in der Systemmedizin “mitOmics” (grant FKZ 01ZX1405A). P.C. was funded by the Advanced Grant TRANSIT of the European Research Council, the DFG, the Volkswagen Foundation, CIMED, and SciLifeLab.
View Abstract

Stay Connected to Science

Navigate This Article