Widespread PERK-dependent repression of ER targets in response to ER stress

The UPR (Unfolded Protein Response) is a well-orchestrated response to ER protein folding and processing overload, integrating both transcriptional and translational outputs. Its three arms in mammalian cells, the PERK translational response arm, together with the ATF6 and IRE1-XBP1-mediated transcriptional arms, have been thoroughly investigated. Using ribosome footprint profiling, we performed a deep characterization of gene expression programs involved in the early and late ER stress responses, within WT or PERK −/− Mouse Embryonic Fibroblasts (MEFs). We found that both repression and activation gene expression programs, affecting hundreds of genes, are significantly hampered in the absence of PERK. Specifically, PERK −/− cells do not show global translational inhibition, nor do they specifically activate early gene expression programs upon short exposure to ER stress. Furthermore, while PERK −/− cells do activate/repress late ER-stress response genes, the response is substantially weaker. Importantly, we highlight a widespread PERK-dependent repression program, consisting of ER targeted proteins, including transmembrane proteins, glycoproteins, and proteins with disulfide bonds. This phenomenon occurs in various different cell types, and has a major translational regulatory component. Moreover, we revealed a novel interplay between PERK and the XBP1-ATF6 arms of the UPR, whereby PERK attenuates the expression of a specific subset of XBP1-ATF6 targets, further illuminating the complexity of the integrated ER stress response.

PERK knockout (PERK −/−) cells have been useful for establishing PERK's function in cellular homeostasis maintenance under ER-stress 10 . Previous genome-wide studies have used mRNA expression profiling to define a transcriptional response following a 6 h ER-stress in PERK −/− and ATF4 −/− cells 11,12 . These experiments have shown PERK-dependent metabolic changes enabling the maintenance of redox potential under ER-stress 12 .
Continuing the wide body of research on the role of PERK in ER stress, we sought to understand the early and sustained PERK-dependent components of the UPR in a transcriptome-wide manner. While the translational arm of the UPR is fairly immediate, the impact of the transcriptional arms on cellular gene expression takes time to manifest. Thus, the different arms of the UPR generate a complex integrated regulation of gene expression programs in various stages of the response. Furthermore, while PERK is known to elicit an eIF2α phosphorylation-mediated global translational repression in response to ER stress, its role in controlling the translation of specific gene expression programs still remains elusive. We therefore chose to approach these questions in a manner that examines gene expression programs as an integration of transcription and translation.
In this study we examined the PERK-dependent dynamic alterations in gene expression programs following ER-stress using ribosome footprint profiling 13 on Wild-Type (WT) and PERK −/− Mouse Embryonic Fibroblasts (MEFs) 10 treated with ER stress. We used Thapsigargin (Tg), a SERCA inhibitor, for short (1 h or 2 h) or long (5 h or 8 h) treatments, to examine the gene expression programs that govern either early or late responses to ER stress.
We characterized three major gene expression programs in response to ER stress in PERK WT MEFs: Early induction, late induction and repression. We further show that all three programs are markedly compromised in the absence of PERK. We describe a widespread, PERK-dependent repression program consisting of ER target genes, including transmembrane proteins, glycoproteins, and proteins with disulfide bonds. Finally, we reveal that PERK attenuates a specific subset of XBP1-ATF6 target genes, thereby unraveling a complex interplay between the different arms of the UPR.

Results
exploring peRK dependent early and late eR stress gene expression programs. In order to identify gene expression programs that are activated or repressed in a PERK dependent manner during ER stress, we used Mouse Embryonic Fibroblast (MEF) cells that are either WT or were knocked out for PERK (PERK −/−) 10 . We treated these cells with Tg (1 µM), for various durations to assay both the early, immediate, response to ER stress (1 h or 2 h), as well as the late adaptive ER stress response (5 h or 8 h).
We first performed polysome profiling, to look at the global kinetics of the early and late translational responses ER stress. WT MEFs exhibited marked translational repression, which was maximal at 1 h of ER stress, and showed a time dependent relaxation, indicative of gradual adaptation (Fig. 1A). In contrast, PERK −/− MEFs showed no change in their global overall translation levels, and only at the 8 h timepoint a slight repression of translation was observed (Fig. 1B). However, even at 8 h, the degree of translational repression observed in PERK −/− MEFs did not resemble nearly any of the treatments of WT MEFs.
We then performed ribosome footprint profiling, to look at the gene expression programs that are regulated under the different conditions. Ribosome footprint profiling provides relative protein synthesis levels in each sample, allowing us to characterize mRNAs which translation is preferentially enhanced or repressed (see Materials and Methods).
Principle Component Analysis (PCA) of the data showed a clear separation between WT and PERK −/− MEFs, as well as a gradual temporal separation of ER stress durations (Fig. S1A). Classical ER stress target genes showed one of three expression patterns: an ER stress-mediated PERK independent upregulation (Fig. 1C,D), an ER stress-mediated upregulation that was partially PERK-dependent, such as in the cases of BiP (HSPA5) and XBP1 (Fig. 1E,F), or complete PERK dependence, with the examples of CHOP (DDIT3) and HERPUD1 (Fig. 1G,H). Finally, ATF4 translation showed the expected pattern of uORF (upstream Open Reading Frame) translation 14 , with a marked temporal increase in ribosome occupancy on its main ORF following ER stress in WT MEFs (Fig. S1B). PERK −/− MEFs did not show any ATF4 translation prior to, or upon ER stress (Fig. S1B), consistent with the lack of eIF2α-phosphorylation in these cells upon ER stress 10 .
Interestingly, examining the positional occupancy of ribosome footprints along ORFs showed a weak but significant accumulation of ribosomes at the 5′ ends of ORFs, until around position 140 bases downstream to the AUG (Fig. S2A). This is consistent with a similar observation from a recent study 9 . We and others have previously observed ribosome accumulation at 5′ ends of ORFs in similar positions following heat shock 15 and another proteotoxic stress condition 16 and showed that they are consistent with ribosomes paused in elongation. In the case of ER stress, however, while highly significant (p-value < 10 −300 in all timepoints relative to control, using a KS-test), the magnitude of pausing was much smaller than observed under heat shock. Specifically, a ~25-35% increase in pausing was observed under ER-stress at different times compared to control, vs. a 3-fold increase observed in heat shock 15 . Notably, PERK −/− MEFs did not show any 5′-end ribosome accumulation (Fig. S2B).
Ribosome footprint profiling revealed three major gene expression programs in response to ER stress. Next, we turned to ask what were the major gene expression programs that were characteristic of the early and late ER stress responses. To this end, we performed a clustering analysis of the Transcripts Per Million (TPM) values of the 1658 genes that changed their expression under ER stress in PERK WT cells (see Materials and Methods). The analysis revealed three major gene expression programs: (1) Early induction (Fig. 2, red cluster), containing 505 genes that were markedly upregulated in either 1 or 2 h following stress, and somewhat relaxed during 5 h and 8 h of stress; (2) Late induction (Fig. 2, blue cluster), including 495 genes that gradually increased in their expression and showed a maximal induction at the 5 h and 8 h timepoints; (3) Repression (Fig. 2, black cluster), containing 658 downregulated genes.
We compared our gene expression programs to a recently published polysome-sequencing data by Guan et al. of MEFs subject to either 1 h or 16 h of Tg treatment 8 , and reassuringly, our results were recapitulated in this www.nature.com/scientificreports www.nature.com/scientificreports/ dataset ( Fig. S3A-C). Furthermore, we compared the gene expression in our PERK −/− cells to MEFs treated with a PERK inhibitor following 16 h of ER stress from Guan et al., and found that genes that were induced or repressed in PERK −/− cells at 8 h of Tg were also significantly induced or repressed respectively at 16 h Tg-stressed cells treated with a PERK inhibitor (Fig. S3D).
To further substantiate the robustness of the major gene expression programs that we have identified in WT MEFs, we performed an additional experiment in NIH3T3 cells. NIH3T3 cells were exposed to ER stress (Tg) for either 2 h or 7 h, followed by polysome profiling and ribosome footprint profiling. Polysome profiling showed a similar trend of overall translational inhibition at 2 h of ER stress, which relaxed at the 7 h timepoint (Fig. 3A). Using ribosome footprint profiling, we found that the three major gene expression programs were largely recapitulated in NIH3T3 cells (Figs 3B-D and S4), demonstrating the generality of our findings.
Both early and late induction gene expression programs are peRK dependent. We next characterized which pathways were enriched within the induction gene expression programs. Functional analysis of the early induction program showed an enrichment for specific classes of histones as well as other DNA-binding proteins (Table S1A), as has been reported previously 9 . The late induction program was significantly enriched with genes annotated as "response to ER stress", as expected (Table S1B). Additionally, the late induction program showed an enrichment for amino acid biosynthetic processes, in agreement with previous reports by Harding et al. 12 using a different ER stressing agent, Tunicamycin (TM).
Examination of the early induction cluster showed that the increase in expression at the 1 h and 2 h timepoints was completely PERK dependent, while at the 5 h and 8 h timepoints, a very weak yet significant induction was observed in PERK −/− cells (Figs 4A,B and S5A,B). Whereas the early induction program genes showed a median induction of 80% and 89% at the 5 h and 8 h of ER stress respectively, their median fold change was 12% and 11% in the PERK −/− cells at these timepoints.
Late induction cluster genes showed a slight aberrant repression in the early timepoints in PERK −/− MEFs (Fig. S5C,D), while at the later timepoints, a small yet significant induction was observed (Fig. 4C,D). Whereas the median fold change of the late induction genes was 128% and 136% in the PERK WT cells following 5 h and www.nature.com/scientificreports www.nature.com/scientificreports/ 8 h of Tg respectively, only 9% and 26% median induction values were observed in PERK −/− cells at these timepoints.
Taken together, we observed that PERK −/− cells were highly impaired in eliciting the desired gene expression program in the early response, (Figs 4 and S5). Interestingly, at the late timepoints, cells were still unable to produce the full extent of activation of the required gene expression programs and showed a poor level of induction ( Fig. 4A-D).

PERK-mediated repression of ER targets during ER stress is a major component of the UPR.
We next turned to characterize which pathways were enriched within the repression gene expression program. The repression program showed a significant enrichment for many functional annotation categories, which mainly converged upon two main themes (Fig. 5A, Table S1C): membrane, transmembrane, signal peptide-containing proteins and cell surface proteins, which represent a major class of ER-targeted proteins; and disulfide-bond containing proteins and glycoproteins, representing proteins undergoing post-translational modifications inside the ER.
Next, we examined whether the repression of gene expression was dependent on PERK. In sharp contrast to the PERK WT cells (Fig. 5B), no downregulation was observed for the repression program genes in the PERK −/− cells at the early timepoints ( Fig. 5C). At the later stage, a significant, yet minimal, repression was observed www.nature.com/scientificreports www.nature.com/scientificreports/ for these genes; with a median repression of 18% and 12% in the 5 h and 8 h timepoints respectively in PERK −/− cells, compared to a median of 50% and 51% repression in PERK WT cells (Fig. 4E,F).
These results indicate that, in addition to the global repression of protein synthesis caused by the PERK-induced eIF2α phosphorylation, ER targets are significantly and selectively further repressed during ER stress, in a PERK dependent manner.
Translational regulation plays a major role in the repression of ER targets during ER stress. As ribosome footprint profiling provides an integrated output of mRNA-level and translation-level changes, we further asked whether the PERK-dependent repression of ER targets that we identified occurs at the mRNA level or at the level of translation. Guan et al. 8 have reported that ER-translated mRNAs were repressed in ER stress, with different subsets regulated at the mRNA level and at the level of translation. To further characterize the relationship between mRNA-level and translation-level regulation for the case of ER target repression, we separated  www.nature.com/scientificreports www.nature.com/scientificreports/ genes that were repressed at 1 h in the Guan et al. dataset into two groups: genes that showed repression both at the mRNA and the translational levels, and genes whose repression was mainly at the level of translation (Fig. S6A, see Materials and Methods). We note that, overall, the Guan et al. dataset showed a strong translational response at 1 h, while at 16 h of ER stress, mRNA levels and translation levels were highly correlated (Pearson R 2 of 0.54, compared to 0.22 at 1 h, Fig. S6A-C). Importantly, functional enrichment analysis of both groups showed . CDF plots of the difference between the LFC at the translation level (Ribo-seq) and the LFC at the mRNA level (RNA-seq) are markedly shifted showing that mRNAs are largely repressed at the level of translation. Repressed mRNAs were defined as mRNAs with 1.5 fold or more translational repression. Both groups were highly enriched for ER targets (Table S2). This analysis was robust to various choices of repression parameters (Fig. S6J-M).
www.nature.com/scientificreports www.nature.com/scientificreports/ enrichment for ER targets (Table S2A,B). Further examination of these groups at 16 h of ER stress showed that their mode of regulation was largely maintained: while the first group was still repressed both at the mRNA and the translational levels, the second group remained mostly translationally repressed (Fig. S6B-E).
Next, we analyzed two additional datasets of mRNA-seq and ribosome footprint profiling (Ribo-seq), generated in NIH3T3 (Kwak lab, GSE103667) and HEK293T cells 17 treated with a short ER stress, 1.5 h and 2 h respectively, using Tg. We first similarly defined groups of genes that were either repressed both at the mRNA level and translationally, or only translationally (see Materials and Methods). Importantly, in both datasets, a much smaller fraction of genes showed mRNA-level repression (15.3% and 6.6% of the repressed genes at a 1.5 fold cutoff, respectively), while the vast majority were translationally repressed (Figs 5D,E and S6F-I, see Materials and Methods). Additionally, in both datasets, both translationally repressed genes and genes inhibited at the mRNAand translation levels, were highly enriched for ER targets (Table S2C-F). We note that this analysis was robust to the choice of the repression cutoff (see Materials and Methods, Fig. S6J-M, Table S2G-N).
These results further support the notion that the repression of ER targets upon ER stress is a widespread phenomenon, occurring at the level of translation for a substantial part of the genes, with mRNA-level contribution for a subset of the genes.
peRK-mediated repression of eR targets during eR stress occurs in other cell types. As we have shown above, the repression of ER targets occurs at the level of translation, as well as at the mRNA level for a subset of genes. This fact allowed us to examine whether we can observe any downregulation of our repression gene expression program, which was highly enriched with ER targets, while examining other mRNA expression profiling datasets.
We first wanted to verify that the repression gene expression program we identified above was general to ER stress, and is not Thapsigargin specific. Analysis of two additional mRNA expression profiling datasets of MEFs treated with Tunicamycin (TM), for short, intermediate 18 or long 11 treatments, showed a mild yet highly significant downregulation of our identified repression program ( Fig. S7A-C). Furthermore, analysis of mRNA expression profiling from DTT treated MEFs (Kaufman lab, GSE84450) also recapitulated this repression (Fig. S7D).
Therefore, our identified gene repression program was not Tg specific, but rather a general characteristic of the response to ER stress.
We then analyzed several additional mRNA expression profiling datasets to see if the repression program can be generalized to other specialized cell types, beyond mouse fibroblasts and human HEK293T cells (Figs 5E and S6I,K,M). Indeed, we found that significant inhibition of the repression program was recapitulated in intestinal stem cells subject to 4 h of TM 19 (Fig. S7E), and in TM treated glioma-derived stem cells, and glioblastoma cells (Dorsey lab, GSE102505) (Fig. S7F,G), but not in livers treated with 6 h of TM 20 (Fig. S7H).
These analyses indicate that overall, the repression gene expression program is general, and occurs across many different, although not all, cell types.
PERK-mediated repression of ER targets during ER stress largely involves eIF2α phosphorylation, and is partly dependent on ATF4. We then turned to explore the potential role of regulators downstream to PERK in the mediation of ER target repression. To that end, we first analyzed a dataset of mRNA expression from a long ER stress treatment (TM, 12 h), performed in either WT, PERK −/− or eIF2α-S51A MEFs (Oyadomari lab, GSE49598). We observed a mild but significant repression at the level of mRNA in WT MEFs, both for our defined repression program (Fig. S8A), as well as for ER targets (signal peptide containing proteins, Fig. S8D). This repression was completely abrogated in PERK −/− MEFs (Fig. S8B,E), consistent with our ribosome footprint profiling data above. Furthermore, eIF2α-S51A MEFs also showed no repression of ER targets, as well as of the repression program (Fig. S8C,F).
Interestingly, our analysis of a dataset of Arginine deprivation in HCT116 and HEK293T cells 21 , a condition that is known to lead to phosphorylation of eIF2α through GCN2, also found a mild yet significant downregulation of the repression program ( Fig. S8G-I). However, this repression did not depend on GCN2 (Fig. S8J), since it remained slight but significant in GCN2 knockout HEK293T cells.
These results indicate that eIF2α phosphorylation downstream of PERK plays a major role in the preferential repression of ER targets (see discussion below).
The ATF4 transcription factor is an important effector in the ER stress response downstream to PERK, in an eIF2α phosphorylation dependent manner (Fig. S1B). Since we observed the repression of ER targets as early as 1 h of stress, and a large part of it was translational, it is less likely to be dependent on ATF4. However, since mRNA-level repression of some ER targets does occur (Fig. S6), we sought to further examine whether ATF4 plays a role in this repression. Analysis of a dataset of 8 h ER stress (TM treatment) in either WT or ATF4−/− MEFs 11 showed a significant downregulation of our identified repression program in both WT and ATF4−/− cells (Fig. S8K,L). These results indicate that the PERK-mediated repression program is largely independent on ATF4. Nevertheless, when we examined the set of ER targets (signal peptide encoding mRNAs), we observed a partial relief in the repression of a subset of them in ATF4−/− MEFs (Fig. S8M,N). Thus, the repression of a subset of ER targets is partially dependent on ATF4.
Additional pathways subject to PERK-mediated repression. To ask whether the repression program contains additional pathways, we removed all 265 ER targets from this cluster, leaving a gene set of 393 repressed genes, and repeated the functional enrichment analysis. We found that this repressed set is enriched for cyclins (Table S1D). This enrichment in repressed cyclins fits well with the cell cycle arrest that occurs upon ER-stress 22 . Further analysis showed that cyclins were repressed as a group (Fig. S9A), and this repression did not occur in PERK −/− cells. Additionally, we observed an enrichment for LSM domain proteins (Table S1D), a family of RNA-binding proteins. For this group as well, repression was PERK dependent (Fig. S9B).
www.nature.com/scientificreports www.nature.com/scientificreports/ PERK attenuates a subset of XBP1 and ATF6 targets, revealing a complex interplay between the UpR arms. At the late stages, transcriptional responses mediated by the ATF6 and IRE1-XBP1 arms of the UPR are expected to take effect. We found that the expression levels of both XBP1 and ATF6 were induced by ER stress, and this induction was partly diminished (40% less induction) in PERK −/− cells (Figs 1E and S10), in agreement with previous reports 20 . We wanted to better understand the potential interplay between the XBP1-and ATF6-mediated response pathways and PERK. We therefore examined the expression pattern of bona-fide XBP1 and ATF6 target gene sets, which were originally defined by Shoulders et al. 23 without eliciting ER stress, within our data. In the early timepoints, the expression of XBP1 and ATF6 targets was largely unchanged (Fig. S11A-C). This is consistent with the delayed accumulation in their transcriptional output. In the late timepoints, however, we found that while 74% of XBP1-ATF6 targets were induced in PERK −/− cells following ER stress, only 56% were induced in PERK WT cells (at 8 h, Fig. 6A,B). Similar trends were observed when we considered XBP1 targets and ATF6 targets separately (Fig. S11D,E). Thus, it seems that PERK might have an inhibitory effect on a subset of XBP1 and ATF6 target genes.
We further examined the expression of these two gene sets in the Guan et al. polysome-seq data 8 , where MEFs have been ER stressed for 16 h, with or without the presence of a PERK inhibitor for the last 4 h of the stress. Our analysis revealed that, consistent with our data above, PERK-independent XBP1-ATF6 targets were induced at 16 h of ER stress irrespective of PERK inhibition (Fig. 6C,D, red curve). The PERK-attenuated XBP1-ATF6 targets, however, showed no induction at 16 h post stress (Fig. 6C, blue curve). Nevertheless, the expression pattern of these genes in the PERK inhibitor treated cells demonstrated that PERK inhibition alleviated this attenuation (Fig. 6D, blue curve), further substantiating that this subset of XBP1-ATF6 targets is subject to PERK-mediated attenuation.
Similar trends were observed when we examined XBP1 targets separately (Fig. S11F-H).
To understand if this PERK-mediated attenuation occurs at the level of transcription or translation, we analyzed mRNA expression data from TM treated WT or PERK −/− MEFs (12 h, Oyadomari lab, GSE49598). www.nature.com/scientificreports www.nature.com/scientificreports/ Here too, we observed the same relief in PERK-mediated inhibition of the subset of PERK-attenuated genes we identified above (Fig. S12A,B), indicating that transcriptional regulation is involved in this attenuating effect. Moreover, our analyses found that similar trends of induction vs. attenuation for the PERK-independent and PERK-attenuated XBP1-ATF6 targets defined above, were evident and highly significant in many mRNA expression datasets, including response to TM 18 (Fig. S12D) and DTT (Kaufman lab, GSE84450, Fig. S12E), as well as in other cell types, including Tg treated intestinal stem cells 19 (Fig. S12F), and human derived glial cells (24 h of TM treatment, Dorsey lab, GSE102505, Fig. S12G).
Nevertheless, examination of mRNA expression data from Guan et al. MEFs 8 , where PERK inhibitor was added only during the last 4 h of a 16 h Tg treatment, showed a more complex picture. While these cells showed alleviation of attenuation at the translation level (Fig. 6D), no relief was observed at the mRNA level (Fig. S12H,I).
Taken together, our analyses suggest that the PERK-mediated attenuation of a subset of XBP1-ATF6 targets contains both translational and transcriptional components, as initial relief in attenuation is observed at 4 h of PERK inhibitor only at the level of translation. At longer durations of ER stress, transcriptional regulation downstream of PERK seems to play a major role in the attenuation of this subset of XBP1-ATF6 targets.
Analysis of the potential contribution of PERK pathway effectors showed that eIF2α-S51 A MEFs (Oyadomari lab, GSE49598) largely recapitulate the relief of the PERK-mediated attenuation of this subset of XBP1-ATF6 targets (Fig. S12C) observed in PERK −/− MEFs. Furthermore, in ATF4 −/− MEFs 11 this attenuation also showed significant alleviation (Fig. S12J,K). These results indicate that eIF2α phosphorylation is largely involved in mediating the PERK-dependent attenuation of this subset of XBP1-ATF6 targets, with a significant contribution of the ATF4 branch to the phenomenon, in line with the combined transcriptional-translational regulatory effect downstream of PERK we found above.
The fact that this phenomenon had similar characteristics to the PERK-mediated repression of ER targets, prompted us to speculate that these two phenomena might be related. Indeed, functional enrichment analysis showed that many of these PERK-attenuated XBP1-ATF6 targets are ER targets (Table S3B,C). Therefore, it is possible that this subset of targets is subject to the PERK-dependent ER targets repression program described above. The PERK-attenuated XBP1-ATF6 target set was selectively enriched with sterol and lipid biosynthesis genes (Table S3B), which are ER residents, and their repression following ER stress has been previously reported 24 . ER targets were also enriched within the PERK-independent set (Table S3A,C). There, however, selectively enriched ER targets included ER stress response pathway genes (Table S3A), indicating that classical ER stress response regulators are not subject to attenuating regulation via PERK.
Together, our analyses revealed a major gene expression program including PERK-mediated repression of ER targets. This program also attenuates a subset of XBP1-ATF6 targets, unravelling another level of the complex interplay between the different arms of the ER stress response.

Discussion
The ribosome profiling high-throughput data presented in our study provides a genome-wide temporal view of ER stress gene expression programs, integrating transcriptional and translational outputs, in both WT and PERK −/− MEFs. Our data and analyses revealed composite dynamics of the integrated expression response to ER stress, orchestrated by PERK, and showed a complex combination of gene repression and induction under ER stress.
Our data suggests that the impairment in specific gene activation, as well as repression, in the absence of PERK are very widespread, and unravel a complex interplay between the three arms of the ER stress response.
Wang et al. 25 has reported a strong and immediate global translational repression within ten minutes after pharmacological ER stress induction, using nascent chain translation live imaging. The polysome profiling above (Fig. 1A,B) demonstrated that this repression was markedly noticeable 1 h post Tg treatment and was gradually relieved.
Our findings capture a strong, selectively enhanced, repressive gene expression program following ER stress and demonstrate its impact on ER targeted proteins (Fig. 5), which seems to be an additional hallmark of the ER stress response. More specifically, we found that proteins undergoing ER-dependent post-translational modifications such as glycosylation and disulfide-bonds, as well as transmembrane and membrane proteins are enriched among genes repressed under ER-stress. Importantly, we have shown that this repression occurs largely at the level of translation, with a small contribution of mRNA-level repression (Figs 5D,E and S6).
Several microRNAs have been shown to be induced in response to ER stress 26,27 . However, examination of the putative targets of these miRNAs did not explain the repression of ER targets observed (data not shown). Rapid changes in translation-dependent compartmentalization of ER target protein synthesis has been previously demonstrated, with an increase in the proportional translation of ER targets by cytosolic ribosomes 9 . This change in translation compartmentalization was shown to be transient, with a robust change at 30 minutes post Tg treatment, and an almost complete reversal by 1 h following ER stress induction 9 . Significant repression in ER stress has been later reported for these genes 8 . Here we show that there is major translational repression of ER targets, which is highly PERK dependent. Importantly, even though it has an mRNA-level component, translation is the major regulatory level for this phenomenon, and thus it was strongly evident from our ribosome footprint profiling data analysis. It is possible that the re-localization of ER targets 9 is an initiating step in the cellular response to ER stress, with a second, yet to be characterized, PERK-mediated mechanism of repression. Moreover, our analyses point to eIF2α phosphorylation as critical for this repression. One possible model that would be consistent with the data is that PERK-mediated phosphorylation of eIF2α factors primarily affects the translational machinery pool at the vicinity of the ER, thus affecting mainly ER targets. Our analysis of ribosome footprint profiling data from Arginine deprived cells 21 has shown that while this stress elicits mild ER target repression, this repression was GCN2-independent. It is therefore possible that in the case of GCN2, eIF2α phosphorylation, (2019) 9:4330 | https://doi.org/10.1038/s41598-019-38705-5 www.nature.com/scientificreports www.nature.com/scientificreports/ which occurs throughout the cytoplasm, does not lead to a GCN2-dependent repression of ER targets. However, whether indeed this model underlies PERK-mediated repression of ER targets remains to be explored.
Past research has indicated Cyclin D1 to be translationally repressed during ER stress in a PERK dependent manner 22 . Interestingly, we observed a wider PERK-mediated repression of cyclins expression, including cyclins D,E,F,G and I, thereby expanding existing knowledge regarding PERK mediated cell-cycle arrest. The cell-cycle arrest has immense implications, which may be essential for cell survival under chronic stress 22 . Furthermore, we found a PERK-dependent LSM domain protein repression, pointing towards a possible regulation at the level of RNA processing. LSM-proteins are known to aggregate within processing-bodies under various stress conditions altering mRNA homeostasis 28 . In addition, we see an early transient PERK dependent induction of DNA-binding proteins, mainly zinc-finger and histones, in agreement with previous reports 9 . These yet to be explored changes may be related to ER stress-mediated chromatin changes. For example, ATF4 target gene activation has been linked to the histone lysine demethylase KDM4C epigenetic rewiring in cancer 29 .
Previous studies have underlined PERK reinforcement of both the IRE1-XBP1 and ATF6 UPR arms 20,30 . It has been shown previously that XBP1-spliced accumulation is PERK dependent 30 . Furthermore, the transcription of selected ATF6 targets was augmented 6 h after ER stress initiation in the liver, and the processing of ATF6 has been shown to be dependent on ATF4 20 . Recent research has reported PERK-mediated IRE1 inhibition via the RPAP2 phosphatase, which has a role in diverting adaptive ER stress response towards apoptosis after 16 h of stress 31 . Consistently, here we observed that both XBP1 and ATF6 expression are induced in a temporal manner following ER stress, and their induction is partially reduced in the absence of PERK (Figs 1E and S10). Thus, several reinforcing positive interactions between PERK and other UPR arms have been demonstrated. Nevertheless, we report here that PERK mediates an attenuating effect on a subset of XBP1 and ATF6 target genes (Figs 6 and S11). As many of these genes are ER targets (Table S3), they could be subject to the same mechanism of PERK-mediated repression of ER targets we highlight above. Notably, sterol biosynthetic genes, which are directly upregulated by XBP1-ATF6 activation, were specifically attenuated in a PERK-dependent manner (Table S3B). Indeed, previous findings have linked the UPR to down-regulation of the sterol biogenesis pathway, and subsequently cholesterol levels, via suppression of SREBP 24 . This PERK-mediated attenuation that we found was not gradually relieved at 8 h (Fig. 6A), nor by 16 h of ER stress, as seen from our re-analysis of the dataset by Guan et al. 8 (Fig. 6C). The repression gene expression program shows a similar trend (Figs 5B and S3C). Furthermore, both phenomena show similar characteristics of eIF2α phosphorylation dependence and a small yet significant contribution of ATF4. This raises the possibility that these two phenomena might indeed be related, such that ER target repression explains the PERK-mediated attenuation of the subset of XBP1-ATF6 targets. While repression of ER targets is primarily translational and starts early, transcriptional effects kick in later on, and thus attenuation of this subset of the transcriptional UPR arms is evident at the late time points. While future studies will unravel the purpose of this attenuating affect, and its role in the adaptation to ER stress, it is clear that the complex cross-talk between the three arms of the UPR, the interplay between transcriptional and translational regulation, and the function of PERK in this interplay, play a critical role in shaping the cellular gene expression program in response to ER stress.

Materials and Methods
Cell culture and stress. Mouse fibroblast NIH3T3 cells were grown in DMEM with 10% FBS and pen-strep.
MEFs were grown in DMEM with 10% FBS supplemented with 0.1 mM Non-Essential Amino Acids (NEAA) and 0.05 mM beta-Mercaptoethanol. NIH3T3 cells or MEFs were plated at low density (2 million cells in 15 cm plate one day prior to harvesting), and treated with Thapsigargin (at 1μM or 200 nM for MEFs and NIH3T3 cells respectively) for short (1 h or 2 h) or long (5 h and 8 h for MEFs, 7 h for NIH3T3) durations, to induce early or late ER stress responses.
Ribosome footprint profiling and Polysome profiling. Ribosome footprint profiling was performed as in Shalgi et al. 15 , with the following modification: cells were grown and treated with Tg for the indicated times, then trypsinized, pelleted by 5 minutes at 1000RPM centrifugation at 4 °C, and flash frozen in LN2, without prior Cycloheximide treatment. Cycloheximide was also omitted from the buffers. Polysome profiling was performed as in Shalgi et al. 15 .
Data analysis. Ribosome footprint profiling provided quantification of number of ribosomes per transcript in each timepoint, resulting in the quantification of protein synthesis levels. We note that the protein synthesis levels are relative, rather than absolute, since the overall number of translating ribosomes is decreased in response to ER stress in WT MEFs, due to eIF2α phosphorylation 4 , as also observed in Fig. 1. Nevertheless, we were interested in characterizing, given the limited pool of overall translating ribosomes, which genes and which pathways are relatively more enhanced or further repressed under each condition.
Ribosome footprint profiling data was mapped to mm10 version of the genome using RefSeq CDSs. Transcripts shorter that 100 nucleotides were filtered out and 30 nucleotides were clipped form the start and end of each CDS, similarly to Ingolia et al. 32 . Ribosome footprint sequences were trimmed to maximal size of 34 and polyA sequences from the polyA tailing step were removes, such that footprints of lengths 22-34 were considered. FASTQ files were then filtered for rRNAs and tRNAs using STAR 33 . Expression levels were quantified using RSEM 34 after mapping to clipped CDSs using Bowtie2 35 to produce transcript and gene level TPM (Transcript Per Million) values (see Table S5).
Lowly expressed genes with an average TPM below 2 across all samples were filtered out and expression of all other genes was thresholded to a TPM of 4. For each experimental group (PERK WT or −/− MEFs, or NIH3T3 cells), fold changes were calculated as a ratio between each condition and its respective control. Regulated genes were designated as genes which changed at least two-fold relative to their respective controls.
www.nature.com/scientificreports www.nature.com/scientificreports/ PCA was generated using R DESEQ2 36 , using the function plotPCA. Hierarchical clustering of genes was done based on spearman correlation using the clustergram MATLAB function.
Additional RNA-seq ER stress datasets (detailed in Table S4) were downloaded from the GEO database. RNA-Seq reads were filtered for rRNAs using STAR 33 . The remaining reads were mapped to hg19 or mm10 versions of the human and mouse genomes respectively, using STAR 33 . Expression levels were quantified using RSEM 34 . TPM values were averaged between sample replicates.
Ribosome footprint profiling data were processed as described above (mapped to hg19 or mm10 versions of the human and mouse genomes, respectively).
Microarray ER stress datasets (detailed in Table S4) were downloaded from GEO. Fold changes for the dataset GSE54581 were derived using the Affymetrix transcriptome analysis console, with CEL files as input. Fold changes for the datasets GSE49598 and GSE84989 were generated using normalized gene expression signals as provided by GEO. Replicates from each condition of interest were averaged, and gene expression fold changes were calculated.
To characterize the relationship between mRNA-level and translation-level regulation for the case of ER target repression, we separated the repressed genes (1.5 fold) in three datasets (Guan et al. 8 , Woo et al. 17 , and Kwak lab NIH3T3 cells) into two groups: genes that showed repression both at the mRNA and the translational level (mRNA-level repression and translation-level repression 1.5 fold or more), and genes whose repression was mainly at the level of translation (mRNA-level repression less than 1.5 fold and translation-level repression more than 1.5 fold). We examined the number of genes in each group as well as the difference between the log2 fold change in the translation-level and the log2 fold change in the mRNA-level. We repeated the above analyses using a range of threshold values from 1.3-2.1 and verified that the results were indeed robust to the choice of threshold (Fig. S6J-M). Additionally, we repeated the functional enrichment analysis using sets with different thresholds for these datasets, and all translationally repressed gene groups, as well as mRNA and translationally repressed gene groups, were found to be significantly enriched for ER targets (Table S2).
CDF plots were used to compare the cumulative distributions of log2 fold changes (LFCs) in the expression of gene subsets compared to their respective control, to the overall distribution of LFCs of the entire transcriptome (expressed genes, defined as above). Significance for different gene groups was calculated using the student t-test statistical test.
Pathway enrichment analyses were conducted using DAVID 6.8 functional annotation tool 37 , using all expressed genes as background. Pathways with Benjamini corrected p-value < 0.05 were designated as enriched.
Signal peptide encoding mRNAs were defined using the signalP program 38 . XBP1-ATF6 targets were defined by Shoulders et al. 23 following specific activation of XBP1, ATF6, or both of them together (termed XBP1-ATF6 targets). We further filtered these sets for genes that were not expressed in our cells (TPM < 2). PERK-independent XBP1-ATF6 targets were defined as XBP1-ATF6 targets that were upregulated in both PERK WT and PERK −/− following 8 hours of Tg, and PERK-attenuated XBP1-ATF6 targets were defined as XBP1-ATF6 targets up-regulated in PERK −/− and down-regulated in PERK WT. These sets were further subject to functional enrichment analysis using DAVID, and the resulting enriched terms (Benjamini corrected p < 0.05) were examined for the uniqueness in only one of the sets (Tables S3A,B) or intersected in both sets (Table S3C).

Data Availability
Original and processed data were deposited in GEO under GSE118660. Processed data is also available in Table S5. Additional datasets used are listed in Table S4.