TOR and heat shock pathways control peroxisome creation during protein stress

Kar2, the yeast homolog of mammalian BiP, belongs to the Hsp70 family of chaperone proteins residing in the ER lumen. During protein misfolding stress in the ER lumen, Kar2 binds to unfolded proteins thereby dissociating from its binding partner, Ire1²⁴. Inactivation of Kar2 induces the unfolded protein response (UPR, monitored by induction of UPRE-GFP) (Supplementary Fig. 1a, b), as well as the cytosolic heat shock response (HSR, monitored by induction of HSE-GFP) (Supplementary Fig. 1c, d), when compared to WT cells at 37 °C²⁴. To understand how loss of Kar2 function is connected to peroxisome homeostasis, we investigated the effects of UPR and HSR activation on peroxisome proliferation.

After confirming the induction of UPRE-GFP by tunicamycin (Supplementary Fig. 1e, f) (Tm), a nucleoside antibiotic that interferes with N-glycosylation of proteins in the ER lumen, we tested its effect on peroxisome number by counting Pex3-GFP and GFP-ePTS1 puncta. Pex3 is a peroxisomal membrane protein (PMP), whereas GFP-ePTS1 is imported into the peroxisome lumen. Tunicamycin (non-lethal doses, Supplementary Fig. 1i) significantly increased the number of peroxisomes as visualized by both markers (Fig. 1c–f), as well as the cumulative area occupied by all peroxisomes per cell (Supplementary Fig. 1j, k), within 2–5 h of treatment (Supplementary Fig. 1n). We also verified that all GFP-ePTS1 and Pex3-GFP puncta represent bona fide peroxisomes in both DMSO and tunicamycin-treatment conditions (Supplementary Fig. 1l, m), and that ER stress induced by non-lethal doses of DTT also increased peroxisome number (Fig. 1g and Supplementary Fig. 1i). We confirmed a recent report, which suggested that peroxisome numbers (monitored indirectly by total Pex11-mNeonGreen fluorescence)⁹ increase as an adaptation to supply additional acetyl-CoA via peroxisomal β-oxidation to support mitochondrial respiration during ER stress. To directly explore this idea, we examined whether ER stress induced peroxisomal β-oxidation enzymes, which are repressed during growth in glucose. Endogenous levels of GFP-tagged 3-ketoacyl-CoA thiolase (Pot1)⁸, which catalyzes the final step in peroxisomal β-oxidation, increased following tunicamycin treatment (Fig. 1h, i). Collectively, our observations indicate that ER stress resulting from misfolded proteins, caused either by loss of Kar2 function or by tunicamycin, results in peroxisome proliferation (Fig. 1j and Supplementary Fig. 1).

Analogous to the UPR in the ER, the heat shock response pathway is activated in response to the presence of misfolded and aggregated proteins in the cytosol. Signaling through this pathway is primarily mediated by Hsf1, which binds to the Heat Shock Elements (HSEs) in gene promoters and activates transcription of chaperones and components of the Ubiquitin-proteasome system (UPS) involved in protein degradation³³. We induced the heat shock response pathway by deleting SSA1, which encodes the Hsp70 chaperone that binds and represses Hsf1³⁴. Under non-stress conditions, Ssa1 and its paralog Ssa2 maintain Hsf1 in an inactive state; however, during stress, unfolded or misfolded proteins titrate the Hsp70s away from Hsf1, thereby allowing transcriptional activation at HSEs. The ssa1∆ cells exhibited not only elevated HSE-GFP levels (Fig. 3e and Supplementary Fig. 3b), but also a marked increase in peroxisome number compared to WT (Fig. 3b, c). The median number of peroxisomes in ssa1∆ was 7–8, versus 3 in WT, with ~70% of ssa1∆ cells containing >5 peroxisomes, compared to ~5% in WT. Loss of Ssa1 caused an increase in the number of peroxisomes, comparable to that of tunicamycin-treated wild-type cells (Fig. 3f, g). Surprisingly, tunicamycin further increased the number of peroxisomes in ssa1∆ cells, with the median rising to 9–12 after tunicamycin treatment, compared to 6-7 in DMSO-treated cells (Fig. 3f, g). These results suggest that tunicamycin and Hsf1 activation may act synergistically in stimulating peroxisome proliferation.

The adaptive response to ER stress occurs via activation of SNF1-, HOG- and RTG-signaling pathways and inactivation of TOR1 and PKA^27,36–40; however, among these, the specific pathways that mediate the signaling response to trigger peroxisome proliferation in response to tunicamycin remain unclear (Fig. 4a). We systematically abrogated signaling through each pathway and investigated the effect of tunicamycin on peroxisome proliferation in mutants compared to WT using our GLM analysis (Fig. 4). Given published findings on the involvement of these pathways in yeast peroxisome proliferation and regulation of PEX gene expression, we were surprised to discover that blocking the retrograde signaling (RTG) (Fig. 4b, g)^9,11, SNF1^4,8,41 (Fig. 4c, g), or HOG pathways⁷ (Fig. 4c, g) did not significantly impair peroxisome proliferation upon tunicamycin treatment. Similarly, the loss of both Msn2 and Msn4 (Fig. 4d, g), which mediate the transcriptional response upon activation of the HOG pathway, as well as the heat shock response pathway⁷, also did not prevent tunicamycin-induced peroxisome proliferation. Furthermore, we also found that cells lacking Mig1, the peroxisome biogenesis repressor that gets phosphorylated and subsequently inactivated by active Snf1, or Mig2, which shares overlapping functions with Mig1^42–44, showed a comparable number of peroxisomes to WT, after both DMSO as well as tunicamycin treatments (Fig. 4e, g). Consistent with our findings that peroxisome proliferation in ER stress is independent of Snf1, we observed that Snf1 was not activated after tunicamycin treatment under our assay conditions; however, Hog1 was activated after 4 h (Supplementary Fig. 4).

Response to ER stress also involves the inhibition of PKA⁴⁵ and TOR1³⁸ (Supplementary Fig. 4). We tested the role of PKA by using a strain overexpressing TPK1 (TPK1oe) using a multicopy plasmid⁴⁶, which constitutively activates PKA, relative to a control strain carrying an empty vector. Our GLM analysis indicated that the relative increase in peroxisome number in response to tunicamycin treatment was reduced in TPK1oe cells compared to that in cells carrying the empty vector (Fig. 4f, g). This result suggests that high levels of PKA activity attenuate the tunicamycin-induced peroxisome proliferation.

We tested the significance of peroxisomes in adaptation to ER stress by measuring the survival of WT and pex3∆ cells in rich media containing 2% or 0.5% dextrose (where respiration is more prevalent) after overnight (~14 h) treatment with a low dose of 0.5 µg/ml tunicamycin followed by a high dose of 20 µg/ml tunicamycin for 6 h (Fig. 5f and Supplementary Fig. 5b). Corroborating similar findings in S. cerevisiae⁹, K. phaffii cells lacking Pex3, which is essential for peroxisome biogenesis², consistently showed <50% survival in the tunicamycin-viability assay, which was reduced to <10% under reduced glucose availability compared to WT cells, underscoring the significance of peroxisome biogenesis in the adaptation.

Peroxisome contribution to cellular adaptation during ER stress was analyzed by comparing the survival of human primary fibroblasts derived from a Zellweger Spectrum Disorder patient carrying homozygous PEX16 R176* mutations (hereafter referred to as pex16KO) to WT⁴⁷. Fibroblasts were first pretreated for 12 h with primary doses of tunicamycin (range 0 to 20 µg/ml), then challenged with a secondary tunicamycin treatment (20 µg/ml) or DMSO for 12 h and then assayed for viability. The sensitivity to primary tunicamycin treatment alone was more blunted in pex16KO cells compared to WT (Fig. 6c, "DMSO" panel), possibly the result of higher metabolic activity in the latter. When primary treatment was followed with secondary tunicamycin challenge, the prior exposure to tunicamycin conferred a protective effect against subsequent stress in both cell lines, but this benefit was diminished in pex16KO cells (Fig. 6c, "Tm" panel). A GLM analysis of this data revealed a significant interaction between genotype, primary tunicamycin dose, and secondary tunicamycin challenge (likelihood ratio test P value = 0.004), indicating that prior exposure to tunicamycin confers a greater adaptation to secondary challenge in WT versus pex16KO cells. These findings suggest that peroxisome biogenesis supports adaptive resilience during ER stress, although the residual benefit observed in pex16KO cells suggests the possibility of additional compensatory pathways.

We next tested whether the increase in peroxisome number during ER stress depended on peroxisome division or on their de novo synthesis. Vps1 is the primary mediator of peroxisome fission in cells grown in the presence of glucose⁴⁹, and its deletion resulted in only a partial reduction in tunicamycin-driven peroxisome proliferation compared to WT cells (Fig. 7d, e). We also quantified this tunicamycin response in dnm1∆ single and vps1∆ dnm1∆ double mutants. While the loss of DNM1 alone had an insignificant effect, the absence of both DNM1 and VPS1 significantly, but still only partially, impaired peroxisome proliferation after tunicamycin treatment (Fig. 7d, e), suggesting additional peroxisome formation occurs by the de novo pathway. In DMSO-treated cells, the median number of peroxisomes was 1 in both vps1∆ and vps1∆ dnm1∆ mutants, while tunicamycin-treated cells had a median number of 2 peroxisomes per cell. This suggests that tunicamycin treatment led to the de novo formation of at least one additional peroxisome. In comparison, WT cells had a median number of 3 peroxisomes in DMSO treatment, and 4.5 peroxisomes after tunicamycin treatment. Thus, of the ~1.5 additional peroxisomes formed in WT in response to tunicamycin, at least one could be attributed to de novo biogenesis. Overall, our results demonstrate that de novo biogenesis accounts for approximately two-thirds of the peroxisome formation in response to ER stress.

Further evidence supporting the role of de novo biogenesis in increased peroxisome production during ER stress comes from analyzing the proportion of cells containing peroxisomes in the inheritance mutants inp2∆ and inp1∆ (Fig. 7f, g). In the absence of Inp2, most daughter cells fail to inherit peroxisomes from mother cells and therefore rely exclusively on de novo biogenesis to produce their first peroxisome^50,51. Conversely, most mother cells fail to retain peroxisomes in the absence of Inp1⁵², resulting in the total transfer of peroxisomes to their daughter cells. We observed ~45% of inp2∆ cells lacked peroxisomes under DMSO-treated conditions; however, tunicamycin treatment reduced this proportion to ~22% (Fig. 7g), strongly suggesting that ER stress stimulated de novo biogenesis. Corroborating this observation, we found ~25% of inp1∆ cells lacked peroxisomes after tunicamycin treatment, compared to ~40% under DMSO treatment conditions (Fig. 7g). Since de novo peroxisome biogenesis in S. cerevisiae is slow, requiring ~4–5 h^50,53, and our analyses were performed within this same timeframe, most of the new peroxisomes observed in response to ER stress are likely due to de novo biogenesis, with the remaining, comprising about one-third, likely derived from the growth and division.

Lipid droplets play important roles in energy metabolism and lipid homeostasis by serving as lipid storage sites in cells. Moreover, not only are lipid droplets and peroxisomes generated from a common site in the ER membrane, but the two organelles also maintain crosstalk through contact sites². Furthermore, ER stress stimulates lipid droplet formation⁵⁶, as well as microlipophagy⁵⁷, and one possibility is that lipid droplets could serve as the source of fatty acids for β-oxidation by peroxisomes. We examined how peroxisome proliferation in response to ER stress is coordinated with lipid droplet homeostasis. We found that upon impairing lipid droplet formation by using the quadruple deletion (∆QD) strain, are1∆ are2∆ lro1∆ dga1∆⁵⁸, peroxisome number was comparable to WT (Fig. 8d, e). The absence of lipid droplets also did not abrogate the tunicamycin-driven increase in peroxisome number (Figs. 4g and 8e) indicating peroxisome proliferation occurs independent of lipid droplet biogenesis.

Yeast strains were generated using lithium acetate-PEG-mediated transformation and are listed in Supplementary Data 1. Strains expressing P_TDH3-GFP-ePTS1-T_PGK1 (i.e., GFP-ePTS1), 4×UPRE-P_CYC1-GFP-T_ACT1 (i.e., UPRE-GFP) or 4×HSE-P_CYC1-GFP-T_ACT1 (i.e., HSE-GFP) were generated by transforming linear constructs obtained by double digesting the plasmids pNS31, pNS32 or pLT1, respectively, with PciI and NotI. All three constructs were targeted after the STOP codon of the UBC9 ORF. Except where noted, parental knock-out strains were obtained from the MATα deletion library (Invitrogen)⁶⁸. SNF1, INP2, and DNM1 genes were knocked out by transforming the corresponding deletion constructs, which in turn, were generated by fusing hphMX6 or Zeo cassettes on either side with ~600-1000 bp sequences upstream and downstream of the respective ORFs. Appropriate replacement of these ORFs by the selection marker was confirmed by primers shown in Supplementary Data 2. 4×UPRE-P_CYC1-GFP-T_ACT1 was amplified from pRH1209⁶⁹. For generating strains overexpressing TPK1 (TPK1oe), plasmid pPHY2056⁴⁶ was transformed in yeast; the corresponding control strain was made by transforming the empty vector (EV), pRS426 in yeast.

Except when noted, S. cerevisiae log phase cells grown in synthetic defined [SD; 6.7 g/L Yeast nitrogen base (YNB without amino acids, BD/Difco; 291920) + 0.79 g CSM (Sunrise Science Products; 1001-100) + 2% dextrose] at 30 °C and shaken at 250 RPM were used. For experiments with HSE-GFP and kar2-159, cultures were grown at 25 °C and 250 RPM. RTG pathway mutants are glutamate auxotrophs, and hence those mutants and their corresponding WT were grown in SD + 0.02% glutamic acid. Cultures were resuspended in fresh media at an OD ~0.05–0.1 prior to performing temperature shifts (to 30 °C or 37 °C for 3 h as indicated in figure legends) or treatments with tunicamycin (Sigma, T7765; 1 mg/ml stock in DMSO) or rapamycin (Sigma, R0395; 1 mg/ml stock in DMSO) or phytosphingosine (Sigma, Cas No. 554-62-1; 10 mM stock in DMSO) or DMSO or DTT (Roche, 10708984001; 1 M stock in H₂O). TPK1 overexpression experiments were performed in SD-Ura media to maintain selection for the plasmids. For all the experiments with K. phaffii, cell growth and treatments were done in YPD media (10 g/L Yeast extract + 20 g/L Bactopeptone + 2% Dextrose). SD-Ino media was prepared using YNB without inositol (US Biological Life Sciences; Y2030-01). The effect of tunicamycin on cell viability was measured by quantifying the percentage of cells stained by propidium iodide (PI). For this, yeast cultures were treated with 4 µg/ml PI (Sigma, P4170; 1 mg/ml stock in DMSO) for 20 min at room temperature, followed by 1× wash and resuspension in media for subsequent imaging.

Cells were imaged using a Zeiss Plan Apochromat 63×/1.4 Oil DIC objective mounted on an Axioskop 2 mot plus microscope (Zeiss) equipped with Axio Cam HRm camera and HBO 100 mercury lamp. For 3D imaging to count the number of peroxisomes per cell, Z stacks consisting of 8 slices (for S. cerevisiae) or 6 slices (for K. phaffii) were acquired with a Z spacing of 1 µm. For UPRE-GFP, HSE-GFP, Pot1-GFP and PI staining only a single slice at the focal plane was imaged. Exposure times were kept identical across all experiments for each fluorescent marker, whereas DIC images were acquired with auto exposure.

Western blotting was performed from cell extracts generated after TCA-precipitation using the following antibodies: Phospho-S6 Ribosomal Protein (Ser235/236) (D57.2.2E) XP Rabbit mAb, Cell Signaling Technology (#4858, 1:2000); Phospho-AMPKα (Thr172) (40H9) Rabbit mAb, Cell Signaling Technology (#2535, 1:2000); Phospho-p38 MAPK (Thr180/Tyr182) (D3F9) XP Rabbit mAb, Cell Signaling Technology (#4511, 1:2000); Phospho-PKA Substrate (RRXS*/T*) (100G7E) Rabbit mAb, Cell Signaling Technology (#9624, 1:2000); Anti-ScActin Rabbit pAb²⁰, Gift from Michael Yaffe (1:5000); Goat Anti-Rabbit IgG (H + L)-HRP Conjugate, Bio-RAD; (#1721019, 1:5000).

Human primary fibroblasts were cultured at 37 °C with 5% CO₂ in high-glucose Dulbecco's Modified Eagle Medium (DMEM; high glucose with sodium bicarbonate and non-essential amino acids; Gibco) supplemented with 15% (v/v) heat-inactivated fetal bovine serum (FBS; VWR). Primary human fibroblasts were purchased from Coriell Institute (Cat. No. GM06231 and GM00969). All experiments were conducted using fibroblasts with fewer than ten passages. For imaging, fibroblasts were seeded in glass-bottom 96-well plates (Cellvis) and treated with the indicated concentrations of tunicamycin (Sigma), rapamycin (Sigma), or 4-phenyl butyric acid (4-PBA) (Sigma). After a 12 h incubation, cells were fixed with formaldehyde and stained with NucBlue Stain mixed with ProLong Glass Antifade Mountant (Thermo Fisher Scientific, P36985), PMP70-AF488 [anti-PMP70 antibody, Abcam (EPR5614; 1:1000); donkey anti-Rabbit IgG (H + L)-Alexa Fluor 488, Thermo Fisher Scientific (A-21206, 1:2000)], and Phalloidin-AF750 (Thermo Fisher Scientific, A30105, 165 nM working conc.). Each treatment condition included three biological replicates. For imaging, 6–8 3D confocal z-stacks per well were acquired using a Zeiss 980 confocal microscope equipped with a 63×/1.2 NA water immersion objective lens and NIR detector. Each z-stack had dimensions of 4096 × 4096 × 31 voxels, corresponding to a physical volume of approximately 210 × 210 × 6 µm³.

To quantify peroxisome abundance in individual human fibroblasts, we developed a custom 3D image processing pipeline implemented in Python (Clarifi3D), optimized for high-resolution confocal z-stacks and using NVIDIA A100 GPUs on the Seattle Children's Research Institute High-Performance Compute Cluster, Sasquatch. Raw images (4096 × 4096 × 31 voxels; 0.0509 × 0.0509 × 0.2 µm voxel size) were loaded and preprocessed using GPU-accelerated normalization and filtering routines. Cell and nucleus segmentation were performed using pre-trained 3D U-Net models applied to the phalloidin and DAPI channels, respectively, followed by seed-based Delaunay watershed segmentation. Peroxisomes were segmented from PMP70-AF488 images using a classical spot detection pipeline comprising 3D Laplacian-of-Gaussian morphological filtering, global Otsu thresholding, and seeded 3D Delaunay watershed segmentation. For each field-of-view, the numerical peroxisome density was calculated as the number of segmented peroxisomes divided by the total cytoplasmic volume of fully segmented cells. Results from six to eight stacks per well were averaged across three biological replicates per condition.

Human primary dermal fibroblasts from a healthy donor (GM00969; wild-type, WT) and a ZSD-patient cell line for PEX16 (GM06231; pex16KO) line were cultured in DMEM (high glucose, with sodium bicarbonate and non-essential amino acids; Gibco) supplemented with 15% FBS. Cells were maintained at 37 °C in a humidified atmosphere containing 5% CO₂. For drug treatments, cells were seeded at 2000 cells per well in black-wall, clear-bottom 384-well plates (Costar) in 40 µL media per well and allowed to adhere overnight. Drug treatments were performed in a randomized matrix format using tunicamycin (dissolved in DMSO) in a 1:2 serial dilution. Tunicamycin was titrated to a top concentration of 20 µg/mL. The viability assay was conducted on two independent 384-well plates, with 12 biological replicates per concentration on each plate. After 12 h of drug exposure, half of the wells received a secondary tunicamycin challenge at 20 µg/mL to induce acute ER stress, while the other half received a DMSO mock treatment. Viability was assessed 12 h later using CellTiter-Blue (Promega) per the manufacturer's instructions. CellTiter-Blue fluorescence was quantified using a Spectromax I3 automated plate imager (Molecular Devices). Raw fluorescence intensity was extracted from each well.

Figure images for depicting the number of peroxisomes in cells were made by generating maximum intensity Z projections for the peroxisome channel (GFP) and average intensity Z projection for cell outline (DIC) channel. Figure images for depicting UPRE-GFP, Pot1-GFP, and HSE-GFP were made with single-Z slice images. Identical brightness and contrast settings were applied across all images within a figure panel to ensure fair comparisons, using FIJI/ImageJ software. Scale bars in all figures represent 5 µm, except for the human fibroblasts, where a 50 µm scale bar was used. Quantification for the number of peroxisomes per cell (for GFP-ePTS1 and Pex3-GFP) as well as for the intensity of GFP signal per cell (for Pot1-GFP, UPRE-GFP and HSE-GFP) was performed using the perox-per-cell software²³ using a minimum peroxisome area threshold of 1 pixel and the following values for the software's peroxisome detection sensitivity: 0.0064 for GFP-ePTS1, Pot1-GFP, 4xUPRE-GFP, and 4xHSE-GFP; 0.003 for Pex3-GFP. For quantification of Pot1-GFP, UPRE-GFP and HSE-GFP levels, the values for GFP signal intensities for each cell were divided by 10,000× the value of their respective cell areas and to obtain the "normalized GFP signal" used in the plots. GraphPad Prism (version 10.4.0, Boston, Massachusetts) and R (version 4.4.0) were used to generate the graphs.

No statistical method was used to predetermine sample size. Except for experiments involving kar2-159 and hac1∆ (rationale explained below), no data were excluded from the analyses. The experiments were not randomized except for experiments in Fig. 6c, where for the human fibroblasts grown in 384-well plates, drug treatments were performed in a randomized manner. The investigators were not blinded to allocation during experiments and outcome assessment. All yeast experiments were done twice, in separate experimental batches and all numerical data from these repeats organized as batches are available in Supplementary Data 2. A batch consisted of a group of 1–4 mutants that were tested together with a common WT control. Two different batches for a mutant thus represented two independent treatments, each with at least 100 cells. While the graphs in the figures depict data from a single repeat/batch, GLM analyses were performed on data pooled across batches and thus included both repeats, to rule out discrepancies due to batch-specific effects. In the box plots, the bound of the box indicates the interquartile range (IQR), whereas the center line indicates the medians. The whiskers in the box plots represent 1–25 and 75–99 percentiles of the data, whereas the top and bottom 1 percentiles are represented as dots outside these bounds. In the histograms, the X axis was cropped mostly at 15 or sometimes at 20 peroxisomes/cell, since any higher bin values accounted for <1% of the frequencies; however, these higher data points were not excluded from analyses, and all numerical data are included in Source Data/Supplementary Data 2. We observed significant cell death for kar2-159 ts mutant cells at 37 °C, and hac1∆ after tunicamycin treatment, therefore we manually excluded dead cells from the histograms for those experiments. For this, any atypical cell that showed shriveling, the presence of large vacuoles, or obvious damage and abnormal shape as visualized on the DIC channel was marked as dead. Our empirical observations indicated that almost all cells with zero peroxisomes were dead or severely damaged, except in inp1∆ and inp2∆, where the lack of peroxisomes can be attributed to inheritance defects. We also did not exclude dead cells from the data for GLM analyses (see below), hence many of the mutants that appear to have significant changes in zero peroxisome counts are those where we observed more death after tunicamycin treatment (viz. hac1∆, hog1∆). Using GraphPad Prism, after confirming that the data did not fit normal distribution using D'Agostino & Pearson, Anderson–Darling, Shapiro–Wilk, and Kolmogorov–Smirnov tests, we used the Two-tailed Mann–Whitney test and computed a P value for each replicate to evaluate statistical significance between control and treatment groups. An unpaired t test with Welch's correction was used to evaluate statistical significance for comparison of peroxisome area and cumulative peroxisome area per cell (n = 3 experiments). Statistical significance for peroxisomal numerical density in human fibroblasts was evaluated using a nested t test (n = 3 experiments).

We used generalized linear modeling (GLM) in R (version 4.4.0) to model the distribution of peroxisome (PO) counts in cells across treatment conditions and strains. Due to the presence of a high number of cells with zero counts in some strains, we modeled distributions with hurdle models³⁵ using the "pscl" R package (version 1.5.9), where zero and non-zero (positive) peroxisome counts are captured using separate functions. A binomial distribution was used to model zero counts, and a negative binomial distribution was used to model non-zero counts. We found that, based on the Akaike Information Criterion, hurdle models fit the data better than zero-inflation and simple negative binomial models.

To compare a mutant strain's tunicamycin response to WT, we first combined that strain's peroxisome count data for both the DMSO and tunicamycin treatment conditions with the full set of WT count data collected across all experiments and strains. This included data from both DMSO- and tunicamycin-treated WT cells. We then used GLM to fit the data using a model formula that included strain, treatment, experimental batch, and an interaction term to quantify strain-specific responses to tunicamycin. 1 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\rm{FULL}}}}\; {{{\rm{MODEL}}}}\!\!:{PO\; count} \sim \,{Strain}+{Treatment}+{Batch}+{Strain}\!\!:{Treatment}$$\end{document} FULL MODEL P O c o u n t S t r a i n T r e a t m e n t B a t c h S t r a i n T r e a t m e n t : ~ + + + :

We tested whether the mutant strain's response to tunicamycin was significantly different than WT by performing a likelihood ratio test where GLM results using the following reduced model were compared to those from the full model. 2 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\rm{REDUCED\; MODEL}}}}\!\!:{PO\; count} \sim {Strain}+{Treatment}+{Batch}$$\end{document} REDUCED MODEL P O c o u n t S t r a i n T r e a t m e n t B a t c h : ~ + +

P values collected from these tests across strains were then adjusted for the false discovery rate using the Benjamini–Hochberg method⁷⁰. We considered mutant strains with P values < 0.05 as those with tunicamycin responses that were significantly different from WT.

To determine whether a mutant strain's response to tunicamycin was stronger or weaker compared to WT, we examined the value of the coefficients for the Strain:Treatment interaction term in the mutant strain's full GLM model. Because hurdle models have two modeling components, there are two interaction term coefficients to consider: one from the zero-count model and one from the non-zero count model. For heatmap visualization, the coefficients from a given hurdle modeling component were z-score normalized across the values obtained from all strains (i.e., normalized by column for the heatmaps). For consistency, positive values in these visualizations indicate that the modeling component predicts a stronger tunicamycin response in a mutant strain compared to WT. That is, it predicts a response resulting in higher peroxisome counts. As mentioned earlier, we observed a significantly large proportion of dead cells in hac1∆ and hog1∆ strains following tunicamycin treatment. The zero-count model interaction term coefficients in the GLM analyses on those strains reflected this, resulting in higher predicted zero-count instances. To ensure that the presence of dead cells did not bias our analyses, we made inferences on the relative effects of tunicamycin-induced peroxisome proliferation for mutants vs. WT based only on the interaction term coefficients of the positive-count model.

To quantitatively characterize the contribution of peroxisomes to the adaptive response of human fibroblasts to ER stress, we performed a GLM analysis on the viability assay data from WT and pex16KO human fibroblast cells described above ("Viability assay in human fibroblasts"). Using the CellTiter-Blue fluorescence readout as an indicator of viability, we used R to fit our data using the following formula.3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Viability} \sim \,{Genotype} \, * \,{PrimaryTmConcentration}\, * \,{SecondaryTreatment} +{Row}+{Column}+{Timepoint}$$\end{document}Viability~Genotype*PrimaryTmConcentration*SecondaryTreatment+Row+Column+Timepoint

Here, PrimaryTmConcentration represents the concentration of tunicamycin used to pretreat cells prior to secondary challenge. The SecondaryTreatment term indicates the secondary challenge applied to the cells (DMSO or tunicamycin). The Row and Column terms were included to account for spatial artifacts on experimental plates. Two post-challenge measurements were performed, one at 12 h and one at 12.5 h, and the Timepoint term is included to quantify viability differences between those two timepoints.

To determine whether there was a difference between strains in the way that tunicamycin pretreatment affected viability following secondary challenge, we performed a likelihood ratio test comparing results from the full model formula above to results from a reduced model where the interaction term Genotype:PrimaryTmConcentration:SecondaryTreatment was removed. This is the term in the model that captures how the strains respond to the secondary treatment stress with respect to the primary concentration of tunicamycin they received. The value of this term's coefficient indicates the degree to which primary tunicamycin treatment concentration alters WT cells' viability following secondary tunicamycin challenge compared to pex16KO cells. For scatterplots used to visualize differences in strain responses to secondary challenge, we used our fitted GLM to correct raw experimental fluorescence values for spatial artifacts on experimental plates, then computed the mean values across replicates. All mean values for a specific combination of genotype and secondary treatment were then normalized by subtracting the value observed at a primary tunicamycin concentration of 0 ng/mL. Trendlines in plots were computed using linear regression.

Further information on research design is available in the linked to this article. Nature Portfolio Reporting Summary