Deconvoluting Stress-Responsive Proteostasis Signaling Pathways for Pharmacologic Activation using Targeted RNA-sequencing

Cellular proteostasis is maintained by stress-responsive signaling pathways such as the heat shock response (HSR), the oxidative stress response (OSR), and the unfolded protein response (UPR). Activation of these pathways results in the transcriptional upregulation of select subsets of stress-responsive genes that restore proteostasis and adapt cellular physiology to promote recovery following various types of acute insult. The capacity for these pathways to regulate cellular proteostasis makes them attractive therapeutic targets to correct proteostasis defects associated with diverse diseases. High-throughput screening (HTS) using cell-based reporter assays is highly effective for identifying putative activators of stress-responsive signaling pathways. However, the development of these compounds is hampered by the lack of medium-throughput assays to define compound potency and selectivity for a given pathway. Here, we describe a targeted RNA sequencing (RNAseq) assay that allows cost effective, medium-throughput screening of stress-responsive signaling pathway activation. We demonstrate that this assay allows deconvolution of stress-responsive signaling activated by chemical genetic or pharmacologic agents. Furthermore, we use this assay to define the selectivity of putative OSR and HSR activating compounds previously identified by HTS. Our results demonstrate the potential for integrating this adaptable targeted RNAseq assay into screening programs focused on developing pharmacologic activators of stress-responsive signaling pathways.

The development of pharmacologic activators of stress-responsive signaling pathways has primarily been pursued using high-throughput screening (HTS) approaches that employ cellular transcriptional reporters of target genes activated downstream of specific stress pathways including the ATF6 signaling arm of the UPR and the HSF1-dependent HSR 20,[25][26][27][28][29][30] . While this approach has effectively identified many putative activators of these pathways, the further development and characterization of these HTS hits is often hampered by complications including reporter interference, lack of compound selectivity for a given pathway, or reporter constructs not reliably reporting on activation of the entire protective transcriptional program 24,28,31,32 . Without proper tools to assess selectivity across broad stress signaling pathways, it is difficult to determine whether previous HTS have identified effective compounds that selectively activate these pathways.
One strategy to increase the efficiency of identifying specific pathway activators from many screening hits is to incorporate upstream transcriptional profiling to first define the activation spectrum among stress responsive signaling pathways. The benefits of this approach have been demonstrated with the recent establishment of compounds that preferentially activate the ATF6 signaling arm of the UPR, where multiplex gene expression (MGE) profiling was integrated into a screening pipeline centered on cell-based transcriptional reporters 25 . However, despite the evidence highlighting the benefit of incorporating transcriptional profiling into screening platforms, cost effective strategies to profile stress-responsive signaling pathway activation in a medium-throughput format are currently lacking.
Defining the magnitude and repertoire of activation among stress-responsive signaling pathways for a given stimulus is complicated by multiple challenges. Stress-responsive genes can be regulated by multiple signaling pathways, making it difficult to discern pathway activation by tracking the expression of a single gene.
For example, the OSR target gene HMOX1 can be regulated by multiple stress-responsive transcription factors including NRF2 (OSR), HSF1 (HSR), and NF-kB 33 . Furthermore, many stress-responsive signaling pathways have overlapping sets of target genes, challenging the ability to define selective activation of a certain pathway.
For example, the majority of genes regulated by the UPR-associated transcription factor ATF6 are also activated, albeit to lower extents, by the alternative UPR-associated transcription factor XBP1s, thus making it difficult to deconvolute specific activation of these pathways by monitoring expression of a single gene 34 .
Additionally, different stress-responsive signaling pathways induce target genes to varying extents. For example, HSF1 (HSR) target genes can be induced >10-fold higher than UPR target genes 17,34 . These properties of stress-responsive signaling challenge the ability to monitor activation of specific pathways using strategies such as geneset enrichment analysis (GSEA), which is biased towards pathways that elicit larger transcriptional responses and does not easily deconvolute overlapping stress-responsive transcriptional programs. Furthermore, GSEA requires whole transcriptome RNA sequencing (RNAseq) profiling to define pathway activation, limiting its application as a medium-throughput screening approach. One potential strategy to address the above challenges is to monitor activation of specified sets of stress-responsive genes regulated downstream of different stress-responsive signaling pathways, wherein pathway activation is defined by the grouped behavior for all relevant target genes. This strategy requires measuring multiple genes activated downstream of different stress-responsive signaling pathways in a cost-effective assay.
Recent advances in RNA sequencing have demonstrated the potential for this approach to be integrated into drug discovery pipelines. For example, the DRUG-seq platform established a cost-effective strategy to profile compounds in a high-throughput format, providing a powerful approach to define compound mechanism of activation and selectivity 35 . However, this approach requires specialized equipment that would make it difficult to implement in most research laboratories. In contrast, the targeted RNAseq platform, described in the last 5 years (previously described as Capture-seq 36 ), provides a unique opportunity to screen expression of 100-1000 genes in a cost-effective, medium-throughput format. Since this approach uses targetspecific generation of sequencing libraries, targeted RNAseq demonstrates improved sensitivity for low-copy transcripts, potentially providing a larger dynamic range to track changes in both low and highly-expressed genes. Furthermore, targeted RNAseq avoids background issues caused by non-specific probe binding or probe cross-hybridization found in such technologies as microarrays 37 . As such, this approach has been used in diverse contexts including measuring expression of alleles in plant populations 38 , detection of gene fusions in solid tumors 39 , and monitoring activation of cell death pathways 40 .
Here, we describe a targeted RNAseq assay designed to define activation of stress-responsive proteostasis pathways in a medium-throughput format. We show that this approach allows accurate deconvolution of stress-responsive pathway activation induced by diverse chemical genetic and pharmacologic agents. Furthermore, we demonstrate the potential for this approach to define the selectivity of pharmacologic activators of stress-responsive signaling pathways by profiling the selectivity of compounds identified by highthroughput reporter screening to activate the OSR-associated transcription factor NRF2 or the HSR-associated transcription factor HSF1 20,26 . Ultimately, our results show that targeted RNAseq profiling is a highly adaptable strategy that can be efficiently incorporated into HTS pipelines and downstream compound development to improve the establishment of pharmacologic activators of specific stress-responsive signaling pathways.

A Targeted RNAseq Assay to Monitor Activation of Stress-Responsive Proteostasis Pathways
To establish a targeted RNAseq assay for monitoring activation of stress-responsive proteostasis pathways, we first defined genesets predicted to accurately report on the activation of the predominant proteostasis pathways: the HSR, OSR, the three signaling arms of the UPR, and the ISR (Fig. 1A). We used published transcriptional profiles from human cells subjected to pharmacologic, chemical genetic, and genetic manipulation of specific stress-responsive signaling pathways to identify sets of proteostasis genes induced by each pathway. From this data, we selected 10-20 reporter genes activated downstream of the HSR 17 , OSR 41 , and the three arms of the UPR regulated by the ER stress sensing proteins IRE1, ATF6, and PERK 34, 42-44 ( Fig.   1A and Table 1). Importantly, the geneset that reports on activation of the PERK signaling arm of the UPR can also be used to monitor activation of the ISR, as both are activated through a process involving phosphorylation of the α subunit of eukaryotic initiation factor 2 (eIF2α) and the activity of the ATF4 transcription factor 18,44 . In addition, we include stress-responsive genes activated downstream of other stressresponsive signaling pathways including the Hypoxic Stress Response 45 , NFκB signaling 41,46 , and the poorly defined mitochondrial unfolded protein response (UPR mt ) [47][48][49] , providing additional readouts of stress-pathway activation. Our final gene panel consists of 150 target genes shown in Table 1.
We used the established targeted RNAseq profiling approach with this custom gene panel to define the activation of stress-responsive signaling pathways in HEK293T cells grown in 96-well plates subjected to multiple conditions predicted to activate the different stress-responsive proteostasis pathways shown in Fig. 1A (see Table 2 and Table S1). Briefly, we isolated RNA from these cells and generated cDNA libraries using a standard reverse transcriptase reaction. We then amplified our genes of interest for sequencing using targeted primer sets directed to the 150 genes in the panel (Fig. 1B). Amplicons were then isolated and pooled for sequencing using the MiSeq desktop sequencer at a target depth of 50 million paired-end reads for all pooled samples. Overall alignment of reads reflected the specific nature of this approach with over 93% of reads aligning to target regions, which is significantly greater than that observed in conventional whole transcriptome RNAseq experiments (Fig. S1A). Our desktop MiSeq sequencing run yielded a median of ~580,000 reads per sample (19 conditions in triplicate, 57 samples total), which is approximately 1% of the number of reads aligned per sample with whole transcriptome RNAseq on the same platform (approximated to 44-50x10 6 reads per sample). Additionally, per gene target, our targeted RNaseq assay yielded an median of 721.7 aligned reads (median total >41,000 reads) across all treatment conditions included in these analyses ( Fig. S1B-C). All of the data from this targeted RNAseq assay is included in Table S2.
Across replicate samples, we observed high reproducibility with all treatments demonstrating an R 2value above 0.7 and most having an R 2 >0.85 ( Fig. S1D-E). From aligned count data, we performed unbiased clustering across all treatment conditions to determine our ability to accurately define different stress-signaling pathways (Fig. 1C,D). This analysis shows that genes regulated by the HSR, OSR, and the three arms of the UPR (IRE1/XBP1s, ATF6 and PERK/ISR) generally cluster together, reflecting their similar regulation across the various stress conditions. However, despite this clustering, we observe significant overlap between genesets, reflecting their integrated activation in response to diverse types of stimuli. For example, the IRE1/XBP1s (red in Fig. 1D) and ATF6 (blue in Fig. 1D) genesets show significant overlap, reflecting the coordinated activation of these two pathways in response to ER stress. Furthermore, certain genes such as SOD1, activated downstream of the OSR-associated transcription factor NRF2 50 , separate from the OSR cluster (purple in Fig. 1C-D), reflecting the ability for this gene to be regulated by multiple stress-responsive signaling pathways apart from the OSR 51 . The PERK/ISR target MTHFR is also regulated by other UPR signaling pathways, as well as OSR 52,53 , and is similarly found to separate from the larger cluster of PERK/ISR targets (green in Fig. 1C-D). The overlap of genesets and promiscuity for specific genes to report on multiple pathways highlights the importance of tracking sets of stress-responsive genes for defining overall pathway activation. However, the general clustering of our stress-pathway genesets indicates that this targeted RNAseq assay is capable of tracking changes in stress responsive genes to accurately define activation of specific stress-responsive signaling pathways.

Defining stress-independent HSR and UPR signaling pathways through targeted RNAseq profiling.
We initially validated the ability for our targeted RNAseq assay to report on activation of specific stresssignaling proteostasis pathways using chemical genetic approaches that allow activation of specific pathways independent of stress. First, we defined activation of the HSR-regulated proteostasis genes in HEK293 TREX cells following doxycycline (dox)-dependent activation of a constitutively active cHSF1 17 -the primary transcription factor regulated by the HSR [6][7][8] . In order to define the induction of specific target proteostasis genes in our targeted RNAseq data, we first median-normalized aligned counts per target gene across all treatment conditions. We took average normalized count values across sample replicates and performed a Log 2 transformation, yielding the "Log 2 normalized counts" used for relative expression analysis. To compare chemical genetic and pharmacologic activating conditions versus vehicle control samples, we conducted a correlation analysis of Log 2 normalized counts to yield a line of best fit ( Fig. 2A), reflecting baseline expression levels for the majority of genes not affected by a given treatment. We then calculated the deviation of each target gene for the experimental condition from the line of best fit, or "residual value", which was used to quantify up/downregulation of that gene. (Fig. 2A). Finally, we define pathway activation by plotting the residual values of each gene from this analysis, grouped according to assigned stress-responsive pathway, and monitoring the overall behavior of the geneset (Fig. 2B). This allows us to normalize variability in gene induction across different treatments and minimize challenges associated with lowly expressed genes that show high levels of induction. From this analysis, we demonstrate that dox-dependent cHSF1 activation robustly and selectively activates the entire target HSR-regulated proteostasis program, thus confirming the ability for our targeted RNAseq assay to define activation of this pathway ( Fig. 2B and Table S3). Interestingly, the activation of this pathway is identical to that observed when we perform the same analysis using published RNAseq transcriptional profiles for dox-dependent cHSF1 activation 17 , demonstrating that our RNAseq assay accurately quantifies the induction of HSR-regulated proteostasis target genes ( Fig. S2A-C).
A significant challenge in monitoring activation of stress-responsive signaling pathways is the overlap between closely related pathways. For example, the IRE1/XBP1s UPR pathway induces expression of multiple genes also regulated by the ATF6 UPR signaling pathway, albeit to lower extents 34 . Furthermore, other XBP1s target genes are often induced to lower levels than that observed for ATF6-selective target genes 34 . To define the potential for targeted RNAseq to discern selective activation of these two UPR signaling pathways, we performed this assay in HEK293 DAX cells subjected to stress-independent XBP1s and/or ATF6 activation.
HEK293 DAX cells express dox-inducible XBP1s and trimethoprim (TMP)-regulated DHFR.ATF6, allowing activation of these two transcription factors in the same cell independent of ER stress through administration of dox and/or TMP 34 . As a control, we also monitored gene expression in response to global ER stress induced by treating HEK293 DAX cells with the SERCA inhibitor, thapsigargin (Tg). As expected, Tg treatment showed robust activation of all three UPR signaling pathways (IRE1/XBP1s, ATF6, and PERK/ISR), confirming global UPR activation (Fig. S3A-B). In contrast, TMP-dependent DHFR-ATF6 activation showed significant increases in the ATF6 target geneset, consistent with selective ATF6 activation (Fig. 3A). However, dox-inducible XBP1s increased expression of both the IRE1/XBP1s and ATF6 target genesets, although ATF6 target genes were induced less than that observed following ATF6 activation (Fig. 3B), which is consistent with previous work 34 .
Previous reports indicate that the overlap between XBP1s and ATF6 target gene expression observed following stress-independent activation could be deconvoluted by normalizing the expression of genes to that observed with Tg treatment, providing a way to sensitively define the extent of pathway activation 25 . Performing this normalization shows that TMP-dependent DHFR.ATF6 activation selectively induces expression of ATF6 target genes to levels similar to those observed for Tg-dependent ER stress (Fig. 3D). Importantly, doxdependent XBP1s activation selectively induces expression of IRE1/XBP1s target genes to levels similar to that observed with Tg by this analysis, while only moderately affecting ATF6 target gene expression (Fig. 3E).
This profile is distinct from that observed in cells where XBP1s and ATF6 are co-activated, which shows significantly higher induction of both genesets (Fig. 3F). Importantly, when residual values from our targeted RNAseq analysis are compared to transcriptional changes from whole-transcriptome RNAseq collected from HEK293 DAX cells subjected to XBP1s and/or ATF6 activation, there is a clear correlation between the two data sets, especially in upregulated targets (Fig. S3C-E). These results demonstrate that our targeted RNAseq assay can sensitively deconvolute the complex integration of stress-responsive signaling pathways involved in UPR signaling.

Targeted RNAseq profiling defines stress-pathway activation induced by cellular toxins.
We next used our targeted RNAseq assay to profile activation of stress-responsive signaling pathways induced by chemical toxins including tunicamycin (Tm; an ER stressor that inhibits N-linked glycosylation), the environmental toxin arsenite (As(III)), the mitochondrial ATP synthase inhibitor oligomycin, and the ROS generating compound paraquat (PQ). As predicted, our assay demonstrates that these compounds induce unique activation profiles of different stress-responsive signaling pathways. Consistent with the selective induction of ER stress, Tm treatment activates the three arms of the UPR without globally impacting other stress-responsive signaling pathways (Fig. 4A). In contrast, As(III) induces robust activation of the cytosolic HSR, OSR, and ISR signaling pathways (Fig. 4B), highlighting the promiscuous nature of this toxin for cytosolic proteostasis pathway activation 54 . Oligomycin treatment only significantly activated the ISR geneset, reflecting emerging evidence showing that mitochondrial stress promotes signaling through this pathway ( Fig.   4C) 8,48,49,55 . PQ treatment also showed modest increases in ISR genes, although the entire pathway was not significantly activated (Fig. S4). However, while our genesets report on activation of whole pathways, numerous individual stress-responsive genes from multiple pathways were induced upon these different treatments. For example, the OSR target gene HMOX1 is induced in cells treated with mitochondrial toxins oligomycin and paraquat, although we do not observe induction of other OSR target genes. Since HMOX1 can be regulated by multiple stress-responsive signaling pathways 33 , these results suggest that administration of these toxins induce pleiotropic effects on multiple stress-responsive signaling pathways outside of the four primary proteostasis pathways profiled in our targeted RNAseq platform. Regardless, it is clear that our targeted RNAseq assay does accurately reflect predicted toxin-induced activation of proteostasis pathways, further validating the benefit of this approach for profiling pharmacologic activators of stress-responsive proteostasis pathways.

Defining selectivity of pharmacologic NRF2 activating compounds through Targeted RNAseq transcriptional profiling.
We next employed our targeted RNAseq assay to define the selectivity of two putative NRF2 activating compounds: bardoxolone and the recently described CBR-470-1 (Fig. S5A) 26 . Bardoxolone is an antiinflammatory compound currently in clinical trials for chronic kidney disease. This compound is reported to induce protective benefits through activation of the OSR-associated transcription factor NRF2 22,56 . However, it also covalently modifies hundreds of proteins 57 and displays additional cellular activities including inhibition of the mitochondrial protease LON 58 , suggesting that, apart from NRF2, bardoxolone could also activate other stress-responsive signaling pathways. Interestingly, we show in HEK293T cells that bardoxolone significantly induces expression of the OSR-target gene HMOX1, but not other OSR target genes, suggesting that this compound does not robustly activate the NRF2 transcriptional response in these cells (Fig. 5A). However, this compound does induce both the HSR and the ISR genesets, indicating promiscuous activity for this pharmacologic agent. Furthermore, we see strong upregulation of the ATF6 target gene HSPA5 (also known as BiP), without complete activation of the ATF6 pathway. These results indicate that bardoxolone induces pleiotropic effects on stress-responsive genes outside of NRF2 activation. In contrast, CBR-470-1 showed selective activation of the OSR geneset with no significant induction of other stress pathways, suggesting improved selectivity of CBR-470-1 for OSR activation (Fig. 5B). Consistent with this, qPCR analysis of BAG3 (an HSR target) and HSPA5 shows that bardoxolone promiscuously induces these non-NRF2 target genes, while CBR-470-1 does not (Fig. 5C,D). However, both compounds induce activation of the OSR target gene HMOX1 (Fig. 5E). This result is identical to that observed by our targeted RNAseq analysis (Fig. S5B-D).
These results show that CBR-470-1 shows increased selectivity for OSR activation relative to bardoxolone and demonstrates the utility for our targeted RNAseq assay to profile selectivity of putative OSR activating compounds in clinical development.

Identification of selective HSR activating compounds by targeted RNAseq
Previous high-throughput screening identified numerous compounds, including compounds A3, C1, D1, and F1 ( Fig. S6A), that activate a cell based reporter of the HSR-associated transcription factor HSF1 in HeLa cells 20 .
However, the selectivity of these compounds for the HSR remains to be fully defined. We used our targeted RNAseq assay to define the selectivity for these putative HSF1 activating compounds for specific HSR activation. Our results show that compound A3 strongly induced the HSR geneset (Fig. 6A) to a level comparable to that observed with dox-dependent cHSF1 activation (Fig. 2B). Compounds C1, D1, and F1 also significantly induced the HSR geneset, albeit to a lower extent (Fig. 6B-D). Interestingly, administration of these compounds also induced expression of other stress-responsive genes. This was most evident with A3, which showed robust activation of select ISR and OSR target genes such as ATF3 and HMOX1, respectively, without global activation of these pathways (Fig. 6A). Similar results were observed for the other three compounds to lesser extents (Fig. 6B-D). Interestingly, both ATF3 and HMOX1 have been shown to be transcriptionally induced following stress-independent activation of the HSR-associated transcription factor HSF1 17 , suggesting that their increased expression in response compound treatment could, in part, reflect HSF1 activity.
To further define the selectivity of these HSR activating compounds for the HSR proteostasis transcriptional program in HEK293T cells, we performed whole transcriptome RNAseq ( Table S4). Analysis of the top 100 most significantly altered transcripts in this whole transcriptome RNAseq data demonstrated that compound A3 induced the greatest effects on gene expression, consistent with our targeted RNAseq results (Fig. 6E). Furthermore, performing the same correlation-based geneset analysis used for targeted RNAseq revealed an identical preferential activation of the HSR in this whole transcriptome dataset (Fig. S6B-I).
Interestingly, comparing genes induced by A3 with those induced by a 43 o C heat shock 6 or dox-dependent cHSF1 17 activation demonstrated an overlap of ~100 genes (Fig. 6F), including many classical HSR proteostasis target genes such as HSPA1A, DNAJB1, and CryAB (Table S4). While this supports an A3dependent induction of the HSR, there are a large number of genes upregulated in the whole-transcriptome data that are not upregulated by these other HSR-activating insults. GO analysis reveals that, apart from proteostasis factors, the majority of targets induced by treatment with A3 are involved in RNA polymerase IIdependent transcription (Fig. 6G). This finding is consistent with recent studies indicating that apart from direct transcriptional upregulation, HSF1 may recruit factors that increase the rate of Pol II release from its paused state in transcript promoter regions 6 . Thus, the altered expression of Pol II regulatory factors suggests that A3 could influence HSR activity by targeting RNA Pol II pause-release. While the impact of A3 on RNA Pol II could limit the further development and usage of this compound as a chemical tool to define HSR-dependent regulation of cellular proteostasis, these results demonstrate the potential for our targeted RNAseq assay to define the selectivity of prioritized compounds identified through reporter based HTS for activating specific stress-responsive proteostasis pathways.

CONCLUDING REMARKS
The establishment of highly selective activators of stress-responsive signaling pathways provides unique opportunities to identify new roles for these pathways in regulating cellular physiology and defining their potential for correcting pathologic defects associated with diverse diseases. Here, we establish a mediumthroughput targeted RNAseq assay that reports on the activation of predominant stress-responsive proteostasis pathways such as the HSR, OSR, ISR, and UPR. We demonstrate the potential for this approach to deconvolute the complex integration of stress-responsive signaling pathways in cells treated with chemical genetic or pharmacologic perturbations. Furthermore, we show that this approach is suitable for defining the selectivity of putative activators of different stress-responsive signaling pathways. These results demonstrate that this assay provides new opportunities to improve screening efforts focused on establishing pharmacologic activators of stress-responsive signaling pathways by identifying compounds or classes with high selectivity earlier in the screening pipeline (i.e., after reporter based HTS). Furthermore, this assay can be paired with medicinal chemistry to establish next generation compounds with improved selectivity and/or potency through monitoring activation of specific pathway reporter genesets. While we specifically focus on four primary stressresponsive proteostasis pathways (see Fig. 1A), this assay can be easily adapted to include new or improved genesets reporting on other signaling pathways, further improving the ability for this approach to define compound selectivity at early stages of compound development. Similarly, this assay can be adapted to other organisms, providing a unique opportunity to monitor activation of stress pathways in different tissues and define compound pharmacodynamics. Ultimately, the targeted RNAseq assay and approach described herein will improve the development of pharmacologic activators of stress-responsive signaling.

Cell lines and treatments
Stable cell lines expressing inducible cHSF1, ATF6, and XBP1s were used as previously described to activate cHSF1, ATF6, and XBP1s transcription factors 17,34 . Cells types as listed in Table 2 were cultured in DMEM with 10% FBS, pen/strep, and glutamine at 37C, 5% CO2 in Corning 96-well tissue culture plates. Cells were treated for the indicated durations (Supplementary Table S2) with either chemical genetic activators or pharmacologics solubilized in DMSO, treatments were performed in triplicate .

RNA extraction
RNA was extracted from HEK293T, HEK293 Dax , and HEK293 TREX using the Zymo Research ZR-96 Quick-RNA isolation kit as per manufacturer's instructions. Briefly, cells were rinsed with DMSO prior to lysis, samples were then cleared of particulate matter through centrifugation. The supernatant was then subject to standard column purification steps including an on-column DNase treatment, prior to RNA elution in DNase/RNase-Free water.

Alignment and Expression Analysis
Alignment of targeted RNA-seq reads was performed using the Illumina TruSeq Targeted RNA-seq software using the custom target manifest containing sequences of targeted region sequences. Alignment of wholetranscriptome RNA-seq data was done using DNAstar Lasergene SeqManPro to the GRCh37.p13 human genome reference assembly. Aligned counts from Targeted RNAseq were median normalized and Log2 transformed prior to correlation analysis.

Whole-transcriptome RNA-seq
HEK293T cells were treated with 10uM A3, C1, D1, or F1 for 6 hours prior to RNA isolation using the ZymoPure RNA-mini kit as per manufacturer's instructions, including on-column DNase I treatment to remove contaminating genomic DNA. RNA was quantified by NanoDrop. Conventional RNA-seq was conducted via BGI Americas on the BGI Proprietary platform, providing single-end 50bp reads at 20 million reads per sample.
Each condition was performed in triplicate.

Gene expression correlation network graph and hierarchical clustering
Raw counts data from the triplicate RNAseq experiments for each condition were averaged and Pearson correlation coefficients were calculated for each pair of genes. We created the gene expression correlation graph by representing each gene as a vertex and connecting the vertices for the genes that had correlation coefficients ≥ 0.6. There were a few genes whose expression levels did not correlate with those of any other genes at this level. These genes were connected only to the gene with which they had the highest correlation coefficient to ensure that the network graph was fully connected. The hierarchical clustering of genes by expression pattern shown in Figure 1 was performed using the Euclidean distance between each genes' expression level correlation coefficients with all other genes as the distance metric and single-linkage clustering as the linkage criterion. Thus, two genes that have similar sets of correlation coefficients with all other genes were most likely to cluster together. The network graph and dendrogram in Figure 1 were produced using Mathematica 11.3.

Statistical analysis
Statistical significance of residual values from targeted RNAseq correlation analysis were calculated using standard one-way ANOVA, with the lowest p-values presented. Full ANOVA tables from our targeted RNAseq assay are included in Table S3. qPCR data was analyzed using one-tailed Student's t-test.  Stresses that activate each pathway and specific transcription factors activated downstream of these pathways are also shown.

TABLES
B. Schematic of the general protocol used for our targeted RNAseq assay. Briefly, RNA is isolated from cells following a given treatment. This RNA is then converted into cDNA libraries that are probed using oligos targeted to specific stress-responsive genes (red) for sequencing library generation.
C. Dendrogram of individual target genes from our targeted RNAseq panel (see Table 1) grouped by hierarchical clustering using the Euclidean distance between each genes' expression level correlation coefficients over all treatment conditions (see Table 2 and Table S1). Genes are colored by assignment to specific stress-responsive signaling pathways to report on activation of the HSR (orange), the OSR (purple), the ATF6 UPR signaling pathway (blue), the IRE1/XBP1s UPR signaling pathway (red), the PERK/ISR signaling pathway (green), or other pathways (grey). The asterisks identify SOD1 (purple) and MTHFR (green).
D. Network graph of individual target genes from our targeted RNAseq panel showing the clustering of genes into defined stress-responsive signaling pathways. This graph is derived by representing each gene as a vertex and connecting the vertices for genes whose expression level changes correlate with Pearson R > 0.6. Genes that do not correlate at this level with any other genes are connected only to the gene with which they have the highest correlation coefficient. Pathways are colored using the same scheme described above in Fig. 1C. SOD1 and MTHFR are identified by name. vehicle. Calculation of residuals was performed as described in Fig. 2A. Statistics were calculated using one-way ANOVA. Significance shown reflects comparison to "Other" target transcript set. **p<0.01, ***p<0.001, ****p<0.0001. See Table S3 for full ANOVA table. Calculation of residuals was performed as described in Fig. 2A. Statistics were calculated using one-way ANOVA. Significance shown reflects comparison to "Other" target transcript set. ****p<0.0001. See Table S3 for full ANOVA table.
B. Plot showing residuals calculated by comparing the expression of our panel of stress-responsive genes between HEK293 DAX cells following treatment with dox (1 µg/mL 4 h; activates dox-inducible XBP1s) or vehicle.
Calculation of residuals was performed as described in Fig. 2A. Statistics were calculated using one-way ANOVA. Significance shown reflects comparison to "Other" target transcript set. ****p<0.0001. See Table S3 for full ANOVA table.
C. Plot showing residuals calculated by comparing the expression of our panel of stress-responsive genes between HEK293 DAX cells following treatment with both trimethoprim (10 µM, 4 h; activates DHFR.ATF6) and dox (1 µg/mL 4 h; activates dox-inducible XBP1s) or vehicle. Calculation of residuals was performed as described in Fig. 2A. Statistics were calculated using one-way ANOVA. Significance shown reflects comparison to "Other" target transcript set. ****p<0.0001. See Table S3 for full ANOVA table. D. Graph showing normalized residuals for genesets regulated by the ATF6 (blue), XBP1s (red) or PERK (green) UPR signaling pathways in HEK293 DAX cells following treatment with TMP (10 µM, 4 h; activates DHFR.ATF6). The residuals for each gene were normalized to those observed for thapsigargin (Tg)-induced ER stress in HEK293 DAX cells (Fig. S3A,B). Normalized data was subjected to ROUT outlier testing. Statistics from one-way ANOVA **p<0.01, ***p<0.001.

E.
Graph showing normalized residuals for genesets regulated by the ATF6 (blue), XBP1s (red) or PERK (green) UPR signaling pathways in HEK293 DAX cells following treatment with dox (1µg/mL, 4 h; activates doxinducible XBP1s). The residuals for each gene were normalized to those observed for thapsigargin (Tg)induced ER stress in HEK293 DAX cells (Fig. S3A,B). Normalized data was subjected to ROUT outlier testing.
Normalized data was subjected to ROUT outlier testing. Statistics from one-way ANOVA. *p<0.05, **p<0.01.  Fig. 2A. Genes are grouped by target stress-responsive signaling pathway. Statistics were calculated using one-way ANOVA, Significance shown reflects comparison to "Other" target transcript set. ***p<0.001, ****p<0.0001. See Table S3 for full ANOVA table. B. Plot showing residuals calculated by comparing the expression of our stress-responsive gene panel between HEK293T cells following treatment with arsenite (As(III); 25 µM, 16 h) or vehicle. Calculation of residuals was performed as described in Fig. 2A. Genes are grouped by target stress-responsive signaling pathway. Statistics were calculated using one-way ANOVA, Significance shown reflects comparison to "Other" target transcript set. **p<0.01, ***p<0.001. See Table S3 for full ANOVA table. C. Plot showing residuals calculated by comparing the expression of our stress-responsive gene panel between HEK293T cells following treatment with oligomycin A (Oligo; 100 nM, 24 h) or vehicle. Calculation of residuals was performed as described in Fig. 2A. Genes are grouped by target stress-responsive signaling pathway. Statistics were calculated using one-way ANOVA, Significance shown reflects comparison to "Other" target transcript set. **p<0.01. See Table S3 for full ANOVA table.  Fig. 2A. Genes are grouped by target stress-responsive signaling pathway.
Statistics were calculated using one-way ANOVA, Significance shown reflects comparison to "Other" target transcript set. **p<0.01, ****p<0.0001. See Table S3 for full ANOVA table. B. Plot showing residuals calculated by comparing the expression of our stress-responsive gene panel between HEK293T cells treated with CBR-470-1 (10 µM, 24 h) or vehicle. Calculation of residuals was performed as described in Fig. 2A. Genes are grouped by target stress-responsive signaling pathway.
Statistics were calculated using one-way ANOVA, Significance shown reflects comparison to "Other" target transcript set. **p<0.01. See Table S3 for full ANOVA table.  Statistics were calculated using one-way ANOVA. Significance shown reflects comparison to "Other" target transcript set. ****p<0.0001. See Table S3 for full ANOVA table.