Translation inhibition, but not disruption of the EJC-PYM1 interaction, results in widespread changes in EJC footprints

To investigate the relative contributions of PYM1 and the elongating ribosome to EJC disassembly in human cells, we measured EJC occupancy after either blocking the EJC-PYM1 interaction or translation elongation (Fig. 1a). To block EJC-PYM1 interaction, we took advantage of a previously reported mutant of MAGOH (E117R, MAGOH), which can still assemble into the EJC but does not interact with PYM1 due to mutation of a key side chain (Gehring et al., 2009, Supplementary Fig. 1a, b). We generated a HEK293 cell line that stably expresses FLAG-tagged MAGOH along with the endogenous MAGOH protein (Supplementary Fig. 1c). As expected, like stably expressed wild-type (WT) FLAG-MAGOH, FLAG-MAGOH interacts with EIF4A3 and RBM8A, but unlike FLAG-MAGOH, it shows no detectable interaction with PYM1 (Fig. 1b). Moreover, the immunoprecipitation (IP) of FLAG-tagged wild-type and mutant MAGOH proteins shows that the FLAG-MAGOH EJC core does not show any change in stability under increasing salt concentrations (Supplementary Fig. 1d) or in its ability to co-immunoprecipitate several key EJC co-factors tested (Supplementary Fig. 1e). Thus, the EJC containing FLAG-MAGOH is likely to maintain all its known activities except for its interaction with PYM1.

To analyze changes in EJC occupancy resulting from any defects in PYM1-mediated EJC disassembly, we performed RIPiT-Seq via tandem IP of FLAG-MAGOH (WT or E117R) and endogenous CASC3. We chose CASC3 for the second IP to enrich cytoplasmic/later-stage EJCs which are likely to be more susceptible to defects in EJC disassembly. To analyze the effect of translation elongation inhibition on EJC occupancy, both WT and the mutant MAGOH expressing cells were treated with cycloheximide (CHX) for 3 h prior to their harvesting for RIPiT-Seq (Supplementary Fig. 2a, b). Notably, 3 h of CHX treatment does not lower nuclear levels of EJC subunits or alter PYM1-MAGOH interaction (Supplementary Fig. 2c, d). The two RIPiT-Seq replicates from each condition were well-correlated for gene-level counts for uniquely mapping EJC footprint reads and clustered with their respective experimental condition in principal component analysis (Supplementary Fig. 2e-i). Therefore, individual replicates were combined for further analysis.

To quantify EJC occupancy for each protein-coding gene in the reference human transcriptome (Ensembl release 100, hg38), we estimated total EJC footprint signal across the entire length of a representative transcript (as per the MANE select annotation). As expected based on a previous report that translating ribosomes remove EJCs from mRNAs, footprints of FLAG-MAGOH and CASC3 containing EJCs (MAGOH EJC from hereon) show dramatic changes between CHX-treated and untreated cells with more than 1/4 of all transcripts showing significant changes in EJC footprints (Fig. 1c, 988 genes with significantly increased and 882 genes with significantly decreased EJC footprints as defined by padj < 0.05; Supplementary Data 1). In contrast, when MAGOH EJCs were compared to those containing FLAG-MAGOH and CASC3 (MAGOH EJC from hereon), overall EJC occupancy changes significantly for only 26 genes (Fig. 1d; Supplementary Data 2). Thus, translocating ribosomes have a profound effect on EJC occupancy in HEK293 cells whereas the effect of EJC-PYM1 interaction on the EJC landscape is more subtle (also see Supplementary Fig. 2j). This finding agrees with a recent study by Bensaude et al. estimating that 85% of EJCs are disassembled via the translation machinery.

We predicted that EJCs bound to efficiently translated, cytoplasmic mRNAs would be the most sensitive to the dramatic impact of translation elongation inhibition. Consistently, transcripts with significantly increased EJC footprints upon CHX treatment have the highest relative cytoplasmic levels (Fig. 1e; median cytoplasmic rank percentile = 71.95) whereas mRNAs that show a relative increase in EJC occupancy in untreated cells have the lowest cytoplasmic abundance (Fig. 1e; median cytoplasmic rank percentile = 52.28; p = 3.78 × 10). Thus, localization of mRNAs from the nucleus to the cytoplasm promotes translation-dependent EJC disassembly. Next, transcripts with higher EJC occupancy upon CHX treatment appear to be translated more efficiently as they have a significantly higher normalized ribosome occupancy (median logTE = -1.07) as compared to mRNAs with elevated EJC occupancy in untreated cells (Fig. 1f; median logTE = -1.27; p = 2.32 × 10). Finally, mRNAs with increased EJC occupancy upon CHX treatment have significantly shorter median half-lives (median log transformed t = 0.25) as compared to mRNAs enriched in EJC under normal conditions (Fig. 1g; median log transformed t = 1.55; p = 2.07 × 10). Presumably, translation inhibition also leads to mRNA stabilization, which in turn can contribute to the observed higher EJC signal on mRNAs upon CHX treatment. Gene Ontology (GO)-term enrichment analysis of gene sets with increased EJC footprint signal upon CHX treatment identifies keywords including nucleus, mRNA processing and mRNA splicing (Supplementary Fig. 2k), functional groups of genes that encode unstable mRNAs. Thus, cytoplasmic mRNA metabolic processes are critical in shaping EJC occupancy on mRNAs.

Paradoxically, a subset of mRNAs appears to have decreased EJC occupancy in CHX-treated cells (or increased occupancy in untreated cells). Multiple factors could be responsible for this effect. Among this set of genes, the most prominently enriched GO terms are for cytoplasmic ribosomes and ribosomal proteins (Supplementary Fig. 2k), consistent with earlier reports that translation inhibition lowers EJC occupancy on this set of mRNAs, perhaps due to their increased translation under these conditions. Additionally, some of this effect could also be an artifact of normalization by DESeq2. Within the CHX-treated versus untreated comparison, a proportional increase in the number of reads from transcripts that are highly translated would result in a decrease in the proportion of reads from transcripts with unchanged EJC footprints, which may consequently appear down regulated. Thus, transcripts with an apparent decrease in EJC footprints upon CHX treatment can have biological underpinnings or appear so due to technical issues.

We next determined if disruption of EJC-PYM1 interaction influences EJC occupancy on distinct groups of mRNAs. Although FLAG-MAGOH EJC footprints significantly change only a few transcripts, we compared properties of genes in the top and bottom quartiles based on fold change in footprints of the mutant versus the wild-type EJC. While the EJC occupancy upon loss of EJC-PYM1 interaction shows no correlation with relative nucleocytoplasmic mRNA distribution or translation efficiency (Fig. 1h, i), the quartile of transcripts with the most MAGOH EJC footprint signal have lower half-lives as compared to those in the bottom quartile (Fig. 1j; median log transformed t= -0.49 and 0.62 for the top and bottom quartiles, respectively; p = 2.45 × 10). The genes with the highest mutant EJC signal are enriched in GO term keywords related to transcriptional regulation, DNA binding and ribonucleoprotein, classes that are characterized by unstable mRNAs (Supplementary Fig. 2k). On the other hand, genes with comparatively lower signal for the mutant EJCs (and hence more association with the wild-type EJC) are enriched in membrane-related GO terms (Supplementary Fig. 2k). Similar membrane-related GO terms are also observed in genes with reduced EJC signal upon CHX treatment. Thus, even though the disruption of EJC-PYM1 interaction generally has only a modest effect on gene-level EJC occupancy, distinct functionally related classes of genes appear to be more dependent on this interaction.

We next investigated the effects of EJC-PYM1 interaction on EJC occupancy in a translation-dependent and -independent manner. A direct comparison of fold changes in EJC footprints associated with translation inhibition and disruption of EJC-PYM1 interaction for each gene shows a very weak correlation between the two conditions (R = 0.12; Fig. 2a). Conversely, a strong correlation is observed between fold changes of the mutant and wild-type EJCs upon CHX treatment (R = 0.86; Fig. 2b). Thus, the loss of PYM1-MAGOH interaction does not alter the inherent ability of the EJC to disassemble during translation. Next, we hypothesized that the effect of PYM1 on EJC occupancy on long non-coding RNAs (lncRNAs) would reflect its translation-independent function. We find that, as compared to the wild-type EJC, the MAGOH EJC shows a significantly higher occupancy on lncRNAs (Fig. 2c, p = 2.61 × 10). To rule out any potential effect of translation on lncRNA EJC occupancy, we compared the differences in MAGOH EJC versus MAGOH EJC binding on the top quartile of nuclear-localized lncRNAs (based on their relative nucleocytoplasmic distribution) to all other lncRNAs. Notably, the top quartile contains lncRNAs with >95 percentile rank for their nuclear localization. We observe a significantly higher MAGOH EJC occupancy on the nuclear lncRNAs (Fig. 2d, p = 0.0263), thereby suggesting that EJC-PYM1 interaction influences EJC binding in a translation-independent manner even though it has no apparent effect on translation-dependent changes in global EJC occupancy.

To investigate if the inhibition of EJC-PYM1 interaction or translation elongation affects EJC positioning on exons, EJC footprint reads from all four conditions (MAGOH and MAGOH EJCs with and without CHX treatment) were plotted relative to exon-exon junctions (Fig. 3a and Supplementary Fig. 3a). Notably, in these meta-exon distributions, all genes were equally weighted, regardless of their expression levels, to prevent bias toward highly expressed genes. Also, the area under the curve for each condition remains the same, so that the shape of the read distribution, and hence relative EJC positioning along exons, can be compared across the conditions. As expected, EJCs from all conditions show a major peak at the canonical EJC position, ~24 nucleotides (nt) upstream of exon 3' ends (Fig. 3a and Supplementary Fig. 3a). Both MAGOH and MAGOH EJCs show a similar footprint distribution from untreated and CHX-treated cells suggesting that, even though CHX dramatically changes the overall EJC occupancy at the transcript level, the relative EJC position on exons remains unchanged. Interestingly, contrary to the expectations if PYM1 was a major disassembly factor, MAGOH EJC footprints are reduced at the canonical EJC position as compared to MAGOH EJCs (Fig. 3a). In turn, the mutant EJC footprints show a modest increase in regions away from the -24 position (Fig. 3a). Previous studies have shown that EJCs can be detected at such non-canonical positions albeit at a lower level. Our observation here suggests that the PYM1 interaction-deficient EJC is more predisposed to persist in non-canonical regions than its wild-type counterpart. As we observed no differences in the number of novel splice junctions detected in the RIPiT-Seq datasets from the wild-type and mutant MAGOH expressing cells (Supplementary Fig. 3b), the observed non-canonical EJC binding is unlikely to be due to novel or altered splicing events in the annotated exonic regions.

To further examine the binding of the mutant and wild-type EJCs at non-canonical positions of each spliced protein-coding MANE select transcript, we separately counted MAGOH and MAGOH EJC footprint reads in (i) canonical regions, (ii) non-canonical regions of all exons with a downstream intron, and (iii) last exons. For non-canonical regions of the first to the penultimate exon, 100 nt regions from each exon's 5'-end and 150 nt immediately upstream of exon 3' -ends were excluded to minimize the contribution of canonical EJC footprints to non-canonical signal (see Methods for details). Also excluded from the analysis were any last exons that overlapped with known introns (based on a list compiled by Kovalak et al.). Differential expression analysis between the MAGOH and MAGOH EJC footprint signal in these three regions showed that as compared to the wild-type EJC, the mutant EJC was most enriched in non-canonical regions (Fig. 3b). The mutant EJC binding was also higher in the last exons as compared to the canonical regions. This increased propensity of the mutant EJC to bind to non-canonical regions is also supported by the enrichment of MAGOH EJC over the wild-type EJC on genes with fewer than average number of exons (i.e., lower intron density and hence longer non-canonical regions). Strikingly, we observe a progressive increase in mutant over wild-type EJC enrichment as exon count in genes decreases, particularly from five to two (Fig. 3c). These observations suggest that the loss of PYM1 interaction allows the EJC to persist longer or bind more often to non-canonical locations within long exonic stretches of RNA, thereby reducing the EJC binding specificity. Importantly, translation elongation inhibition does not cause such an increase in non-canonical EJC occupancy (Supplementary Fig. 3c, d).

Motivated by the above findings, we systematically examined the relationship between PYM1's influence on EJC occupancy and 34 different gene sequence and architectural features (Fig. 3d). We find that the top two features directly correlated with higher mutant EJC binding pertain to exon length (e.g., exon length geometric mean has Kendall's τ = 0.108, p = 1.35 × 10) (Fig. 3d). On the other hand, the top two features anti-correlated to the preferred mutant EJC binding are exon number (τ = -0.115, p = 2.02 × 10) and GC content (τ =-0.1, p = 9.01 × 10) (Fig. 3d). Next, we took an unbiased machine learning based approach to assess the relative importance of the 34 gene features in predicting changes in mutant versus wild-type EJC binding. We permuted the values of one feature at a time and calculated the loss in the model's ability to predict fold change in mutant versus wild-type EJC binding (see Methods for details). This analysis suggested that GC content is the most important feature for accurate prediction of the change in mutant EJC binding (Fig. 3e). Features pertaining to exon length and number are also observed among the top ten features in this analysis although their contribution to the fold change prediction is not significant. It is noteworthy that GC content was previously shown to be negatively correlated with exon number, which underscores their individual correlation to mutant EJC binding predicted by our model. Overall, these results confirm that disruption of PYM1-EJC interaction increases EJC occupancy on genes with fewer but longer exons that also have a higher GC content. While GC content also comes out as an even better predictor of increased EJC binding upon CHX treatment (Supplementary Fig. 4), both the mutant and wild-type EJC show similar relationships to various gene features in this condition, consistent with our earlier results that the loss of PYM1 binding does not affect translation-dependent EJC disassembly.

By definition, any EJC binding on unspliced RNAs will be considered non-canonical binding. Although the EJC is not expected to be deposited on RNAs transcribed from intronless genes, which are not acted upon by the spliceosome, our observation of increased MAGOH EJC binding in other intron lacking RNA segments like 3'UTRs prompted us to test for occupancy of the mutant and wild-type EJCs on single exon mRNAs. Importantly, we filtered single-exon genes to limit our analysis to only those that do not contain any known introns. Strikingly, when single-exon mRNAs are compared as a group to multi-exon mRNAs, a highly significant enrichment of MAGOH over MAGOH EJC footprints is observed on these unspliced mRNAs (Fig. 4a). Moreover, MAGOH EJC enrichment on single exon mRNAs is even more pronounced than spliced mRNAs containing two exons, which show the highest mutant EJC enrichment among the spliced mRNAs (Fig. 4a, inset). When MAGOH and MAGOH EJC footprint counts on single exon mRNAs are compared to other non-canonical regions (3'UTRs and sites away from -24 position in spliced exons), the highest enrichment of the mutant EJC footprints is again observed on single exon mRNAs (Fig. 4b). These findings suggest that under normal conditions EJC-PYM1 interaction limits non-canonical EJC binding, including on unspliced mRNAs.

We next sought to further characterize the unexpected EJC binding on single-exon mRNAs. At least 127 single exon mRNAs that lack any evidence of splicing show MAGOH EJC footprint signal in HEK293 cells, suggesting that EJC binding to unspliced transcripts is not uncommon. As expected, EJC footprint detection on these mRNAs is correlated with their RNA expression levels (R = 0.66, Fig. 4c). Interestingly, single-exon genes with similar RNA expression levels can exhibit large differences in EJC footprints as exemplified by two histone protein-encoding genes, H1-2 and H3C2 (Fig. 4c). EJC binding in non-canonical regions, including on unspliced mRNAs in RIPiT-seq datasets, could be due to an artificial assembly of EJC-like complexes on RNA post cell lysis. To test whether the non-canonical EJC binding on unspliced mRNAs can be observed in situ, we analyzed a published CLIP-Seq dataset that captures RNA segments directly UV-crosslinked to EIF4A3, the RNA-binding subunit of the EJC. As UV-crosslinking of the RNA-protein interaction surface happens prior to cell lysis, CLIP-Seq crosslinking sites provide more concrete evidence for RNA-protein interaction inside cells. We find that EJC RIPiT-seq read counts from HEK293 cells and EIF4A3 CLIP-Seq read counts from HeLa cells are modestly correlated on single exon mRNAs (Fig. 4d, R = 0.46). Further, EIF4A3 CLIP-Seq read distribution on H1-2 and H3C2 genes exhibits a similar trend (Fig. 4e) with H1-2 showing a much higher coverage of CLIP-Seq reads as compared to H3C2 despite similar expression levels of the two genes in HeLa cells. Importantly, the crosslinking of EIF4A3 with unspliced transcripts strengthens our conclusion that EJCs can bind to these unspliced RNAs. Although mechanisms that may dictate binding and specificity of non-canonical EJCs on unspliced mRNAs remain unknown, it is possible that splicing-independent association of EJC subunits to active transcription regions could contribute to such EJC binding.

Next, we evaluated if EJCs in non-canonical locations can have functional consequences. If such EJCs occur downstream of a stop codon, they can be expected to activate NMD. To test this idea, we identified transcripts without 3'UTR introns and grouped them into quartiles based on expression-normalized 3'UTR (non-canonical) EJC footprint signal. Interestingly, the quartile of genes with the highest 3'UTR non-canonical EJC footprint signal are the most upregulated in a dataset from HEK293 cells where NMD is inhibited via co-depletion of SMG6 and SMG7 (Fig. 5a). Indeed, we observe a significant and dose-dependent ability of 3'UTR non-canonical EJCs to regulate transcript levels via NMD. Notably, this effect is independent of 3'UTR lengths or transcript expression levels (Supplementary Fig. 5a, b). Further, a quantitative PCR-based assay showed that a subset of transcripts with 3'UTR non-canonical EJCs are unregulated in HEK293 cells treated with an inhibitor of the NMD factor SMG1, independently validating that these transcripts can indeed be regulated by NMD (Fig. 5b). Moreover, we took advantage of similar EJC occupancy and UPF3 double-knockout NMD inhibition datasets from our previous work in HCT116 cells, and detected robust, dose-dependent NMD activity of 3'UTR non-canonical EJCs in these cells as well (Supplementary Fig. 5c-e). Thus, we conclude that non-canonical EJCs in 3'UTR are capable of regulating transcript levels via NMD. We also surmise that such NMD activity of non-canonical 3'UTR EJCs is likely to spread to even more transcripts under PYM1-deficient conditions. We observe that the quartile of transcripts with lowest non-canonical 3'UTR EJC density in wild-type HEK293 cells shows the highest fold-increase in FLAG-MAGOH vs FLAG-MAGOH EJC (Supplementary Fig. 5f), which is also consistent with higher propensity for non-canonical EJC assembly in otherwise EJC poor regions (Fig. 4).

We next sought to investigate the role of PYM1 in maintaining gene expression in HEK293 cells. A 72-hour siRNA-treatment was used to deplete PYM1 to nearly a fifth of its normal levels in this cell line (Supplementary Fig. 6a). Such PYM1 depletion did not alter EJC-dependent splicing events (Supplementary Fig. 6b), or global translation as measured by puromycin incorporation into nascent polypeptides (Supplementary Fig. 6c, d). Therefore, we next considered the effect of PYM1 depletion on NMD, which was previously shown to be influenced by PYM1 using NMD reporter RNAs. We performed total RNA-Seq after ribosomal RNA depletion from the knockdown and control cells and quantified gene expression at the transcript level. To assess the effect of PYM1 on NMD of endogenous mRNAs, we focused on transcripts that contain premature termination codons (PTC), which are translation stops that occur >50 nucleotides from the downstream exon junction (Supplementary Data 3). As a control for the effect of PYM1 knockdown on other processes (e.g., transcription), these PTC+ transcripts were compared to other isoforms from the same genes that do not contain a PTC (PTC-). When we consider all PTC+ transcripts detected in HEK293 cells, a modest but significantly higher fold change is observed for PTC+ transcripts as compared to their PTC- counterparts (Supplementary Fig. 6e; Supplementary Data 4). When we limited the PTC+ group to transcripts that are > 1.5-fold upregulated in HEK293 cells double-depleted of SMG6 and SMG7 and therefore most susceptible to NMD, the effect of PYM1 depletion on this subset of PTC+ transcripts is more pronounced (Fig. 5c). Still, the global NMD defect upon PYM1 depletion (median log2FC of PTC+ transcripts = 0.05) is not as strong as that observed upon depletion of CASC3 (median log2FC of PTC+ transcripts = 0.38) or SMG6 and SMG7 co-depletion (median log2FC of PTC+ transcripts = 0.65) in HEK293 cells (Supplementary Fig. 6f). Among the select NMD targeted transcripts tested by quantitative RT-PCR based assay, only SRSF3 PTC+ mRNA shows a significant upregulation upon PYM1 knockdown whereas all tested NMD targets are robustly upregulated upon EIF4A3 depletion (Fig. 5d). Thus, NMD efficiency in PYM1 knockdown cells is reduced, but only modestly. Presumably, PYM1 deficiency leads to longer EJC lifetimes causing the complex to accumulate within untranslated RNPs and non-canonical locations thereby reducing canonical EJC deposition on newly spliced RNAs, and hence a minor NMD defect.

A previous study suggested that PYM1 overexpression can inhibit NMD of PTC-containing TCRβ and β-globin reporters by promoting EJC disassembly. To test if elevated PYM1 levels can cause a similar upregulation of endogenous NMD targets, we stably expressed in HEK293 cells N-terminal HA-tagged full-length PYM1 or its EJC interaction domain lacking deletion mutant, PYM1 ( ~ 5-fold and ~3-fold overexpression, respectively, Supplementary Fig. 6g). We confirmed that the deletion of the N-terminus of PYM1 disrupts the interaction of PYM1 with MAGOH (Supplementary Fig. 6h). However, the exogenous expression of full-length PYM1, like that of control PYM1, did not cause any upregulation of SMG6 and 7 sensitive PTC+ transcripts (Supplementary Fig. 6i, j) or all expressed PTC+ transcripts (Supplementary Fig. 6k, l). Thus, even ~5-fold increase in PYM1 levels does not appear to cause sufficient EJC destabilization to alter NMD of endogenous PTC+ mRNAs in HEK293 cells.

It is possible that reduced EJC specificity in PYM1 depleted cells also affects gene expression beyond NMD. To investigate this further, GO term enrichment analysis showed that transcripts up- and downregulated upon PYM1 depletion belong to distinct functional classes (Supplementary Fig. 7a). Strikingly, the top five GO terms enriched among the downregulated transcripts are related to membranes and membrane-linked processes (Supplementary Fig. 7a). The specific downregulation of membrane targeted transcripts in response to PYM1 depletion is also observed upon independent analysis of mRNAs translated at the ER (definition based on Jan et al., 2014), those encoding ER-targeting sequences (signal peptide or transmembrane domain), and the ones targeted to mitochondria and peroxisomes (Supplementary Fig. 7c; Supplementary Data 5-6).

Interestingly, a recent study has suggested a link between mRNA expression at the ER and gene architectural features such as exon number and exon length. mRNAs enriched in the ER were found to have higher than average exon count. On the other hand, mRNAs localized within an ER-linked membrane-less condensate known as TIS-granule (TG) have longer and fewer exons. These findings have intriguing parallels with our observations that membrane targeted mRNAs are downregulated upon PYM1 knockdown and those with long exons are enriched in footprints of PYM1 interaction-deficient EJCs (Fig. 3d). Thus, PYM1 (and hence EJC) function could be linked to gene architecture and mRNA expression. Indeed, we find that in PYM1 depleted HEK293 cells, mRNAs preferentially localized to TGs and the ER are significantly up- and downregulated, respectively (Fig. 6a, b). When we compare these two classes of mRNAs for binding to wild-type versus the PYM1-interaction deficient EJCs, we observe a similar although weaker trend (Fig. 6c). Notably, the increased binding of the mutant EJC on TG+ mRNAs is mainly in the non-canonical regions (Fig. 6d) including their 3´UTRs, which play a key role in regulation of this class of mRNAs.

The relationship between PYM1 function and gene architecture is further confirmed by direct correlations between changes in mRNA expression upon PYM1 knockdown and exon length-based features (Fig. 6e). Conversely, GC content, which is intimately linked to lower exon count, is the most anti-correlated feature to altered gene expression in PYM1 depleted cells (Fig. 6e). Of all the features tested in our machine learning based analysis, high GC content is also the most important predictor of the observed changes in gene expression upon PYM1 depletion (Fig. 6f), as it was for the increased binding to PYM1-deficient EJCs (Fig. 3d). In the PYM1 depleted cells, a small but significant downregulation of mRNAs in the top quartile based on exon number as compared to those in the bottom quartile remains when limiting to only non-ER+ transcripts (Supplementary Fig. 7d). Thus, the link between PYM1 function and gene architecture likely exists independently of ER/membrane localization.

We next sought to determine if the change in ER+ and TG+ mRNA levels upon PYM1 knockdown is due to transcriptional or post-transcriptional mechanism(s). To this end, an estimation of the change in stability of these classes of mRNAs showed that as compared to other mRNAs, ER+ mRNAs are destabilized whereas TG+ mRNAs are stabilized (Fig. 6g, median relative stability: ER + = -0.17, TG + = +0.46, Others = +0.04; p = 8.5 × 10: ER+ vs Others, p = 1.5 × 10: TG+ vs Others). This destabilization of ER+ mRNAs is possibly caused by a mechanism that is unrelated to general NMD (Supplementary Fig. 7e), the NBAS-dependent ER-linked NMD (Supplementary Fig. 7f, i), the ER-linked mRNA downregulatory activity of UPF1 long-loop isoform (Supplementary Fig. 7g, i), or thapsigargin-induced ER stress response (Supplementary Fig. 7h).

Recent biochemical screens for host proteins that interact with flavivirus capsid proteins have identified PYM1 as one of the key targets. This flavivirus capsid-PYM1 interaction reportedly aids in viral replication and causes an enrichment of EJC core proteins in a membrane fraction. Notably, flaviviruses remodel ER membranes to create special compartments where the viral RNA genome is replicated. Our observation that PYM1 knockdown alters stability of ER and TG associated mRNAs raises a hypothesis that, upon infection, the flavivirus capsids sequester free PYM1 molecules to reshape ER-linked mRNA expression, as observed upon PYM1 depletion. To test this prediction, we first confirmed that 2X Strep-tagged WNV and DENV capsid proteins co-immunoprecipitate with FLAG-tagged PYM1 (Fig. 7a and Supplementary Fig. 8a). In support of our hypothesis, transient WNV or DENV capsid protein expression in HEK293 cells led to upregulation of TG+ transcripts as compared to ER+ mRNAs (Fig. 7b and Supplementary Fig. 8b-d). Notably, these trends are similar to although weaker than those observed in PYM1 depleted cells (Fig. 6a). Interestingly, capsid protein expression also caused mildly elevated levels of PTC+ mRNAs as compared to their PTC- counterparts (Fig. 7c and Supplementary Fig. 8e). To examine if flavivirus infection also leads to similar changes in expression of TG+ and ER+ mRNAs, we compared expression of these groups of mRNAs in the available RNA-seq datasets from human hepatoma cells (Huh7) infected with ZIKV or DENV and human neuronal stem cells (hNSC) infected with ZIKV. In agreement with our experiments with transient capsid expression, all three flavivirus-infected samples displayed significant and robust upregulation of TG+ mRNAs as compared to ER+ or other mRNAs (Fig. 7d, f and Supplementary Fig. 8f). The infected cells also exhibit reduced NMD activity as indicated by significantly increased levels of PTC+ as compared to PTC- transcripts, which is consistent with the transient capsid expression experiments (Fig. 7e, g and Supplementary Fig. 8g). Overall, these findings suggest that in flavivirus infected cells, the viral capsid protein-PYM1 interaction may limit PYM1 function in regulating EJC occupancy to thereby alter gene expression in ER-associated compartments to aid viral replication.