Simultaneous proteome localization and turnover analysis reveals spatiotemporal features of protein homeostasis disruptions

Abstract

The functions of proteins depend on their spatial and temporal distributions, which are not directly measured by static protein abundance. Under endoplasmic reticulum (ER) stress, the unfolded protein response (UPR) pathway remediates proteostasis in part by altering the turnover kinetics and spatial distribution of proteins. A global view of these spatiotemporal changes has yet to emerge and it is unknown how they affect different cellular compartments and pathways. Here we describe a mass spectrometry-based proteomics strategy and data analysis pipeline, termed Simultaneous Proteome Localization and Turnover (SPLAT), to measure concurrently the changes in protein turnover and subcellular distribution in the same experiment. Investigating two common UPR models of thapsigargin and tunicamycin challenge in human AC16 cells, we find that the changes in protein turnover kinetics during UPR varies across subcellular localizations, with overall slowdown but an acceleration in endoplasmic reticulum and Golgi proteins involved in stress response. In parallel, the spatial proteomics component of the experiment revealed an externalization of amino acid transporters and ion channels under UPR, as well as the migration of RNA-binding proteins toward an endosome co-sedimenting compartment. The SPLAT experimental design classifies heavy and light SILAC labeled proteins separately, allowing the observation of differential localization of new and old protein pools and capturing a partition of newly synthesized EGFR and ITGAV to the ER under stress that suggests protein trafficking disruptions. Finally, application of SPLAT toward human induced pluripotent stem cell derived cardiomyocytes (iPSC-CM) exposed to the cancer drug carfilzomib, identified a selective disruption of proteostasis in sarcomeric proteins as a potential mechanism of carfilzomib-mediated cardiotoxicity. Taken together, this study provides a global view into the spatiotemporal dynamics of human cardiac cells and demonstrates a method for inferring the coordinations between spatial and temporal proteome regulations in stress and drug response.

Keywords: Unfolded protein response; mass spectrometry; protein turnover; spatial proteomics; spatiotemporal dynamics; stress response; subcellular localization.

Figure 1.

Overview of the SPLAT strategy. A. Experimental workflow. Control, thapsigargin-treated, and tunicamycin-treated human A16 cardiomyocytes were labeled with ¹³C₆¹⁵N₂ L-Lysine and ¹³C₆¹⁵N₄ L-Arginine dynamic SILAC labels. For each condition, 3 biological replicate SPLAT experiments were performed (n=3). After 16 hours, cells were harvested and mechanically disrupted, followed by differential ultracentrifugation steps to pellet proteins across cellular compartments. Proteins from the ultracentrifugation fractions were digested and labeled using tandem mass tag (TMT) followed by mass spectrometry. B. Dynamic SILAC labeling allowed differentiation of pre-existing (unlabeled, i.e., SILAC light) and post-labeling (heavy lysine or arginine, i.e., +R[10.0083]) synthesized peptides at 16 hours. The light and heavy peptides were isolated for fragmentation separately to allow the protein sedimentation profiles containing spatial information to be discerned from TMT channel intensities. C. Computational workflow. Mass spectrometry raw data were converted to mzML format to identify peptides using a database search engine. The mass spectra and identification output were processed using RIANA (left) to quantitate the time dependent change in SILAC labeling intensities and determine the protein half-life, and using pyTMT (right) to extract and correct TMT channel intensities from each light or heavy peptide MS2 spectrum. The TMT data were further processed using pRoloc/Bandle to predict protein subcellular localization via supervised learning. D. Temporal information and spatial information is resolved in MS1 and MS2 levels, respectively. SPLAT allows the subcellular spatial information of the heavy (new) and light (old) subpools of thousands of proteins to be quantified simultaneously in normal and perturbed cells. HL: Half-life.

Figure 2.

Simultaneous measurements of spatial and temporal kinetics under UPR. A. Bar charts showing activation of known ER stress markers upon thapsigargin treatment for 8 hours and 16 hours. X-axis: ER stress markers; y-axis: expression ratio (n=6 normal AC16; n=3 thapsigargin). *: limma adjusted P < 0.01; **: limma adjusted P < 0.001; ***: limma adjusted < 0.0001; error bars: s.d. B. PC1 and PC2 of proteins spatial map showing the localization of confidently allocated proteins in normal and thapsigargin-treated AC16 cells. Each data point represents a protein; color represents classification of subcellular localization. C. Distribution of light (unlabeled) protein features in each of the 12 subcellular compartments (n=3); fill color represents whether the protein is also annotated to the same subcellular compartment in UniProt Gene Ontology Cellular Component terms. D. Histograms of the determined log10 protein turnover rates in control and thapsigargin treated cells (n=3). Text overlay indicates median half-life. E. Boxplot showing the log2 turnover rate ratios in thapsigargin over normal AC16 cells for proteins that are localized to the ER (blue) (T) or not (F); or the Golgi (GA; green). P values: Mann-Whitney test. A Bonferroni corrected threshold of 0.05/13 is considered significant. Center line: median; box limits: interquartile range; whiskers: 1.5x interquartile range. F. Example of best fit curves in the first-order kinetic model at the protein level between normal (gray), and thapsigargin treated (red) AC16 cells showing four known ER stress markers with elevated turnover (HSPA5, RCN3, HSP90B1, and PDIA4). Because the sampling time point is known, the measured relative isotope abundance of a peptide (prior to reaching the asymptote) is sufficient to define the kinetic curve and the parameter of interest (k). G. Turnover rate ratio (thapsigargin vs. normal) of the top proteins with elevated temporal kinetics in UPR within the ER (blue) and Golgi (green); .: Mann-Whitney test FDR adjusted P value < 0.1; *: < 0.05; ** < 0.01; red dashed line: 1:1 ratio; bars: standard error. H. Gene set enrichment analysis (GSEA) of turnover rate ratios in thapsigargin treatment; proteins with faster kinetics are significantly enriched in DNA damage response and UPR pathways. Color: FDR adjusted P values in GSEA; x-axis: GSEA normalized enrichment score (NES). Size: number of proteins in gene set.

Figure 3.

SPLAT captures extensive protein translocation in AC16 cells under UPR. A. Alluvial plot of translocation events (> 0.99 BANDLE translocation probability; estimated FDR < 1%) following thapsigargin treatment showing a cohort of proteins moving from the Golgi apparatus (GA) and lysosome towards the plasma membrane (PM) (n=3). B. Protein spatial map for SLC3A2 (open black circle) in normal (left) and thapsigargin-treated (right) AC16 cells, showing its colocalization with lysosomal proteins in normal cells and in PM proteins in thapsigargin-treated cells. Colors represent allocated subcellular localization. C. Ultracentrifugation fraction profile of SLC3A2 and other amino acid transporters SLC7A5, SLC1A4, SLC1A5 and ion channel proteins SLC30A1, ATP1B1, ATP1B3, and ATP2B1 with similar migration patterns. X-axis: fraction 1 to 10 of ultracentrifugation. Y-axis: relative channel abundance. Bold lines represent the protein in question; light lines represent ultracentrifugation profiles of all proteins classified to a respective localization. Colors correspond to subcellular localization in panel B and for all AC16 data throughout the manuscript; numbers within boxes correspond to BANDLE allocation probability to compartment. D. Immunofluorescence of SLC3A2 (red) against the lysosome marker LAMP2 (green) and DAPI (blue). Numbers in cell boundary: colocalization score per cell. Scale bar: 90 μm. E. Colocalization score (Mander’s correlation coefficient) between SLC3A2 and LAMP2 decreases significantly (two-tailed unpaired t-test P value: 3e–8) following thapsigargin treatment, consistent with movement away from lysosomal fraction (n= 205 normal cells, n = 32 thapsigargin treated cells). Center line: median; box limits: interquartile range; whiskers: 1.5x interquartile range; points: outliers. F. Alluvial plot showing the migration of ER, GA, and nucleus proteins toward the peroxisome/endosome containing fraction in thapsigargin treated cells. G. Ultracentrifugation fraction profile of stress granule proteins UBAP2, USP10, CNOT1, CNOT2, ZC3H7A, and NUFIP2. RNA Granule Score: score from RNA Granule Database (https://rnagranuledb.lunenfeld.ca/). A score of 7 or above is considered a known stress granule protein. Phi: predicted phase separation participation. Circle denotes a prediction of True within the database. RBP: Annotated RNA binding protein on the RNA Granule Database. One circle denotes known RBP in at least one data set; two circles denote known RBP in multiple datasets.

Figure 4.

SPLAT reveals protein-lifetime dependent translocation. A. Histogram showing the similarity in light and heavy proteins in normalized fraction abundance profiles in (left) normal and (right) thapsigargin-treated AC16 cells. X-axis: the spatial distribution distance of two proteins is measured as the average euclidean distance of all TMT channel relative abundance in the ultracentrifugation fraction profiles across 3 replicates; y-axis: count. Blue: distance for 1,614 quantified light-heavy protein pairs (e.g., unlabeled EGFR, heavy SILAC-labeled EGFR). Grey: distribution of each corresponding light protein with another random light protein. P value: Mann-Whitney test. B. Proportion of heavy-light protein pairs with confidently assigned localization that are assigned to the same location (purple) in normal (left; 93%) and thapsigargin-treated (right; 89%) cells. C. Ranked changes in heavy-light pair euclidean distance upon thapsigargin treatment. The difference in heavy-light distances in thapsigargin is adjusted by the average changes in the spatial distance of the light protein with 250 other sampled light proteins to calculate the normalized difference. The majority of proteins show no change (+/− 0.02 in euclidean distance). The positions of EGFR and ITGAV are highlighted. Inset: Z score distribution of all changes. D. Spatial map showing the location of the light and heavy EGFR in normal and thapsigargin-treated AC16 cells. Each data point is a light or heavy protein species. Colors correspond to other AC16 experiments in the manuscript. Numbers correspond to euclidean distance in fraction profiles over 3 replicates. E. Corresponding fraction profiles; x-axis: ultracentrifugation fraction; y-axis: fractional abundance. Post-labeling synthesized EGFR is differentially distributed in thapsigargin and shows ER retention (blue), whereas the preexisting EGFR pool remains to show a likely cell surface localization (pink) after thapsigargin. F-G. As above, for ITGAV. H. Confocal imaging of EGFR immunofluorescence supports a partial relocalization of EGFR from the cell surface toward internal membranes following thapsigargin treatment. Numbers: The mean intensity of the labeled EGFR channel of a 3 pixel border at cell boundaries was divided by mean intensity of the whole cell to estimate the ratio of EGFR at the plasma membrane to the cell interior. Blue: DAPI; Green: EGFR; scale bar: 90 μm. I. Cell areas; Mann-Whitney P: 0.72. Center line: median; box limits: interquartile range; whiskers: 1.5x interquartile range. J. Total EGFR intensity per cell; Mann-Whitney P: 2.2e-05. Center line: median; box limits: interquartile range; whiskers: 1.5x interquartile range. K. Edge/total intensity ratios in normal and thapsigargin-treated AC16 cells (n=71 normal cells; n=93 thapsigargin cells; Mann-Whitney P: 1.7e–4). Center line: median; box limits: interquartile range; whiskers: 1.5x interquartile range.

Figure 5.

Comparison of ER stress induction methods. A. Bar charts showing activation of known ER stress markers upon tunicamycin treatment for 8 hours and 16 hours. X-axis: ER stress markers; y-axis: expression ratio (n=6 normal AC16; n=3 tunicamycin). *: limma adjusted P < 0.01; **: limma adjusted P < 0.001; ***: limma adjusted P < 0.0001; error bars: s.d.. B. Histograms of the determined log10 protein turnover rates in control and tunicamycin treated cells (n=3). C. Boxplot showing the log2 turnover rate ratios in tunicamycin over normal AC16 cells for proteins that are localized to the ER (T) or not (F); Golgi apparatus, or the lysosome. P values: two-tailed t-test. Center line: median; box limits: interquartile range; whiskers: 1.5x interquartile range. A Bonferroni corrected threshold of 0.05/13 is considered significant. D. Gene set enrichment analysis (GSEA) of turnover rate ratios in tunicamycin treatment. Color: FDR adjusted P values in GSEA; x-axis: GSEA normalized enrichment score (NES). Size: number of proteins in the gene set. E. Example of best fit curves in the first-order kinetic model at the protein level between normal (gray), tunicamycin (blue) and thapsigargin (red) treated AC16 cells showing known ER stress markers with elevated turnover in both ER stress inducers (HSPA5, HSP90B1, and PDIA4) as well as stress response proteins with elevated turnover only in tunicamycin (PDIA3, DNAJB11, NIBAN1). F. Alluvial plot showing the migration of ER, GA, and peroxisome/endosome proteins toward the lysosome (left). On the right, the ultracentrifugation fraction profiles of translocating proteins RRBP1, FKBP11, GANAB, MANF, IKBIP, and P3H1 are shown that are targeted toward the lysosome in both tunicamycin and thapsigargin treatment (BANDLE differential localization probability > 0.95). Numbers in boxes are the BANDLE allocation probability in each condition (n=3).

Figure 6.. Applicability in human iPSC-derived cardiomyocytes.

A. Schematic of human iPSC differentiation into cardiomyocytes, carfilzomib treatment, and SPLAT analysis. B. Confocal microscopy images showing sarcomeric disarray in iPSC-CMs upon 48 hrs of 0.5 μM carfilzomib; green: cTNT, red: alpha-actinin; blue: DAPI; scale bar: 20 μm. C–H. Cell viability (%), normalized Seahorse oxygen consumption rate (OCR; pmol/min), basal respiration, maximal respiration, proton leak, and ATP production upon 0 – 5 μM carfilzomib for 24 or 48 hrs; .: adjusted P < 0.1; *: adjusted P < 0.05; **: adjusted P < 0.01, ANOVA with Tukey’s HSD post-hoc at 95% confidence level; n=5. Error bars: s.d. for bar charts in panels C, E, F, G, H; s.em. for the OCR graph in panel D. Colors in panel D: dosage, same as panel C. O: Oligomycin; AA/R: Antimycin A/Rotenone. I. Spatial map with 13 assigned subcellular localizations in iPSC-CMs at the baseline (top) and upon 0.5 μM carfilzomib treatment (n=2). J. Histogram of log10 protein turnover rates (k), with median half-life of 97.4 hours and 100.0 hours in normal and carfilzomib-treated iPSC-CM. K. Proteasome activity in iPSC-CMs treated with 0 (Ctrl) vs. 0.5 μM carfilzomib (Cfzb) for 48 hrs. P value: two-tailed t-test; n = 3. L. Autophagy assay for iPSC-CMs treated with 0 (Ctrl) vs. 0.5 μM carfilzomib (Cfzb) for 48 hrs, and positive control (Pos); data were normalized to DAPI and normal cells. P value: two-tailed t-test; n = 10. M. log2 Turnover rate ratios between carfilzomib-treated and untreated iPSC-CM from the spatiotemporal proteomics data. Proteins assigned the proteasome compartment have significantly increased temporal kinetics; proteins in the lysosome/junction and chromatin/sarcomere compartments have significantly reduced temporal kinetics. P values: Mann-Whitney; with a threshold of 0.05/14 considered significant. Center line: median; box limits: interquartile range; whiskers: 1.5x interquartile range.

Figure 7.

Proteostatic pathways and lesions in carfilzomib mediated cardiotoxicity in iPSC-CMs. A. Changes in protein turnover rates between carfilzomib vs. normal iPSC-CMs across selected cellular compartments; **: P < 0. 01; *: P < 0.05; .: P < 0.1; Mann-Whitney test FDR adjusted P values. error bars: standard error. B. Kinetic curve representations of proteins with accelerated temporal kinetics in carfilzomib including PSMC2 which corresponds to the ratio in panel A, as well as additional ERAD proteins and chaperones; gray: normal iPSC-CM; green: carfilzomib. C. Kinetic curve representations of slowdown of protein kinetics in DSP, DMD, MYH6, and MYH7, corresponding to the ratios in panel A. D–E. Spatial map (PC1 vs. PC2) and ultracentrifugation fraction profiles of D. BAG3 and E. PSME4 in normal and carfilzomib-treated human iPSC-CM, showing a likely differential localisation in conjunction with kinetics changes. White-filled circles: light and heavy BAG3 or PSME4 in each plot. The kinetic curves of BAG3 and PSME4 are in panel B. Numbers at arrows correspond to BANDLE differential localization probability (Diff. Loc. Pr.).