Cellular stress alters 3′UTR landscape through alternative polyadenylation and isoform-specific degradation

We were interested in APA profiles in cells under stress and during recovery from stress. We treated mouse NIH3T3 cells with 250 μM sodium arsenite (NaAsO₂), a commonly used strong stressor, for 1 h, and let cells recover for 4, 8, 12, or 24 h after removal of arsenic stress (AS) (Fig. 1a). As expected, while cell death was not detected, cells after stress displayed significantly slowed growth (Supplementary Fig. 1). In line with the integrated stress response mechanism¹⁴, the level of eIF2α phosphorylation increased by ~7-fold after 1 h of AS and went back to the baseline level after 4 h of recovery (RC, Fig. 1b), indicating successful stress response and recovery. Consistently, SGs were detected by immunocytochemistry using anti-PABPC1 antibody after 1 h of AS, and disappeared after 4 h of recovery (Fig. 1c).

To examine APA, we extracted total cellular RNA from non-treated (NT) control cells, cells with 1 h of AS, and those at different time points of recovery (4–24 h), and subjected them to 3′ region extraction and deep sequencing (3′READS), a method specialized for interrogation of 3′ ends of poly(A) + transcripts^21,40. To simplify APA analysis, we focused on the two most abundant 3′UTR APA isoforms as determined by PAS reads (see Methods for details). Based on their relative locations in the 3′UTR, they were named proximal PAS (pPAS) and distal PAS (dPAS) isoforms, respectively (illustrated in Fig. 1d, top). We found that the number of genes showing shortened 3′UTRs in AS-treated cells (pPAS isoform abundance relative to that of dPAS isoform increased) was 2.3-fold greater than the number of genes showing the opposite trend (247 vs. 108, Fig. 1d, bottom), indicating that AS leads to general 3′UTR shortening in the cell. The global 3′UTR APA changes were also measured by the relative expression difference (RED) score, which was the difference between two samples in ratio (log2) of dPAS isoform abundance (reads per million, RPM) to pPAS isoform abundance. A negative RED value indicates 3′UTR shortening and a positive value 3′UTR lengthening. To robustly compare samples with different sequencing depths, we used the Global Analysis of Alternative Polyadenylation (GAAP) method we previously developed⁴¹, which employs random sampling of reads from comparing samples to derive a median RED and its standard deviation (see Methods for details). Using GAAP, the median RED of cells with AS was −0.11 (Fig. 1e, red bar).

Interestingly, 3′UTR shortening became more conspicuous after 4 h of recovery (median RED = −0.27, P = 1.5 × 10⁻¹¹ vs. 1 h AS, Wilcoxon test, Fig. 1e). The median RED score gradually increased afterwards (median RED = −0.14 and −0.12 for 8 and 12 h of recoveries, respectively, Fig. 1e) and was similar to that of NT cells after 24 h (median RED = 0.03, Fig. 1e).

We found that the extent of 3′UTR shortening of a gene, as measured by the RED score, was a function of the distance between pPAS and dPAS, or size of alternative UTR (aUTR, illustrated in Fig. 1d, top), in both AS and RC cells (Fig. 1f). That is the longer the aUTR, the greater the 3′UTR shortening. The gene group with the longest aUTRs (bin 5 in Fig. 1f) displayed significantly greater shortening than the group with the shortest aUTRs (bin 1 in Fig. 1f) in both AS and RC cells (P= 2.7 × 10⁻²⁰ and 5.1 × 10⁻²⁹, respectively, Wilcoxon test).

An example gene Nmt1 (encoding N-myristoyltransferase 1) is shown in Fig. 1g, which expressed two 3′UTR APA isoforms with 3′UTR lengths of 337 nt and 2.8 kb, respectively. Using reverse transcription-quantitative PCR (RT-qPCR) and primer pairs targeting a region common to both pPAS and dPAS isoforms as well as the aUTR, we validated 3′UTR shortening of Nmt1 in cells with AS and in recovery, as well as those of Calm1, Purb, and Timp2 (Fig. 1h). Notably, 3′UTR shortening could also be detected in cells treated with oxidative stress inducer H₂O₂ (Supplementary Figs. 2a and 2b) or with a lower level of AS (25 µM) (Supplementary Figs. 2c and 2d), although the timing and the degree of shortening were different.

Shortening of 3′UTRs can be caused by increased usage of the PASs of short 3′UTR isoforms or by enhanced degradation of long 3′UTR isoforms, or both. We reasoned that analysis of APA profiles with newly made RNAs would reveal PAS usage changes⁴². To this end, we added 4-thiouridine (4sU) to cell culture media 1 h before harvesting cells in NT, AS (1 h), or RC (4 h) conditions (Fig. 2a), separated 4sU-labeled RNAs from non-labeled, or flow-through (Ft) RNAs and subjected both RNA pools to 3′READS analysis (Supplementary Fig. 3a, see Methods for details). As such, newly made RNAs and pre-existing RNAs were enriched in the 4sU and Ft fractions, respectively.

Similar to the result with steady-state RNAs, 4sU-labeled RNAs displayed 3′UTR shortening in both RC and AS cells as compared with NT cells (Fig. 2b). However, the difference in 3′UTR shortening between the two conditions, although still significant, was much smaller in newly made RNAs (ΔMedian RED = 0.03; P = 4.4 × 10⁻¹⁰, Wilcoxon test, Fig. 2b), as compared to the difference using steady state RNAs (ΔMedian RED = 0.16, Fig. 1e).

We next compared RED scores between 4sU and Ft fractions. In control cells, 3′UTRs were longer in newly made RNAs than in pre-existing RNAs (median RED = 0.10, Fig. 2c), consistent with the notion that short 3′UTR isoforms are generally more stable than long 3′UTR isoforms⁴³. A similar trend was observed with RC cells, but to a greater degree (median RED = 0.18, Fig. 2c). By contrast, newly made RNAs in AS cells had shorter 3′UTRs than pre-existing RNAs (median RED = −0.06, Fig. 2c). These results indicate that 3′UTR shortening in cells with 1 h AS treatment is largely due to newly made RNAs, supporting the notion that AS leads to preferential usage of proximal PASs. While this also takes place in RC cells (Fig. 2b), the much greater shortening of 3′UTRs in those cells appears to be attributable to the pre-existing pool of RNA (Fig. 2c), suggesting differential regulation of mRNA stability between 3′UTR isoforms.

To specifically address how mRNA stability contributed to 3′UTR APA profiles, we calculated for each isoform the ratio of its abundance in the Ft fraction to that in the 4sU fraction, or log2(Ft/4sU). After adjustments for the number of uridines in each transcript to remove its influence on the abundance of 4sU-labeled RNAs (Supplementary Fig. 3b), the log2(Ft/4sU) values in control cells correlated well with mRNA half-lives previously measured in NIH3T3 cells by Spies et al.⁴³ (r = 0.52, Pearson correlation, Supplementary Fig. 3c), supporting the suitability of using log2(Ft/4sU) values to reflect mRNA stability.

To examine the influence of 3′UTR size on stability, we divided genes into five bins based on their aUTR sizes. We found that long 3′UTR isoforms in control cells were generally less stable than short 3′UTR isoforms, as indicated by the negative median Δlog2(Ft/4sU) values in all gene groups (Fig. 2d, gray line). The difference in stability between short and long 3′UTR isoforms was a function of the aUTR size (P = 5.5 × 10⁻⁵ comparing genes in bins 1 and 5, Wilcoxon test, Fig. 2d, gray line), indicating a negative role of 3′UTR size in mRNA stability. Significantly, the difference between long and short 3′UTR isoforms became greater in cells during recovery (blue vs. gray lines, Fig. 2d), indicating enhanced 3′UTR size-based mRNA decay. Consistently, the difference in stability between gene bins 1 and 5 was much greater in RC cells (P = 9.4 × 10⁻⁹, Wilcoxon test, Fig. 2d, blue line). By contrast, in AS cells Δlog2(Ft/4sU) values were positive across all bins and were not correlated with aUTR sizes (red line in Fig. 2d), indicating general stabilization of distal PAS isoforms during stress. The stability difference between short and long 3′UTR isoforms in control cells and in cells after stress was validated by RT-qPCR analysis of mRNA decays of APA isoforms of Nmt1 and Purb (Fig. 2e).

We next asked whether some sequence motifs in 3′UTRs were associated with mRNA decay during recovery from stress. To this end, we identified two groups of genes (Fig. 2f, top), i.e., genes whose long 3′UTR isoforms were less stable than short 3′UTR isoforms (group 1) and those showing the opposite trend (group 2). We then compared tetramer frequencies in the aUTRs of these two gene groups. Interestingly, U-rich tetramers were significantly enriched in the aUTRs of group 1 genes, whereas C-rich elements in those of group 2 genes (Fig. 2f, bottom), indicating involvement of sequence motifs in mRNA decay in cells after stress.

Together, our data indicate that long 3′UTR isoforms are generally less stable than short 3′UTR isoforms in control cells; this trend is blunted in cells under stress but is accentuated during recovery from stress. In addition, U-rich motifs play a prominent role in mRNA decay during recovery from stress.

Sequence motif analysis of unstable mRNAs suggested that U-rich motif-binding RBPs might be responsible for the enhanced 3′UTR-based mRNA decay in cells after stress. We next focused on TIA1 because of its activity in binding U-rich motifs and role in recruiting mRNAs into SGs during stress¹⁶. As reported, TIA1 was mostly in nucleus in control cells and became enriched in SGs during stress (Supplementary Fig. 4a). Four hours of recovery after AS restored the TIA1 subcellular distribution (Supplementary Figure 4a). Notably, TIA1 protein levels were unchanged across these conditions (Supplementary Fig. 4b).

To examine interactions between mRNAs and TIA1, we carried out ribonucleoprotein immunoprecipitation (RIP) using an anti-TIA1 antibody (Supplementary Fig. 4b) to isolate TIA1-bound mRNAs from cytoplasm of AS-treated or control cells, followed by 3′READS (illustrated in Fig. 3a). As such, this approach, named 3′READS + RIP, could distinguish APA isoforms in TIA1 binding.

Using input cytoplasmic RNA as a reference, we calculated TIA1 binding efficiency, as calculated by log2(TIA1/Input), for each transcript with a defined PAS (see Methods for details). Two biological replicates were generated, which were well correlated (r = 0.65 and 0.72 for NT and AS samples, respectively, Pearson correlation, Supplementary Fig. 4c). In addition, log2(TIA1/Input) values were well correlated between AS-treated and control cells (r= 0.75, Pearson correlation, Fig. 3b), indicating that TIA1 binding efficiency with its target RNAs generally does not change by AS. Also notable is that the TIA1 binding efficiency measured by 3′READS + RIP is generally correlated with the number of TIA1 binding sites in 3′UTRs as measured by a recent iCLIP study using B cells⁴⁴ (r = 0.34, Supplementary Fig. 4d), indicating similarities in TIA1 binding across cell types and supporting the suitability of using 3′READS + RIP for quantitative analysis of TIA1 binding.

Interestingly, we found that TIA1-bound transcripts tended to have significantly longer 3′UTRs in both control and AS cells, as indicated by the positive median Δlog2(TIA1/Input) values measuring the difference between dPAS and pPAS isoforms (Fig. 3c) and that genes whose long 3′UTR isoform had a higher log2(TIA1/Input) value than short isoform markedly outnumbered genes with the opposite trend (4.7- and 5.1-fold for NT and AS cells, Supplementary Fig. 4e). Consistently, the Δlog2(TIA1/Input) value between long and short 3′UTR isoforms was a function of the aUTR size in both NT and AS cells (Fig. 3d). In line with these results, the log2(TIA1/Input) values of all transcripts correlated with their 3′UTR size (r = 0.62, Pearson correlation, Supplementary Fig. 4f), and 3′UTR size was a better determinant than other transcript features, such as 5′UTR and CDS sizes, for interaction with TIA1 based on a linear regression model (Supplementary Fig. 4g). These results indicate that the efficiency of TIA1 binding to an mRNA positively correlates with its 3′UTR size.

We next examined sequence motifs in aUTRs of two gene groups, i.e., genes whose long 3′UTR isoforms had a higher TIA1 binding efficiency than short 3′UTR isoforms (group 1, dPAS » pPAS in Fig. 3e), and those with the opposite trend (group 2, dPAS « pPAS in Fig. 3e). Consistent with the known TIA1 binding properties⁴⁴, U-rich motifs, especially UUUU, were significantly enriched in the aUTRs of group 1 genes (Fig. 3e, bottom). An example gene Tmem248 is shown in Fig. 3f, which expressed a 174-nt, short 3′UTR isoform and a 2,211-nt, long 3′UTR isoform. A prominent TIA1 binding site was located in its aUTR, as identified by the iCLIP study⁴⁴. Consistently, the long 3′UTR isoform was markedly enriched by 3′READS + RIP as compared to the short 3′UTR isoform (Fig. 3f). In addition, using the iCLIP data and all mouse PASs in the PolyA_DB database⁴⁵, we found that about 80% of TIA1 binding sites in 3′UTRs were located in aUTRs (Fig. 3g). Taken together, these results indicate that long 3′UTR isoforms preferentially bind TIA1 as compared to short 3′UTR isoforms through U-rich motifs in alternative 3′UTR sequences.

We found that TIA1 binding efficiency of a transcript, as measured by log2(TIA1/Input), was negatively correlated with its mRNA stability as measured by log2(Ft/4sU)(r = −0.56, Pearson correlation, Fig. 4a), suggesting a connection between TIA1 binding and mRNA decay. Importantly, the transcripts with both strong TIA1 binding and a short half-life tended to have significantly longer 3′UTRs than those with both weak TIA1 binding and a long half-life (Fig. 4b). In addition, transcripts with strong TIA1 binding in AS cells tended to have significantly decreased mRNA stability during recovery (Δlog2(Ft/4sU), RC vs. AS) than transcripts with weak TIA1 binding (P = 6.0 × 10⁻²⁵, K–S test, Fig. 4c). More importantly, long 3′UTR isoforms which had higher TIA1 binding efficacies in AS cells than short 3′UTR isoforms were significantly less stable in RC cells than the latter (Fig. 4d). Together, these data indicate that TIA1 binding correlates with enhanced mRNA decay during recovery from stress.

A recent study identified mRNAs in the SG core through isolation of RNAs bound to G3BP1, another resident RBP in SGs⁴⁶. We thus asked whether TIA1 binding efficiency was related to enrichment in the SG core. To this end, we obtained the SG-enrichment data of human transcripts expressed in U-2 OS cells that were generated by Khong et al.⁴⁶ and compared them to TIA1 binding efficiencies of their mouse orthologs calculated in this study (see Methods for details). As shown in Fig. 4e, SG-enrichment was generally correlated with TIA binding in AS cells (r = 0.58, Pearson correlation), supporting the role of TIA1 in recruiting mRNAs to SGs and indicating that the recruitment is conserved between human and mouse orthologs. Importantly, U-rich motifs were significantly enriched in the 3′UTRs of transcripts with both high SG-enrichment scores and high TIA1 binding efficiencies (light blue dots in Fig. 4e). U-rich motifs were also enriched, albeit with a lower significance, for transcripts with high TIA1 binding efficiencies but modest SG-enrichment scores (red dots in Fig. 4e). By contrast, the 3′UTRs of transcripts with high SG-enrichment scores but modest TIA1 binding efficiencies (dark blue dots in Fig. 4e) had no enrichment of U-rich motifs. This result indicates that U-rich motifs facilitate TIA1 binding leading to efficient SG recruitment, and some transcripts can be recruited to SGs through TIA1-independent pathways. Taken together, as summarized in Fig. 4f, our data reveal a TIA1-based SG recruitment mechanism in cells under stress, which leads to enhanced mRNA decay during recovery.

We next reasoned that if short 3′UTR isoforms had a higher mRNA stability than long isoforms in cells after stress, the genes with shortened 3′UTRs in response to stress should have a higher transcript abundance after stress than other genes. To this end, we measured RNA abundances in NT, AS, and RC cells using a strand-specific RNA-seq method (Fig. 5a, see Methods for details). We divided genes into three groups, namely, genes with significantly shortened 3′UTRs through APA by stress (group 1), genes expressing 3′UTR APA isoforms but without significant 3′UTR shortening (group 2), and genes showing only one 3′UTR isoform (no APA, group 3). We defined group 1 genes as those showing 3′UTR shortening in both the 4sU fractions of AS and RC cells (blue dots in Fig. 5b).

All three gene groups had a median log2(ratio of gene expression change) of 0 after 1 h of AS, indicating no global expression changes. By contrast, after 4 h of recovery, group 1 genes showed a global upregulation of expression (median log2(ratio) = 0.26), significantly higher than groups 2 and 3 genes (median log2(ratio) = 0.04 and −0.14, respectively; P = 3.6 × 10⁻⁶ and 7.6 × 10⁻¹⁵, respectively, K–S test, Fig. 5c). This result indicates that AS-induced 3′UTR shortening can indeed preserve transcripts during recovery from stress. Because group 3 genes, but not group 2 genes, showed general downregulation of gene expression, this result also indicates that having APA sites can help genes maintain expression after stress. Interestingly, Gene Ontology (GO) analysis indicated that genes with functions in cell differentiation, cellular component assembly, signal transduction, and cell proliferation were enriched for group 1 genes (Table 1), indicating that these genes are more likely to use APA to evade mRNA degradation in stressed cells. By contrast, no GO terms were found to be significantly associated with genes showing lengthened 3′UTRs.

To further validate APA changes by stress and its impact on gene expression, we resorted to reporter assays using the pRiG vectors we previously constructed²⁶. Due to APA, the pRiG-AD vector could express two isoforms: a long isoform encoding both red fluorescent protein (RFP) and enhanced green fluorescent protein (EGFP), and a short isoform encoding RFP only (Fig. 6a). Twelve hours after transfection, we subjected cells transfected with pRiG-AD to AS for 1 h and measured abundances of two APA isoforms by Northern blot 12 h later (Fig. 6b). Consistent with the global trend of cellular genes, we found that the ratio of short isoform abundance to that of the long isoform increased by 1.4-fold in AS-treated cells (Fig. 6c).

In agreement with the Northern blot data, analysis of cells with fluorescence-activated cell sorting (FACS, Supplementary Fig. 5) showed a 2.8-fold increase of the red to green fluorescence ratio (log2(red/green)) in cells transfected with pRiG-AD and treated with AS, as compared to transfected cells without treatment (Fig. 6d, middle). By contrast, a much milder increase (1.1-fold, Fig. 6d, top) was observed in cells transfected with the pRiG vector, which expressed the long isoform only due to absence of a PAS between the RFP and EGFP sequences⁴⁷.

We also transfected cells with a mixture of pRiG and pRiG-SV40 vectors, the latter of which predominantly expresses the short APA isoform owing to the presence of a strong SV40 early PAS between the RFP and EGFP sequences⁴⁷. As such, while both short and long isoforms were produced in cells transfected with pRiG + pRiG-SV40, the isoforms were expressed from two separate plasmids. We adjusted the ratio of the two plasmids to ensure similar amounts of short and long isoforms expressed in the transfected cells to cells transfected with pRiG-AD. The pRiG + pRiG-SV40 mixture yielded an increase of 1.6-fold in log2(red/green) in stressed vs. control cells (Fig. 6d, bottom), substantially lower than that of pRiG-AD. This result indicates that enhanced usage of proximal PAS contributes to the increased log2(red/green) value in pRiG-AD-transfected cells after stress. On the other hand, the difference between cells transfected with pRiG and those with pRiG + pRiG-SV40 indicates that short 3′UTR isoforms can have higher protein expression potentials than long 3′UTR isoforms in cells after stress.

To further validate the effect of 3′UTR APA on protein expression, we cloned full 3′UTR sequences of Nmt1, Timp2, and Dnajb1 into the psiCHECK2 vector (Fig. 6e). All three genes displayed AS-induced 3′UTR shortening and the full 3′UTRs cloned all contained APA sites. Using dual luciferase assays, we found that after stress constructs containing these 3′UTRs produced ~40–50% higher luciferase activities than the control construct lacking these 3′UTRs (Fig. 6f). In the case of Nmt1, we also found that its full 3′UTR, which had the ability of APA, led to 1.2-fold higher luciferase activities than the short 3′UTR (cUTR) in stressed cells (Fig. 6g), indicating the impact of APA on protein expression in response to stress.

Together, reporter assays confirmed that stress enhances proximal PAS usage in the cell, leading to increased expression of short 3′UTR isoforms that have higher protein expression potentials than long 3′UTR isoforms.

The GO analysis result indicated that genes with functions in cell proliferation and differentiation tended to have their 3′UTRs shortened by stress (Table 1). We next asked how stress impacts 3′UTR lengths in cells in different proliferation and differentiation states, which have different APA programs²⁶. To this end, we analyzed APA by 3’READS in proliferating C2C12 myoblast (MB) cells and differentiated myotube (MT) cells under normal, stressed (250 μM NaAsO₂ for 1 h), and recovery (4 h after AS) conditions (Fig. 7a).

As in NIH3T3 cells, AS elicited significant 3′UTR shortening in both MB and MT cells to a similar extent after 1 h of stress, with median RED = −0.14 and −0.15 for MB and MT cells, respectively (Fig. 7b). Also similar to NIH3T3 cells, both MB and MT cells showed more significant 3′UTR shortening after 4 h of recovery, with MB cells showing a greater extent of 3′UTR shortening (median RED = −0.48) than MT cells (median RED = −0.31) (Fig. 7b). Note that these 3′UTR size changes are greater in scale than the 3′UTR size changes taking place in cell differentiation (median RED = −0.19, MB vs. MT, Supplementary Figs. 6a and 6b), highlighting the substantial impact of stress on 3′UTR landscape in the cell. A modest correlation in 3′UTR changes could be discerned between MT and MB cells during recovery from stress (r = 0.37, Pearson correlation, Supplementary Fig. 6c).

We next asked whether stress-regulated 3′UTR APA events were related to those modulated in cell proliferation/differentiation (P/D). To this end, we compared 3′UTR APA changes between MT and MB cells (representing P/D) with those between RC and NT samples of MT and MB cells. A mild correlation could be discerned between APA changes in MT vs. MB and those by stress (r = 0.13, Pearson correlation, Fig. 7c). We then divided genes into four groups based on 3′UTR shortening in the two conditions: both P/D and stress (Both), P/D only, Stress only, and neither P/D nor stress (Other) (Fig. 7c). Interestingly, genes in the Both group tended to have longer aUTRs (median = 1759 nt, Fig. 7c) as compared to genes in the Stress only group (median = 788) or in the P/D only group (median = 1103). By contrast, the median aUTR size was the shortest (397 nt) for genes in the Other group. This result indicates that long aUTR may confer regulability under both P/D and stress conditions.

We found that genes in the Both and Stress only groups showed significant 3′UTR shortening in stressed NIH3T3 cells as well, as compared to genes in the Other group (P = 7.8 × 10⁻¹⁵ and = 5.3 × 10⁻¹⁵, respectively, Fig. 7d). By contrast, no difference could be discerned between genes in the P/D only group and those in the Other group (Fig. 7d), highlighting the distinction between stress and P/D in 3′UTR length regulation.

Using RT-qPCR, we validated several genes in different groups (Fig. 7e). Nmt1 showed 3′UTR regulation in both P/D and stress. Timp2 also showed 3′UTR regulation in P/D, but its regulation by stress was restricted to MB cells only. Rpl22 and Fam49b showed substantial 3′UTR shortening in P/D but modest shortening by stress, whereas 3′UTR shortening of Dnajb1 and Hspa4l only occurred in stressed cells. Together, these results indicate that stress-induced 3′UTR shortening takes place in both proliferating and differentiated cells, with proliferating cells showing a greater extent of changes. While 3′UTR length controls by P/D and stress are largely distinct, some genes have their 3′UTRs regulated in both conditions.

NIH3T3 cells from ATCC were cultured in high glucose Dulbecco’s modified Eagle's medium (DMEM) supplemented with 10% calf serum. C2C12 cells from ATCC were cultured in DMEM supplemented with 10% fetal bovine serum. C2C12 cell differentiation was induced by changing culture media to DMEM supplemented with 2% horse serum when cells were at >95% confluency. All media contained 100 IU/ml penicillin and 100 μg/ml streptomycin. All cells were incubated at 37 °C with 5% CO₂. To induce stress, cells were incubated in culture medium containing 250 µM or 25 µM of sodium arsenite (NaASO₂), or 1 mM of H₂O₂. For stress recovery, stressed cells were washed twice with phoshate-buffered saline (PBS) to remove NaASO₂ or H₂O_2, followed by culturing in regular media for a period of time before harvest. Cell viability assay and cell counting were performed using Vi-CELL Cell Viability Analyzer (Beckman Coulter).

pRiG, pRiG-SV40 and pRiG-AD were described in ref. ²⁶. psiCHECK2-3′UTR constructs were made by inserting PCR-amplified mouse 3′UTR sequences into the psiCHECK2 plasmid (Promega) linearized with XhoI and NotI. The PCR primer sequences are listed in Supplementary Table 3.

Twelve hours after transfection with reporter plasmids using Lipofectamine 3000 (Thermo Fisher), NIH3T3 cells were treated with 250 µM sodium arsenite for 1 h. Cells were then washed twice with PBS and cultured in regular medium for another 12 h. Cells lysate was prepared, and dual luciferase assay was performed using Dual-Luciferase Reporter Assay System (Promega) following the manufacturer’s protocol. Biological triplicates were performed.

Cells grown on glass coverslips were fixed using 4% paraformaldehyde in PBS for 10 min at room temperature and permeabilized with PBS containing 0.1% Triton X-100 at room temperature for 10 min. After blocking with 1% bovine serum albumin in PBST (PBS + 0.1% Tween 20) for 30 min, SGs were immunostained with rabbit anti-PABPC1 (Abcam, ab21060, 1:500 in PBST) or goat anti-TIA1 (Santa Cruz Biotechnology, sc-1751, 1:500 in PBST) for 60 min. FITC-conjugated goat anti-rabbit antibody (Jackson ImmunoResearch, 111-095-144, 1:500 in PBST) or Texas red-conjugated donkey anti-goat (Jackson ImmunoResearch, 705-075-147, 1:500 in PBST) was applied at RT for 60 min. Cells were mounted on glass slides in SlowFade Gold Antifade reagent (Thermo Fisher) with DAPI. Fluorescence images were collected using the EVOS FL Auto Cell Imaging System (Thermo Fisher).

The protocol for isolation of TIA1-associated RNA was adapted from the RIP protocol described in Keene et al.⁵⁹ Briefly, NIH3T3 cells with or without 1 h of AS were trypsinized, washed twice with pre-chilled PBS, re-suspended and incubated in lysis buffer (20 mM HEPES-KOH pH 7.5, 15 mM MgCl₂, 80 mM KCl, 1% Triton X-100, 2 mM DTT, 100 μg/ml cycloheximide, 1× SIGMAFAST protease inhibitor cocktail, and 200 U/ml SuperaseIn) on ice for 10 min. During the incubation, cells were gently sheared three times through a pre-chilled 26-gauge needle. The cell lysate was centrifuged at 700 × g at 4 °C for 5 min to pellet nuclei. The supernatant was transferred to a new tube and centrifuged at 14,000 × g at 4 °C for 5 min. About 10% of the supernatant was saved for cytoplasmic RNA extraction, while the rest of the supernatant was diluted five times with NT2 buffer (50 mM Tris, pH 7.4, 150 mM NaCl, 1 mM MgCl₂, 0.05% NP-40) supplemented with 100 U/ml SuperaseIn. The diluted lysate was incubated with goat anti-TIA1 IgG (Santa Cruz, sc-1751) conjugated to Dynabeads Protein G at 4 °C for 4 h. After washing the beads five times with NT2 buffer supplemented with 1 M urea, RNA was extracted from the beads and the input cytoplasmic fraction using Trizol, followed by 3′READS analysis. Alternatively, proteins were eluted from the beads following the manufacturer’s protocol and analyzed by western blot analysis using rabbit anti-TIA1 (Proteintech, 12133-2-AP). Uncropped blots can be found in Supplementary Fig. 7d.

The 3′READS procedure (3′READS+ version) was described in ref. ⁴⁰. Briefly, Poly(A)+ RNA in 0.1 or 1 µg of input RNA was captured using oligo(dT)₂₅ magnetic beads (NEB) and fragmented on the beads using RNase III (NEB). After washing away unbound RNA fragments, poly(A)+ RNA fragments were eluted from the beads and precipitated with ethanol, followed by ligation to heat-denatured 5′ adapter (5′-CCUUGGCACCCGAGAAUUCCANNNN) with T4 RNA ligase 1 (NEB). The ligation products were captured by biotin-T₁₅-(+TT)₅ (Exiqon) bound to Dynabeads MyOne Streptavidin C1 (Thermo Fisher). After washing, RNA fragments on the beads were incubated with RNase H to remove bulk of the poly(A) tail and then eluted from the beads. After precipitation with ethanol, RNA fragments were ligated to a 5′ adenylated 3′ adapter (5′-rApp/NNNGATCGTCGGACTGTAGAACTCTGAAC/3ddC (Bioo Scientific) with T4 RNA ligase 2 (truncated KQ version, NEB). The ligation products were then precipitated and reverse transcribed using M-MLV reverse transcriptase (Promega), followed by PCR amplification using Phusion high-fidelity DNA polymerase (NEB) and bar-coded PCR primers for 15–18 cycles. PCR products were size-selected twice with AMPure XP beads (Beckman Coulter). The size and quantity of the libraries were examined on an Agilent Bioanalyzer. Libraries were sequenced on an Illumina NextSeq 500 (1 × 75 bases).

Cells were cultured in medium supplemented with 50 µM of 4-thiouridine (4sU; Sigma) for 1 h before harvest. Total RNA was extracted using TRIzol. Newly made (4sU-labeled) and pre-existing RNA populations were fractionated following the protocol described in ref. ⁴². Briefly, 100 µg of total RNA was biotinylated using biotin-HPDP (1 µg/µl in DMF; Thermo Fisher Scientific), and then extracted with chloroform three times and precipitated with ethanol. The biotinylated RNA was captured by Streptavidin C1 Dynabeads. The beads were washed six times, and the unbound, flow-through (Ft) RNA was collected. Bound, biotinylated RNA was eluted by DTT. All RNA was precipitated with ethanol and then used for RT-qPCR or 3′READS experiments.

Poly(A)+ RNA in 5 µg of total RNA was purified twice using oligo(dT)₂₅ magnetic beads (NEB) following the manufacturer’s protocol. RNA was then fragmented by partial alkaline hydrolysis (94 °C for 2 min), 3′ dephosphorylated with shrimp alkaline phosphatase (NEB), and 5′ phosphorylated with T4 PNK (NEB). The RNA fragments were sequentially ligated to a 5′ adenylated 3′ blocked 3′ adapter (5′-Aden/NNNNNNGATCGTCGGACTGTAGAACTCTGAAC/3ddC, where N is a random nucleotide) with truncated RNA ligase 2 (KQ, 200 U/µl) and to a 5′ adaptor (5′-CCUUGGCACCCGAGAAUUCCA) with T4 RNA ligase I. After reverse transcription of ligation products, cDNA amplification by PCR, and library purification by AMPure beads (Beckman Coulter), libraries were sequenced on an Illumina NextSeq 500 (1 × 75 bases).

Twelve hours after transfection with reporter plasmids using Lipofectamine 3000, NIH3T3 cells were treated with 250 µM sodium arsenite for 1 h. Cells were then washed twice with PBS and cultured in regular medium for another 12 h. Cells were released from culture dishes by Trypsin-EDTA. Green and red fluorescent signals were read at 530 and 585 nm, respectively, in a BD FACScalibur system (BD Biosciences). Un-transfected cells were used to determine background fluorescence.

Northern blot analysis was performed using the NorthernMax^TM kit reagents (Thermo Fisher). Briefly, 12 h after plasmid transfection using Lipofectamine 3000, NIH3T3 cells were treated with or without 250 µM sodium arsenite for 1 h. After washing with PBS twice, cells were cultured in regular medium for another 12 h. Total RNA was extracted using TRIzol reagent (Thermo Fisher). A total of 5 µg of RNA was loaded onto a 1% formaldehyde-agarose gel for electrophoresis. RNA was then transferred onto a nylon membrane (Thermo Fisher). For hybridization, DIG-labeled DNA probes were generated by PCR using RFP-specific primers and the PCR DIG Probe Synthesis Kit (Sigma). The DIG-labeled probes were used for overnight hybridization with the NorthernMax^TM Kit. Immunological detection of the DIG-labeled probes was performed using the DIG Northern Starter Kit (Sigma). Signals on the blot were detected by a ChemiDoc Touch Imaging System (Bio-Rad). Images were quantified using the Fiji program⁶⁰. Uncropped blots can be found in Supplementary Fig. 7b.

cDNA was synthesized using the M-MLV reverse transcriptase (Promega), oligo(dT) primer, and 2 µg of RNA that was pre-treated with Turbo DNase (Thermo Fisher). cDNA was then mixed with primers and Hot Start Taq-based Luna qPCR master mix (NEB). PCR was run on an StepOne Plus Real Time PCR system (Thermo Fisher). For analysis of APA isoforms, primers were designed to amplify a region common to both pPAS and dPAS isoforms (common), and a region between the two PASs (aUTR). Relative expression of dPAS isoform to pPAS isoform was calculated by log2(aUTR/common). Two-tailed t-test (without assuming equal variance) was performed to calculate P-values. qPCR primers are listed in Supplementary Table 1.

Protein concentration was measured using the DC Protein Assay (Bio-Rad). A total of 20 µg of protein per sample was resolved by sodium dodecyl sulfate polyacrylamide gel electrophoresis, followed by transfer to PVDF membrane for immunoblotting. Signals on the blot after adding the ECL reagent (Bio-Rad) were detected by the ChemiDoc Touch Imaging System (Bio-Rad). Quantification was carried out in using the Fiji program⁶⁰. Pre-stained protein molecular weight marker (Thermo Scientific, #26612) was used for protein size calculation. Uncropped blots can be found in Supplementary Figs. 7a, 7c and 7d.

3′READS data were analyzed using methods described previously⁴⁰. Briefly, the 5′ adapter sequence was first removed using Cutadapt⁶¹. Reads with short inserts (<23 nt) were discarded. The retained reads were mapped to the mouse genome (mm9) using bowtie2 (local mode)⁶². The three or six random nucleotides at the 5′ end (from the 3′ adapter) were removed before mapping using the setting “−5 3” or “−5 6” in bowtie2. Reads with a mapping quality score (MAPQ) ≥10 were kept for further analysis. Reads with ≥2 non-genomic 5′Ts after alignment were called PAS reads. PASs within 24 nt from each other were clustered. PAS read counts mapped to genes were normalized by the median ratio method in DESeq⁶³. Only APA isoforms with read count greater than five in the pair of compared samples were used. The two most abundant APA isoforms (based on PAS reads) in the 3′UTR of the 3′-most exon were selected for 3′UTR APA analysis. Significant APA events were those with a relative abundance change >5% and P-value <0.05 (Fisher’s exact test). RED was calculated as the difference in log2(ratio) of abundances of two PAS isoforms (dPAS vs. pPAS) between two samples. To statistically assess the significance of overall APA differences between samples with different sequencing depths, we used the GAAP method we previously developed⁴¹. Briefly, one million PAS reads from each sample were randomly sampled through bootstrapping, and a median RED was calculated each time using all sampled genes. Sampling was carried out 20 times, and the mean of median RED and standard deviation were calculated.

For each transcript with a defined PAS, its expression value (RPM) in the flow-through (Ft) sample was divided by the RPM value in the 4sU sample. The resulting value, log2(Ft/4sU), was adjusted to control for the number of uracils (Us) in each transcript because of its effect on labeling and capture efficiencies of 4sU-labeled RNA. A linear regression model for log2(number of Us) vs. log2(Ft/4sU) was constructed by a python function. The adjusted log2(Ft/4sU) was calculated using the model.

TIA1 binding efficiency for each PAS-defined transcript was calculated by log2 (RPM in the TIA1 RIP sample/RPM in the input sample). The averaged binding efficiency of each transcript based on two replicates was used for analysis. Linear regression models using transcript features, including 5′UTR size, CDS size, 3′UTR size, number of exons and exon density, againt log2(TIA1/input) were constructed with the python program ‘scikit-learn’.

SG-enrichment information of human transcripts was obtained from the supplementary file in ref. ⁴⁶. TIA1 binding efficiency of the most abundant transcript of each mouse gene was compared to the SG-enrichment score of its human orthologue based on NCBI HomoloGene database⁶⁴.

Raw reads from RNA-seq were first trimmed for adapter sequences by Cutadapt⁶¹and then mapped to the mm9 genome assembly using STAR (v2.5.2)⁶⁵ with default parameters. The read count of the coding sequence (CDS) of each gene was summed and then normalized by the median ratio method in DESeq⁶³. Only genes with more than five reads in all samples were used.

The GO terms and their associations with mouse genes were obtained from the Gene Ontology database⁶⁶. The Fisher’s exact test was used to derive P-values to indicate significance of association between a gene set and a GO term.

Genomic positions and corresponding read counts of TIA1 iCLIP sites were obtained from GSE93575↗ (ref. ⁴⁴). iCLIP sites were further clustered within a 6-nt or 30-nt window. iCLIP sites with fewer than 10 reads were not used in analysis.

Assignment of TIA1 binding sites to cUTRs and aUTRs was based on mouse mammal-conserved PASs in PolyA_DB 3 (ref. ⁴⁵).

Tetramer frequencies in aUTRs were calculated and compared between gene sets. P-values based on the Fisher’s exact test were used to indicate the significance of enrichment.

All scripts used for data processing and statistical analysis were written in Python, Perl, or R, and are available upon request.

Sequencing datasets generated in this study, including those by 3’READS and RNA-seq (listed in Supplementary Table 2), have been deposited into the GEO database under the accession number GSE101851↗ [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE101851↗]. All other data supporting the findings in this study are available from the corresponding author on reasonable request.

The authors declare no competing interests.

Sequencing datasets generated in this study, including those by 3’READS and RNA-seq (listed in Supplementary Table 2), have been deposited into the GEO database under the accession number GSE101851↗ [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE101851↗]. All other data supporting the findings in this study are available from the corresponding author on reasonable request.