Tunable protein synthesis by transcript isoforms in human cells

1Eukaryotic genes can produce a staggering diversity of messenger RNA products. In humans, there is a median of five annotated transcript isoforms per gene, with more than 80 annotated isoforms in some cases (Ensembl release 75). These transcript isoforms arise from the combined action of alternative transcription initiation, splicing and cleavage, and polyadenylation, which are often cell type specific. As each of these transcript isoforms may in principle harbor a unique set of translational control elements, gene-level expression, which averages out the translational potential of each individual transcript isoform, may not be an accurate measure of protein levels. Since changes in both translation and splicing are linked to numerous human disorders, it is critical to understand the relationship between the two.

2Eukaryotic mRNAs are decorated with diverse sequence features that can control translation, which can vary between transcript isoforms. Classic examples of translational control include upstream ORFs (uORF)-mediated translational control of yeast GCN4, protein binding such as the iron regulatory protein, and the action of micro-RNAs or DEAD-box proteins such as eIF4A and Ded1. Alternative 5′ leader sequences, uORFs, and select tandem 3′ untranslated region (UTR) isoforms have been demonstrated to influence protein production. Any of these features may in principle be different between transcript isoforms, but the prevalence and dynamic range of isoform-specific translational control across the human genome is currently unknown.

3Previous work measuring genome-wide translation in human cells has focused largely on the relationship between gene-level mRNA abundance and protein levels, which is blind to the contribution of transcript isoforms. Ribosome profiling is not well-suited for measuring transcript isoform-specific translation, primarily due to the short ~30 bp length of ribosome-protected fragments. Prior attempts to characterize isoform-specific translation have measured the effects of 5′ end diversity in yeast and 3′ end diversity in mouse cells, or splicing differences between cytoplasmic and aggregate polysomal mRNAs. However, sequencing just the ends of transcripts cannot distinguish between transcript isoforms of the same gene harboring degenerate termini. In addition, aggregating polysome fractions averages lowly- and highly-ribosome-associated messages. Therefore, a different strategy is required to understand how the diversity of the human transcriptome impacts translational output.

4Here, we adapt a classic approach of polysome profiling coupled with global gene expression analysis to measure transcript-isoform specific translation using deep sequencing, which we term Transcript Isoforms in Polysomes sequencing (TrIP-seq). By using high gradient resolution and sequencing depth, this approach yields polysome profiles for over 60,000 individual transcript isoforms representing almost 14,000 protein coding genes. We observe frequent intron retention on ribosome-associated transcripts, even in high-polysome fractions, identifying a population of retained but not nuclear-detained introns. Properties of 3′ untranslated regions predominate over the 5′ leader sequence as the driving force behind differential polysome association for transcript isoforms of the same gene among the transcript features tested. We show that regulatory sequences differentially included in transcript isoforms of the same gene are modular and can trigger differences in the translation of reporters spanning two orders of magnitude. These findings provide a lens through which to ascribe functional consequences to RNA-seq-generated transcriptomes. Moreover, TrIP-seq analysis uncovers regulatory elements that can be utilized to tune translation of synthetic messages robustly in cells.

5We determined the ribosomal association of transcript isoforms by sequencing transcripts cofractionating with different numbers of ribosomes with sufficient depth to determine isoform abundances, as was performed at the gene level in yeast. We treated HEK 293T cells with cycloheximide to stall translation and fractionated the cytoplasm into ribosome-containing samples including one to eight or more ribosomes (Figures 1A and Figure 1—figure supplement 1A; see Materials and methods for details). We made RNA sequencing libraries from each fraction in biological duplicate and obtained transcript-level abundances using the Cufflinks suite (Figure 1—source data 1 and 2;). Clustering of the samples recapitulates the gradient order (Figure 1B), indicating the polysome profile was accurately fractionated. Four subgroups emerge from this clustering: the 80S (monosome), low polysomes (two-four ribosomes), high polysomes (five-eight+ ribosomes), and total cytoplasmic RNA. We tested the robustness of the clustering of samples by computing the average Jaccard distance between clusters from data subjected to three different resampling methods, which was ≥0.75, suggesting stable clusters (Figure 1—figure supplement 1G; Materials and methods). This suggests that in cells, transcript isoforms are predominantly poorly- or highly-ribosome associated, causing low polysomes to cluster away from high polysomes, which could be have numerous biological origins including highly abundant short ORFs.

6To confirm the sequencing data, we first analyzed beta-actin (ACTB), which is known to be heavily translated with many ribosomes on each message, and migrates with high polysomes in a gradient (Figure 1C;. However, ACTB only has one transcript isoform, so we examined ATF4, which has one dominant isoform exhibiting low-polysome association in TrIP-seq data and by RT-PCR (Figure 1C), consistent with it being translationally repressed in the absence of cellular stress. We then analyzed eukaryotic elongation factor 1 beta 2 (EEF1B2), which has three isoforms with the same coding sequence but alternative 5′ leaders as a representative gene that exhibits isoform-specific ribosome association (Figure 1—figure supplement 1B). The isoform-specific polysome profiles and RT-PCR conducted from polysome fractions agree qualitatively, validating the accuracy of the isoform-level quantifications (Figure 1D). Weaker amplification of the longer isoform (EEF1B2-003) in high polysomes may be due to PCR bias towards smaller amplicons, and replicate experiments show EEF1B2-003 present in high polysomes (Figure 1—figure supplement 1C). We additionally validated the TrIP-seq data using qRT-PCR from polysome fractions and find that the two measurements agree for both EEF1B2 and SRSF5 (Figure 1—figure supplement 2). Both coding mRNAs and long noncoding RNAs (lncRNAs) associate with ribosomes, so we measured their polysomal abundance. Protein-coding genes and isoforms are found across the polysome, but noncoding and lncRNA genes are predominantly in the low polysome fractions indicating that they are generally weakly ribosome associated or in other large macromolecular complexes (Figure 1E and Figure 1—figure supplement 1D).

7To determine whether TrIP-seq is a measure of translation as opposed to cryptic association with large macromolecular complexes or stalled ribosomes, we compared TrIP-seq data to ribosome profiling and proteomics datasets. First, we reanalyzed a ribosome profiling dataset from HEK 293T cells and compared the number of ribosome-protected fragments to a weighted sum of the TrIP-seq polysome reads (Materials and methods; nribo ≥2). We observe a strong correlation between ribosome profiling and TrIP-seq at the gene level (RS = 0.88; Figure 1F). However, at the isoform level, the correlation decreases (RS = 0.47), which is worse than the correlation between TrIP-seq replicates (RS = 0.90; Figure 1—figure supplement 1E), suggesting the discrepancy is not due to variability in TrIP-seq but instead is likely due to known issues quantifying transcript isoforms in ribosome profiling. We next computed the translation efficiency (TE) of ribosome profiling and TrIP-seq data by dividing by cytoplasmic RNA levels. TE measured by the two methods correlated with RS = 0.24, perhaps due to variability in culture conditions or in either technique. We compared the data of Sidrauski et al to an unpublished ribosome profiling dataset in 293T cells from the Yoon-Jae Cho lab. Ribosome-protected fragments between experiments correlate with RS = 0.78, while TE has a substantially lower correlation of RS = 0.41, suggesting biological or technical variability disproportionately affects TE. We then compared TrIP-seq data to protein abundances from HEK 293T cells, and find that the two datasets correlate (RS = 0.57, Figure 1—figure supplement 1F,, which is better than RNA-seq (RS = 0.45; not shown). The observations that most translational stalling in mammalian cells is transient and translation elongation rates are homogeneous, TrIP-seq appears to underestimate rather than overestimate protein abundance (Figure 1—figure supplement 1F), and TrIP-seq correlates well at the gene-level with ribosome profiling (Figure 1F) suggest that TrIP-seq primarily measures translating ribosome association as opposed RNAs bound to other large complexes.

8To extract global trends in isoform-specific translation, we hierarchically clustered transcript isoform polysome profiles and selected eight clusters that are representative of general trends in the data (Figure 2; Materials and methods). The depth of sequencing (Figure 1—figure supplement 3), augmented by the fractionation strategy, enables detection of 62,703 transcript isoforms in the polysome profile (Figure 2—source data 1). Isoforms in the observed clusters exhibit diverse average patterns across polysomes (Figure 2A,B), from clusters 1 and 2, which contain isoforms primarily in high polysomes, to cluster 3 with isoforms in the middle, to clusters 6 and 7 where isoforms are in low polysomes. Independent clustering of the two biological replicate TrIP-seq datasets shows that the high- and low-polysome clusters (1, 2 and 6) appear more robust when comparing between averaged and individual replicate clusterings (Figure 2—figure supplement 1A). Many clusters have similar total polysome abundance but different distributions, indicating that to obtain accurate measurements of isoform-specific translatability it is crucial to fractionate the polysome profile.

9Surprisingly, the poorly translated cluster 7, which has ~9,000 transcripts, contains a similar number of annotated retained intron or protein-coding isoforms (Figure 2C). We queried the gene ontology (GO) terms associated with each cluster and found that the poorly translated clusters 6 and 7 are enriched for translation and splicing genes (Figure 2—figure supplement 1B), implicating the alternative splicing-nonsense-mediated decay (AS-NMD) pathway. However, not all retained intron transcripts are subject to NMD, and we also observe retained intron transcripts in the best-translated clusters (1 and 2), indicating at least some of the introns may be translated. In sum, clustering of isoform abundance distributions across polysomes immediately provides insight into the diverse patterns of ribosome association by individual isoforms and transcript types in cells, and reveals intron-rich clusters of transcripts associated with polysomes that escape nuclear detention.

10We reasoned we could extract features regulating the translation of a gene by comparing transcript isoforms of the same gene that are well- or poorly-translated. We selected 24 features to explore that are involved in translational control, such as the length of the coding sequence and untranslated regions, predicted secondary structure, and microRNA binding sites. As their average polysome profiles and composition are broadly similar (Figure 2), we merged clusters 1 and 2 to generate a larger pool of high-polysome isoforms (Figure 3A), and compared these to cluster 6, representing poorly translated isoforms in low polysomes. We then compared transcript isoforms of the same gene that are in high- and low-polysome clusters, or gene-linked isoforms, to extract transcript features that influence translation. The number of gene-linked isoforms per feature per set varies between 569 and 6491. We then measured the effect size and calculated statistical significance for all features between gene-linked isoforms (Figures 3B and Figure 3—figure supplement 1A; Materials and methods;).

11We find that longer coding sequences and highly abundant versions of gene-linked isoforms are biased towards high polysomes, likely because shorter coding sequences cannot accommodate as many ribosomes (Figure 3B and Figure 3—figure supplement 1A). However, the length of the coding sequence is not the sole determinant of polysome association, since the ribosome density (measured by dividing the weighted sum of TrIP-seq polysome reads by the number of cytoplasmic reads) is also higher in gene-linked isoforms found in high polysomes (Figure 3B). We also find that transcripts from gene-linked isoforms with more exons also tend to be better translated, as was observed using reporter genes. This is not due to overall transcript length, as there is no significant difference between transcript length in gene-linked isoforms, reflecting a positive influence of splicing on translation. Highly translated isoforms contain fewer rare codons on average, but contain more stretches of rare codons, perhaps indicating a need for translational pausing (Figure 3B). The observation that exons promote translatability has been shown for select genes, and here we show that this extends across the human genome.

12One of the strongest effects seen on polysome association of gene-linked isoforms comes from the length and content of the 3′ UTR (Figure 3B and Figure 3—figure supplement 1A). Specifically, we find that gene-linked 3′ UTRs in low polysomes are considerably longer than those on high polysomes (mean length 1551 nt versus 982 nt). There are numerous regulatory elements contained within 3′ UTRs, including microRNA binding sites, AU-rich elements, and protein-binding sites. The fraction of the 3′ UTR containing AU-rich elements is increased in low polysome gene-linked isoforms, possibly due to translational repression. Conserved predicted binding sites for miRNAs are also more abundant in poorly translated gene-linked isoforms, but this could be due to a correlation between 3′ UTR length and the number of miRNA binding sites. We, therefore, filtered the miRNAs to those bound to AGO1 or AGO2 in HEK 293T cells. In this HEK 293T-expressed miRNA set, the difference for both the number of conserved miRNA binding sites and their score increases between gene-linked isoforms, suggesting that miRNA binding may be functionally relevant in this context (Figure 3B and Figure 3—figure supplement 1A). Increased AU-rich elements and miRNA binding should lead to decreased half-life for low polysome isoforms, which is borne out in half-life comparisons (Figure 3B and Figure 3—figure supplement 1A; [).

13Surprisingly, gene-linked isoforms in high polysomes have increased predicted structure both at the cap and in 75-nucleotide windows in the 5′ leader. There are at least three possibilities that could explain this result. First, RNA structures can both repress or promote translation in a context dependent manner, so it is possible that some isoforms are being driven to high polysomes by recruitment of transacting factors. Second, recent studies of global RNA secondary structure have consistently observed different structure in cells than in vitro. Lastly, HEK 293T cells may translate messages with inhibitory 5′ leaders more efficiently than other cell types. In contrast, we find no significant dependence for the length of the 5′ leader region on polysome association based on this comparison (Figures 3B and Figure 3—figure supplement 1A). However, isoforms containing 5′ leaders over 1000 nucleotides long are poorly ribosome-associated relative to shorter 5′ leaders (Figure 3—figure supplement 1B) and cluster 7, which is associated with few ribosomes, contains longer 5′ leaders (Figure 3—figure supplement 1C). It is likely that more subtle features of the 5′ leader also influence translatability, as in the examples in Figure 1D.

14Upstream open reading frames (uORFs) can positively or negatively influence translation. We, therefore, counted the number of experimentally determined uORFs in each cluster and observe surprisingly complex behavior. Cognate uORFs (those starting with ATG) promote or repress translation, while noncognate uORFs generally promote translation. Along these lines, we find that cluster 1 (high polysomes) is enriched for cognate and noncognate uORFs while cluster 6 (low polysomes) is enriched for cognate uORFs. However, cluster 2 (high polysomes) is enriched for only noncognate uORFs, indicative of the complex and idiosyncratic behavior of uORFs. Taken together, we find predominant influences for the 3′ UTR and the number of introns in determining polysome occupancy of gene-linked transcript isoforms in human cells, while a diversity of other features can influence translatability to a lesser extent.

15RNA-seq characterizes the isoform diversity of a sample, but not the functional consequences of this diversity. We reasoned that TrIP-seq data could be used to predict which isoform changes are likely to lead to translation changes in other systems. We chose a human embryonic dataset containing 124 individual cells from seven preimplantation developmental stages and reprocessed these data to directly compare to TrIP-seq data (Materials and methods;). We mapped transcripts expressed during each developmental stage onto the TrIP-seq clusters, to attempt to gain insight into global translational properties in human embryos (Figure 4A). In embryos, our analysis predicts a shift towards the extremes of translation, with an increase in the percent of transcript isoforms that are both highly (e.g. clusters 1 and 2) and lowly (e.g. cluster 6) translated isoforms in HEK 293T cells (Figure 4A). Localized translation is widespread in development, and it is possible that the observed increase in poorly translated isoforms reflects a greater need for translational control. It is also possible that translation in early embryos is differentially regulated than in HEK 293T cells, which is now testable by applying TrIP-seq to other cell types, yielding cell type-specific translational control programs.

16We then collected genes with isoforms that change between early and late embryonic stages. Clustering of abundances for 45,895 isoforms across embryogenesis yielded seven clusters with varying profiles, which we then analyzed by comparing transcripts of the same gene located in different embryonic clusters as before (Figure 4—figure supplement 1A). These isoforms were filtered by those that move between low- and high-polysome clusters in the TrIP-seq data, yielding 366 isoform pairs belonging to 270 genes that are developmentally regulated and exhibit differential translation in HEK 293T cells. For example, the CSDE1 gene expresses one isoform in oocytes, early development and in human ES cells, which is poorly translated (CSDE1-002), but a second appears following fertilization that is well-translated (CSDE1-007; Figure 4B,C). As global zygotic transcription begins in the two-cell stage, it is possible that the alternative CSDE1 isoform is supplied by the spermatid or that it may be precociously transcribed. We also present three other genes (EIF4A2, RNF170, and TBC1D15) that show stage-specific expression in the embryo data and differential polysome association in HEK 293T cells (Figure 4—figure supplement 2). Therefore, widespread changes at the transcriptome level during human embryogenesis may produce concomitant changes in protein production, which can now be predicted using TrIP-seq data.

17Engineering the translation of a transcript without altering its coding sequence is desirable when introducing exogenous mRNA to cells or patients. We hypothesized that the regulatory features we discovered by comparing highly to lowly translated isoforms (Figure 3B) should be transferrable to an arbitrary reporter. Nine different 5′ leaders and eight different 3′ UTRs derived from gene-linked isoforms (Figure 3) were appended onto Renilla luciferase with a synthetic poly-A60 tail, which were individually in vitro transcribed, capped, and 2′-O-methylated (Figure 5—figure supplement 1A and Materials and methods). We elected to transcribe RNA both because plasmid-encoded transcripts can be heterogeneous and to mimic a scenario where one is delivering RNA to affect cell activity or for therapeutic intervention.

18All tested 5′ leader sequences modulate protein production by the luciferase reporter in concordance with the observed polysome association by TrIP-seq, when transfected into HEK 293T cells (Figure 5A,B). The three 5′ leaders from EEF1B2 alter luciferase production in a stepwise manner by roughly a factor of 20. The TrIP-seq profiles from NAE1 are the most distinct among those tested, and the two luciferase constructs differ in output by two orders of magnitude. Surprisingly, the two 5′ leaders from RICTOR differ by only nine nucleotides, yet still exhibit differential protein production, possibly due to altered local RNA secondary structure near the 5′ cap (Figure 5—figure supplement 1B). Both 5′ leaders from SRSF5 contain two uORFs, but they are close (SRSF5-002) or far (SRSF5-005) from the start codon (Figure 5—figure supplement 1A), suggesting reinitiation following uORF translation may be impacting luciferase production. Not only do these data provide strong evidence that TrIP-seq data can be used to predictably tune the output of heterologous mRNAs using isoform-specific untranslated regions, it additionally validates that the polysome profiles observed are connected to translational output.

19We additionally selected eight 3′ UTRs from four sets of gene-linked isoforms to test their ability to control translation, and find good agreement with the TrIP-seq polysome abundance for two out of four 3′ UTR pairs (Figure 5C,D). In all cases, replacing the 3′ UTRs leads to decreased protein production compared to the short 3′ UTR in the control RNA. The paired 3′ UTRs of NAB1 and RICTOR each differ by two kilobases, and a factor of 35 and 44 in luciferase output, respectively, indicating that 3′ UTRs can strongly modulate protein production. However, the two other paired 3′ UTRs, despite also differing by two kilobases each, are not distinguishable at the protein production level in cells. The two tested CCNE2 isoforms differ in both the 5′ and 3′ UTRs, so it is possible that the isoform-level translational control is occurring via the 5′ UTR, and there is a small but not significant difference between the two NDC1 3′ UTRs, suggesting regulation of this 3′ UTR may be subtle. We also find a positive relationship between the average number of ribosomes on each transcript isoform in TrIP-seq and the luciferase fold change (Figure 5E). In sum, we show that elements found to control translation at the isoform level using TrIP-seq can be grafted on to heterologous coding sequences to control translation over a range of two orders of magnitude in human cells.

20As different cell types contain different macromolecules, they may also translate messages with the same regulatory features differently. We, therefore, tested the panel of gene-specific untranslated regions fused with luciferase (Figure 5) in four additional cell lines: A549 (lung carcinoma), K-562 (chronic myelogenous leukemia), MCF-7 (breast adenocarcinoma), and Hep G2 (hepatocellular carcinoma). We additionally tested HEK 293T cells at a shorter timepoint of 2 hr, to explore the role of RNA stability in the observed luciferase output.

21Remarkably, the luciferase production from all tested 5′ leader reporters was qualitatively similar in all six conditions (Figure 6A). However, the difference between the two RICTOR 5′ leader sequences is decreased. In contrast, the 3′ UTR reporters show considerably more variability (Figure 6B). Calculation of the coefficient of variation of all reporters across all six conditions highlights the greater variability found in 3′ UTR reporters (Figure 6—figure supplement 1A). We additionally compared HEK 293T cells at 2-hr and 8-hr post-transfection to ascertain differences in RNA stability. Most 5′ leader reporters are similar between these two timepoints (Figure 6A), but 3′ UTR reporters show larger changes. Specifically, NAB1-001 and RICTOR-001 generate much more luciferase at the 2-hr timepoint, suggesting the long UTRs of these genes may promote RNA degradation leading to decreased protein at eighteen hours.

22Thus, translational control by the 5′ leader sequence is robust to changes in macromolecular composition across the cell lines tested, perhaps because the abundance of most translation initiation factors is high, which will buffer against small changes in their expression level. However, factors interacting with 3′ UTRs, specifically miRNAs, can vary considerably between cell types, and this variability may be the cause of the diverse behavior of the 3′ UTR reporters across cell lines. Therefore, at least for the panel tested here, transcript 5′ leaders confer robust control of protein production across cell types, while transcript 3′ UTRs may be better suited to tune the production of a protein to a particular cell type.

23In this work, we showed that the dynamic range of transcript-isoform-specific translational control spans at least two orders of magnitude in human cells (Figure 5), indicating that it is crucial to take isoform-level effects into account when assessing translation in organisms with extensive alternative transcript processing (Figure 7A). Globally, we find that alternative 3′ UTRs broadly influence translation, with longer isoforms of the same gene associated with lower protein production (Figure 3B). We further show that regulatory regions are sufficient to control the translation of unrelated coding sequences, enabling predictable tuning of translational output of arbitrary genes (Figures 5 and 6B). We demonstrate that translational control conferred by a panel of 5′ leaders is more robust across cell types than 3′ UTRs, suggesting predictable control of protein production. Future work focused on measuring isoform-specific translation in different cell types will yield cell type-specific regulatory sequences, which could then be used to engineer cell type-specific translation of mRNA, as in derivation of pluripotent stem cells, or mRNA therapeutics.

24Several observations indicate that TrIP-seq is a faithful measure of translation in cells. First, TrIP-seq and ribosome profiling are highly correlated at the gene level (Figure 1F) and both ribosome profiling and TrIP-seq correlate with mass spectrometric measurements of protein abundance (, Figure 1—figure supplement 1F). Second, protein-coding genes are enriched in high ribosomal fractions (Figure 1E), while lncRNAs are found primarily in the low ribosomal fractions (Figure 1E and Figure 1—figure supplement 1D). Third, a systematic investigation of translational stalling in mouse embryonic stem cells found this phenomenon to primarily be transient and not lead to ribosome accumulation on transcripts, although in specific circumstances it can be more widespread. Fourth, known highly- and lowly-translated transcript isoforms are enriched in high- and low-polysome fractions in TrIP-seq data (Figure 1C). Fifth, reads derived from stalled polysomes would artificially inflate the apparent translatability, however, TrIP-seq primarily underestimates rather than overestimates protein abundance (Figure 1—figure supplement 1F). Sixth, cryptic 'pseudo-polysomes' induced by miRNA complexes would preferentially enhance polysome association for long 3′ UTR transcript isoforms, which is the opposite of what is observed (Figure 3B). Lastly, transcript isoform changes found to lead to differential ribosome association using TrIP-seq are sufficient to modulate the protein output from a reporter RNA in a manner predicted by the TrIP-seq data (Figure 5), indicating that transcript-specific ribosome association is a correlate of protein output. Therefore, we conclude that TrIP-seq is a measure of transcript-specific translation.

25What are the mechanistic foundations of the observed transcript-specific translation? Each of the different untranslated regions tested by reporter assays operate using different mechanisms. For example, the three EEF1B2 transcript isoforms all share the same uORF, but alternative splicing changes the distance between the uORF and the start codon, which can influence the efficiency of downstream initiation. Globally, we found a large difference between the 3′ UTR lengths of poorly- versus well-translated isoforms by querying the prevalence of features likely to influence translation (Figures 3 and 5). There are at least three possibilities for how 3′ UTRs could influence translation. First, miRNA-mediated regulation targets the 3′ UTR. Second, numerous RNA binding proteins target the 3′ UTR through AU-rich elements or other specific binding sites and are known to influence translation, mRNA decay and even protein localization. Lastly, it is possible that increasing the length of the 3′ UTR influences the impact of the mRNA closed loop on translation initiation or reinitiation. Targeted work on specific 3′ UTR isoform sets could test the different possibilities for 3′ UTR-mediated translational control. We additionally demonstrated isoform-specific translational control for select 5′ leaders (Figures 1 and 5) as was shown genome-wide in yeast, where 5′ leaders are short enough to directly sequence, and find that globally uORFs can both up- and down-regulate translation (Figure 3). However, we find a surprisingly minor global dependence of other 5′ leader features tested on translation. The mechanisms of translational control of transcript isoforms therefore follow some general trends but are likely highly idiosyncratic; these data empower investigations into the mechanisms underlying transcript-specific translational control for thousands of human genes.

26We predicted the consequences of transcript isoforms observed during human embryonic development on protein production (Figure 4). We show that even small changes to transcripts can dramatically affect protein production (Figure 5), so changes at the isoform level (Figure 4) can considerably affect protein abundances even if the gene-level RNA expression remains similar. The application of TrIP-seq to different cell types should yield cell type-specific translation enhancer and repressor elements, contributing to the understanding of cell type-specific translational control, and augmenting our ability to precisely engineer translation in complex systems. If simpler sample preparation is desired, it is likely sufficient to sequence the cytoplasmic and monosome fractions, as well as pooled two-four ribosome and five-eight+ ribosome fractions, as these cluster together (Figure 1B). Accumulation of TrIP-seq data in additional cell types may enable building a holistic model of translational control with isoform resolution in human cells.

27We showed that the translation of an arbitrary gene, Renilla luciferase, can be controlled in a manner predicted by TrIP-seq data over two orders of magnitude (Figure 5) and in different cell types (Figure 6). These data can thus be used to select regulatory regions to control translation, without redesign of the coding sequence of the message (Figure 7B). It may prove superior to use 3′ UTRs to design cell type specific translation, as the repertoire of miRNAs and RNA-binding proteins that may affect translation through the 3′ UTR vary between cell types (Figure 6B). Indeed, in Caenorhabditis elegans, 3′ UTRs are sufficient to specify germline-specific expression of the attached ORF. Even without engineering of unnatural mRNAs, it may be possible to use antisense oligonucleotides to direct splicing of poorly- or well-translated isoforms to adjust protein expression in situ, for example to downregulate the oncogenic NAE1 gene [Figure 5A, ]), which is currently being targeted by small molecules. We anticipate the ability to tune the translational output of mRNA will facilitate research and therapeutic uses of designed and endogenous mRNA molecules.

28Two independently passaged, biological replicate 15 cm dishes of HEK 293T cells obtained from the University of California, Berkeley cell culture facility were grown to ~70% confluency in DMEM + 10% FBS. The cell line was authenticated by DDC Medical (Fairfield, OH) and were verified to be free of mycoplasma contamination. Cells were actively growing when harvested. The media was aspirated and replaced by PBS + 100 μg/ml cycloheximide and incubated at 37°C for 10 min. We chose cycloheximide because it induces rapid protein synthesis arrest, has been successfully used in ribosome profiling of HEK 293T cells, and 100 μg/ml is ~100 times higher than the concentration required to inhibit protein synthesis in reticulocytes. Each dish was then placed on ice, media aspirated, and replaced by ice cold PBS + 100 μg/ml cycloheximide. Cells were scraped, pelleted at 16,000× g for 30 s, and re-suspended in three pellet-volumes ice cold hypotonic lysis buffer (10 mM HEPES pH 7.9, 1.5 mM MgCl2, 10 mM KCl, 0.5 mM DTT, 1% Triton X-100 and 100 μg/ml cycloheximide). After 10 min, cells were lysed on ice by ten strokes through a 26-gauge needle and nuclei were pelleted at 1,500× g for 5 min. Lysate from ~15 million cells (one dish) was layered on top of triplicate 10–50% (w/v) sucrose gradients (20 mM HEPES:KOH pH 7.6, 100 mM KCl, 5 mM MgCl2, 1 mM DTT and 100 μg/ml cycloheximide) made using a Biocomp Instruments (Canada) gradient master. Gradients were centrifuged for 2 hr at 36,000 RPM in a SW-41 rotor, punctured, and manually peak fractionated using real-time A260 monitoring with a Brandel (Gaithersburg, MD) gradient fractionator and ISCO (Lincoln, NE) UA-6 detector.

29RNA was extracted from pooled technical triplicate sucrose gradient fractions by ethanol precipitation followed by acid phenol:chloroform extraction. Direct phenol:chloroform extraction was precluded by phase inversion in high sucrose fractions. RNA was then DNase treated, acid phenol:chloroform extracted, and ethanol precipitated. Cytoplasmic RNA was TRIzol extracted (Life Technologies, Grand Island, NY) and ethanol precipitated. Total RNA integrity was verified using a BioAnalyzer (Agilent, Santa Clara, CA). Ribosomal RNA was then depleted using Ribo-Zero (Illumina, San Diego, CA) and biological duplicate sequencing libraries were generated using the TruSeq RNA Sample Prep v2 kit (Illumina) without the poly-A selection steps. An equal mass (100 ng) of rRNA-depleted RNA was used as input to each individual library preparation. Libraries were verified using a BioAnalyzer and quantified using a Qubit (Life Technologies) prior to pooling for sequencing. Library insert sizes were typically ~150 ± ~25 bp. Pooled libraries were 75-bp paired-end sequenced on an Illumina HiSeq 2500 and runs of the same library in different lanes or flowcells were merged.

30Adapters were trimmed using Cutadapt v1.5 followed by subtractive alignments against the repeatmasker (RMSK) database (retrieved from UCSC on 2/11/2015) and abundant sequences from the Illumina iGenomes project (e.g. ribosomal RNA and the mitochondrial chromosome) using Bowtie2 v2.2.4. Unaligned reads were then aligned to the Ensembl release 75 transcriptome using Tophat v2.0.13 with parameters "-r 5 --mate-std-dev 50 -g 100 --report-secondary-alignments". Mapping percentages for each pipeline stage are presented in Figure 1—figure supplement 3. Transcript isoform level abundances were calculated using Cuffquant v2.2.1, normalized between samples with Cuffnorm v2.2.1, and transcripts per million (TPM) values were calculated according to:

31Biological duplicate datasets were processed independently. Quantified abundances are in Figure 1—source data 1 and Figure 1—source data 2. TrIP-seq plots are available for all Ensembl GRCh37 isoforms at http://meru.qb3.berkeley.edu/tripseq.

32Cuffnorm counts for each replicate were subjected to a variance stabilizing transformation (VST) using the DESeq2 R package to correct for heteroscedasticity; the VST approaches log2 for large counts but compresses low counts to suppress Poisson noise. Variance stabilized counts were averaged between replicates, filtered such that the mean across all nine samples was greater than one (roughly translating to 100 reads), and then mean subtracted to generate relative expression values. Inter-row distance was computed using Spearman’s rank correlation, and hierarchical clustering was performed using the fastcluster R package with Ward’s agglomeration method. The resulting dendrogram was split at increasing heights into subtrees until clusters contained similar overall trends, generating the large clusters presented in Figure 2. The samples were clustered similarly except with Euclidean distance and complete agglomeration. The R packages magrittr, dendextend, and ggplot2 were used to generate figures shown and are available through CRAN. Transcript isoforms in each cluster are in Figure 2—source data 1.

33The intersample clustering in Figure 1B was performed by hierarchically clustering the VST-transformed isoform-level counts as above between samples. The statistical significance of the resulting clustering was analyzed by measuring the Jaccard distance between clusterings of these data after subjecting the data to subsampling by bootstrap, jittering (adding random noise to each point), or replacing random points by noise using the R package fpc. The mean Jaccard distances of 100 such subsamplings as well as the average are presented in Figure Figure 1—figure supplement 1G S1G; a Jaccard distance ≥.75 is considered a stable cluster by the R package fpc.

34Isoform-specific 5′ and 3′ untranslated regions were amplified from anchored oligo-dT primed cDNA libraries from HEK 293T cells and Gibson cloned into a vector based on pUC57 containing Renilla luciferase and a synthetic polyA60 tail (pA60; ). Gibson cloning was performed such that untranslated regions were precisely cloned next to the ATG but contained two guanosines as the 5′-most nucleotides for T7 transcription for 5′ leaders, or the stop codon used by the isoform, and the six nucleotides CTGCAG at the 3′ end of the 3′ UTR immediately preceding the polyA60 tail for 3′ UTRs. The differences between cloned isoforms are shown in Figure 5—figure supplement 1, which was generated using IGV. The entirety of all untranslated regions was verified using dideoxy sequencing. Transcription templates were generated by PCR using Phusion polymerase (NEB, Ipswich, MA), size verified by gel electrophoresis, and gel purified. Transcription was performed using T7 polymerase with 1 μg template in a buffer containing 7.5 mM each NTP, 1 μg pyrophosphatase (Roche, Pleasanton, CA), 30 mM DTT, 35 mM MgCl2, 2 mM spermidine, 0.01% Triton X-100 and 30 mM Tris pH 8.1 for four hours at 37°C followed by DNase treatment with RQ1 RNase-free DNase (Promega, Madison, WI) for 30 min. Transcription products were purified by ethanol precipitation and a Zymo (Irvine, CA) Clean & Concentrator column, followed by simultaneous capping using Vaccinia capping enzyme (NEB) and 2′-O-methylation (NEB). Capped products were purified using a Zymo Clean & Concentrator column, followed by size verification on glyoxylated samples (Ambion, Foster City, CA) using an agarose gel, and full-length 3′ UTR-containing RNAs were gel purified from an agarose gel (Zymo).

35The concentration of RNAs containing 5′ untranslated regions was determined using A260, which was normalized by the intensity of the full-length product on an agarose gel, and then by the molar ratio of the construct to empty pA60. Molar-adjusted amounts of each RNA relative to 100–200 ng of pA60 were transfected into three technical triplicate wells of ~50% confluent cells in a 96-well plate using the TransIT-mRNA reagent (Mirus, Madison, WI). The concentration of RNAs containing 3′ UTRs was determined using a Qubit RNA HS assay (Life Technologies), normalized for the molar ratio of the construct to empty pA60, and molar-adjusted amounts of RNA relative to 7 ng of pA60 were transfected using TransIT-mRNA. A reduced amount of the 3′ UTR RNAs was used due to decreased yield of the longest 3′ UTRs. HEK 293T cells were grown in DMEM + 10% FBS, Hep G2 cells were grown in EMEM + 10% FBS, MCF7 cells were grown in DMEM:F12 + 10% FBS, A549 cells were grown in F12-K media + 10% FBS, and K-562 cells were grown in RPMI media + 10% FBS. Cells were harvested after ~18 hr (Figure 5) or 2 hr (Figure 6) and Renilla luminescence was measured (Promega).

36The Ensembl release 75 annotation set from the Illumina iGenomes project was first decomposed into 5′ leader, start codon, CDS, 3′ UTR, and whole-transcript regions. Length, GC-content, and number of exons were computed directly. Cytoplasmic expression and the median expression from the 80S through the eight+ ribosome fraction were extracted from TrIP-seq data. Transcript halflife data were derived from HeLa cell measurements. The structure was computed using RNALfold from the ViennaRNA package in a 75-nt window. Codon usage statistics were downloaded from http://www.kazusa.or.jp/codon/ on 11/20/2014 and the minimum codon frequency is the average of codon usage across a five-codon window. The fraction of AU-elements is calculated as the percentage of the 3′ UTR that is of repeating A or U nucleotides for more than 5nt in a row. TargetScan 6.2 scores were downloaded from http://targetscan.org and parsed for the properties indicated, and an identical comparison was performed after filtering the miRNA list by those expressed in HEK 293T cells. All feature tabulation was performed using custom Python programs, which are available through GitHub at https://github.com/stephenfloor/tripseq-analysis.

37Isoforms belonging to the same gene present in different clusters were compiled, yielding gene-linked isoforms. For example, if the gene A has isoform 001 in cluster one and isoforms 002 and 003 in cluster two, features of isoform 001 are added to the cluster one set and features of 002 and 003 are added to the cluster two set and 001–002 as well as 001–003 would be gene-linked isoforms. Features in each set were then compared for statistical significance using the Mann-Whitney U test and visualized as empirical cumulative distribution functions (Figure 3—figure supplement 1A). The effect size was then computed between distributions using Cliff’s d, which is a nonparametric, dimensionless measure of the distance between distributions. Cliff’s d is a measure of the number of times that a point xi in one distribution is greater than all points xj in the second distribution, or

38where # denotes the number of times, the two distributions are of sizes and, and and are items of the two sets. Cliff’s d is also related to the Mann-Whitney U statistic, by

39Cliff’s d and the boostrap confidence intervals shown in Figure 3 were computed using the R package orddom, which is available through CRAN.

40Ribosome profiling was performed in HEK 293T cells. These data were downloaded from the NCBI (GEO: GSE65778) and reprocessed by subtractive Bowtie 2 alignment to rRNA and RMSK sequences, and then mapped onto the human transcriptome using Tophat, as for the TrIP-seq data. Aligned reads were quantified at the gene and isoform level using Cuffquant onto the Ensembl release 75 transcriptome and normalized with Cuffnorm as above, to facilitate direct comparison to the TrIP-seq data. The number of reads derived from polysomes for TrIP-seq data was computed as the sum of read counts for fractions containing two to eight+ ribosomes multiplied by the number of ribosomes in each fraction:

41where is the number of reads in the th polysome fraction averaged for both biological replicates. The factor of 24 is applied to the last fraction based on a ceiling of ~40 ribosomes per transcript which then leads to an average of 24 ribosomes in the eighth peak (40 + 8 / 2 = 24). The number of ribosomes in the final peak is not directly measured by the sucrose gradients used Figure 1—figure supplement 1A).

42RNA sequencing datasets of human preimplantation embryos were downloaded from the NCBI (GEO: GSE36552). Reads were converted to FASTQ using fastq-dump (NCBI) and then processed as for TrIP-seq data by adapter trimming with Cutadapt, subtractive alignment to the RMSK and abundant sequences, and transcriptome alignment with Tophat to Ensembl release 75. Transcript abundances were calculated using Cuffquant and normalized with Cuffnorm. Data from individual cells were processed independently and averaged for each stage to generate the data shown in Figure 4 and Figure 4—figure supplement 1. The data were reprocessed for consistency to facilitate direct abundance comparisons.

43Cytoplasmic lysate was fractionated using a sucrose gradient as for the RNA-seq libraries and RNA was extracted from each fraction. Libraries of cDNA were then synthesized from these fractions using random primers (Applied Biosystems) from an equal amount of RNA as measured using a Qubit (Life Technologies; RNA HS assay). Template RNA was removed using RNase H treatment (NEB) and cDNAs were purified over an oligo cleanup column (Zymo). PCR was then performed for 25 (ACTB) or 30 (EEF1B2, ATF4) cycles with gene-specific primers (below) using Taq Titanium (Clontech) using an equal amount of cDNA input into each PCR, as measured using a Qubit (ssDNA assay). Reactions were then run on a 1% agarose gel and stained with SYBR-Gold (Life Technologies). Gels are representative of three biological replicates.

44Primers used for RT-PCR are:

45ACTB-forward AGAGCTACGAGCTGCCTGAC

46ACTB-reverse AGCACTGTGTTGGCGTACAG

47EEF1B2-forward TTCCCGTCATCTTCGGGAGCCGT

48EEF1B2-reverse CTTTTCAGGTCTCCGAAACCCATGG

49ATF4-forward GGCTCTGCAGCGGCAACCCC

50ATF4-reverse CGACTGGTCGAAGGGGGACA

51RNA was extracted from polysome fractions as for TrIP-seq. qRT-PCR was performed using the SuperScript III Platinum SYBR Green One-Step kit (Life Technologies) with 1ng input RNA in a 20 ul reaction volume. CT values were converted to fold changes over the cytoplasmic abundance and plotted.

52Primers used for qRT-PCR are:

53EEF1B2-001-fwd AGCCGTGGAGCTCTCGGATA

54EEF1B2-001-rev AGCTCTTGTCCGCCAGGTAA

55EEF1B2-003-fwd TCCAAACGCAACGAAAGGTCC

56EEF1B2-003-rev CGCCGGACCGAAGGTTAAAG

57EEF1B2-201-fwd AGCCGTGGAGCGTGGG

58EEF1B2-201-rev TCGGCTGTATCCGAGAGCTG

59SRSF5-002-fwd GACCCCGTCCGGTAGGAAGTACTAGCC

60SRSF5-002-rev CAATATCTCTTATCCGTCCATATCCC

61SRSF5-005-fwd GGTGAGTGGCTCACTTTGAGGGCAAG

62SRSF5-005-rev CGAATCAACTGCGCTCATTAGACGC

63The DAVID server was used to calculate gene ontology (GO) terms for biological processes (BP) associated with the individual clusters. Transcripts associated with each cluster were input into DAVID and GO BP terms with a Benjamini-corrected p-value of <0.05 were tabulated.

64Raw sequencing reads for all samples are available through the NCBI via the GEO Accession ID GSE69352.