Proteomics and phosphoproteomics datasets of a muscle-specific STIM1 loss-of-function mouse model.

STIM1 is an ER/SR transmembrane protein that interacts with ORAI1 to activate store operated Ca2+ entry (SOCE) upon ER/SR depletion of calcium. Normally highly expressed in skeletal muscle, STIM1 deficiency causes significant changes to mitochondrial ultrastructure that do not occur with loss of ORAI1 or other components of SOCE. The datasets in this article are from large-scale proteomics and phosphoproteomics experiments in an inducible mouse model of skeletal muscle-specific STIM1 knock out (KO). These data reveal statistically significant changes in the relative abundance of specific proteins and sites of protein phosphorylation in STIM1 KO gastrocnemius. Protein samples from five biological replicates of each condition (+/- STIM1) were enzymatically digested, the resulting peptides labeled with tandem mass tag (TMT) reagents, mixed, and fractionated. Phosphopeptides were enriched and a small amount of each input retained for protein abundance analysis. All phosphopeptide and input fractions were analyzed by nano LC-MS/MS on a Q Exactive Plus Orbitrap mass spectrometer, searched with Proteome Discoverer software, and processed with in-house R-scripts for data normalization and statistical analysis. Article published in Molecular Metabolism [1].

Specifications Table

Value of the Data

These proteomics and phosphoproteomics datasets serve as a resource for gaining insights into the signaling perturbations stemming from disrupted Ca2+ homeostasis in skeletal muscle that significantly impact mitochondrial function, proteostasis, and lean mass.

The custom R scripts could be useful for proteomics researchers interested in normalizing PTM and protein measurements from isobaric tag-based quantitative proteomics experiments, using the current data as an example.

These data and methods serve as an example workflow for applying kinase-substrate enrichment analysis (KSEA) to an organism other than human by updating the database of experimentally-determined kinase-substrate pairs to those derived from the species of interest.

Data Description

The deposited raw data was generated from a large-scale proteomics and phosphoproteomics experiment aimed at better understanding the function of STIM1, a key ER/SR protein required for proper Ca2+ homeostasis and mitochondrial function in skeletal muscle [1]. Fig. 1 illustrates the experimental workflow. A model of inducible STIM1 deficiency (iSTIM−/−) was generated to circumvent developmental defects. Briefly, STIM1fl/fl mice were subjected to three consecutive treatments with tamoxifen (TMX), which activated the transcription of Cre recombinase, resulting in skeletal muscle-specific loss of STIM1 protein. After 8 weeks, mice were subjected to a 16 hour fast before gastrocnemius samples were collected and processed as described below in Experimental Design. Raw data were searched using Proteome Discoverer (PD), filtering peptide and protein identifications to a confidence threshold of 1% FDR, and calculating quantitative values for each TMT channel—keeping the quantification for the input and phosphopeptide fractions separate using the Study Factors feature in PD. The Proteins and Peptide Isoforms output were exported from PD as .txt files and were further processed with custom R scripts, which have been made available as part of this Data in Brief (doi: 10.5281/zenodo.6226588) along with our analyzed data [dataset] [2]. Fig. 2 demonstrates how the R script normalized the TMT11plex quantitative proteomics and phosphoproteomics data for any subtle (and unintended) differences in total sample “loading” (i.e., µg of peptide labeled) on each TMT channel. First, to calculate sample-specific loading normalization factors for each TMT channel, the sum of the quantitative values for all peptides in the input fraction was divided by the mean of the TMT channel sums across all 11 samples. Next, all protein and phosphopeptide quantitative values were divided by the sample-specific loading normalization factors. Loading-normalized quantitative values were log2 transformed to approximate a normal distribution. Relative phosphopeptide occupancy (i.e., phosphorylation stoichiometry fold change) was obtained by taking the log2 abundance of each phosphopeptide and subtracting (division in the log space) the log2 abundance of the respective protein (if measured). Lastly, the protein abundance, phosphopeptide abundance, and phosphopeptide relative occupancy values were assessed for significance using the limma R package, correcting the p.values for multiple hypothesis testing using the Benjamini-Hochberg approach. In large-scale datasets, a phosphopeptide with striking changes in abundance could result from a large change in the abundance of a highly regulated protein with a fixed PTM stoichiometry; however, phosphopeptides exhibiting significant abundance changes with loss of STIM1 generally do not appear to be such cases. Fig. 3 illustrates the log2 fold change values from the phosphopeptide abundance on the x-axis and the log2 fold change values for the relative occupancy (phosphopeptide change normalized to the protein's change) on the y-axis. The strong correlation of phosphopeptide abundance and occupancy for significantly changing phosphopeptides is suggestive of altered protein kinase/phosphatase activity in this model. These data are associated with a research article published in Molecular Metabolism [1], wherein kinase-substrate enrichment analysis (KSEA) [3] was used to predict changes in the activities of protein kinases resulting from the loss of STIM1, a subset of which were validated with biochemical methods.

Fig. 1

Proteomics workflow and key data metrics. Mice with skeletal muscle-specific loss of STIM1 (iStim−/−, KO) and control animals (STIM1fl/fl) were sacrificed after an overnight fast and gastrocnemius muscles were collected. Proteins were extracted from five biological replicate samples, enzymatically digested, the resulting peptides labeled with unique tandem mass tag (TMT) reagents. Samples were mixed and fractionated by high-pH reversed phase (HPRP) chromatography. Phosphopeptides were enriched by immobilized metal affinity chromatography (IMAC) after a small amount of the input material retained from each fraction for assessment of protein abundance. All phosphopeptide and input fractions were analyzed by nano LC-MS/MS on a Q Exactive Plus (QE+) Orbitrap mass spectrometer, searched with Proteome Discoverer software, and processed with in-house R-scripts for normalization of quantitative values and analysis of statistical significance. For each measurement type, the numbers quantified and statistically significant are indicated. Article published in Molecular Metabolism[1].

Fig. 2

Sample loading normalization with custom R scripts. Log2 intensity of each TMT channel for all protein relative abundance measurements before (A-B) and after (C-D) normalization for any subtle (and unintended) deviations in equal peptide input. The data are displayed as both density plots (A, C) and box-and-whisker plots (B, D). Phosphopeptide quantitative values (not shown) were corrected with the same loading normalization factors.

Fig. 3

Correlation between phosphopeptide abundance and relative occupancy. Scatter plot of phosphopeptide relative abundance on the vertical axis and the phosphopeptide relative occupancy on the y axis for the same protein phosphorylation sites (both measurements expressed as log2 KO/WT). Red dots represent phosphosites that only exhibited a statistically significant (10% FDR) change when considering phosphopeptide abundance, whereas blue indicates sites only changing by phosphopeptide relative occupancy, and site changing by both metrics are shown in purple. Sites that were not significant (N.S) are shown in grey (also slightly transparent to better visualize statically significant sites behind them). Phosphopeptides for which a corresponding protein was not measured are not shown.

Experimental Design, Materials and Methods

Experimental design

Tissue lysis and protein digestion

Gastrocnemius tissue was prepared for proteomics following a published protocol [4] (Fig. 1). After pulverizing tissue into a powder using a mortar and pestle under liquid nitrogen, ∼15 mg of each sample was solubilized in ice-cold urea Lysis buffer (8 M urea in 50 mM Tris, pH 8.0, 40 mM NaCl, 2 mM MgCl2, 1x Roche cOmplete ULTRA EDTA-free protease inhibitor mini tablet, 1x PhosStop phosphatase inhibitor tablet, 5 mM NaF and 10 mM Na pyrophosphate). Samples were lysed with a TissueLyzer (Qiagen) for 30 sec at 30 Hz, subjected to three freeze-thaw cycles, and further disrupted with a probe sonicator (three 5-sec bursts, power setting of 3). Lysates were centrifuged at 10,000 × g for 10 min at 4 °C and supernatants retained. Protein concentration was determined by BCA, and a 400 µg aliquot of each sample was adjusted to 2.5 mg/mL with Urea Lysis Buffer. To reduce and alkylate cysteine residues, samples were incubated with 5 mM DTT at 32 °C for 30 min, cooled to RT, and incubated with 15 mM iodoacetamide for 30 min in the dark at RT. Following quenching of unreacted iodoacetamide by the addition of equimolar DTT (up to 15 mM), samples were enzymatically digested with 5 µg LysC (100:1 w/w, protein to enzyme) at 32 °C for 4 hr. After dilution to 1.5 M urea with 50 mM Tris (pH 8.0), 5 mM CaCl2, samples were digested with trypsin (50:1 w/w, protein:enzyme) overnight at 32 °C. Samples were acidified to 0.5% TFA and centrifuged at 4000 × g for 10 min at RT to pellet insoluble material. Supernatant containing soluble peptides was subjected to solid phase extraction (SPE) with a Waters 50 mg tC18 SEP-PAK column to remove salt, eluting once with 500 µL 25% acetonitrile/0.1% TFA and twice with 500 µL 50% acetonitrile/0.1% TFA. The 1.5 mL eluate was frozen and dried by speed vac.

Peptide isobaric labeling and fractionation

Samples were resuspended in 110 µL of 200 mM triethylammonium bicarbonate (TEAB), vortexed, and 10 µL of each removed for making a “pool” mixture representing all 10 samples. Each sample and the pool were mixed with a different 11-plex Tandem Mass Tag (TMT) reagent (0.8 mg resuspended in 50 µL 100% acetonitrile) and shaken for 4 hr at RT. As a QC measure, 1 µL of each was removed (freezing the remainder at -80 °C), mixed, and subjected to the nLC-MS/MS workflow described below but searching with TMT as a variable modification on the peptide N-terminus. After calculating 91% TMT labeling efficiency, the remainder of each sample was thawed at RT and 0.8 µL 50% hydroxylamine was added to quench unreacted TMT reagent, shaking for 15 min at RT. All samples were mixed, frozen, and dried by speed vac. The mixture was resuspended in 1 mL 0.5% TFA and subjected to SPE again with a Waters 100 mg tC18 SEP-PAK SPE column as described above and the eluate was frozen and dried by speed vac. The 4.4 mg TMT-labeled peptide mixture was fractionated by high pH reversed-phase (HPRP) chromatography as described previously [5]. After each of 12 fractions was lyophilized and resuspended in 1 mL 80% acetonitrile/0.15% TFA, 50 µL was removed for analysis of unmodified peptides, frozen, and dried by speed vac. The remaining 950 µL aliquot was subjected to immobilized metal affinity chromatography (IMAC) to enrich phosphopeptides based on a previous protocol (Grimsrud et al., 2012). Ni-NTA Magnetic Agarose Beads were washed 3X with water, incubated in 40 mM EDTA, pH 8.0 for 30 min while shaking, washed 3X with water, incubated with 100 mM FeCl3 for 30 min while shaking, and washed 4X with 80% acetonitrile/0.15% TFA. Samples were resuspended in 1 mL 80% acetonitrile/0.15% TFA, added to the beads, and incubated for 30 min at RT while shaking. Flow-through was retained and subjected to a second round of IMAC, in both cases washing the beads 3X with 80% acetonitrile/0.15% TFA and eluting for 1 minute by vortexing in 100 µL of 50% acetonitrile, 0.7% NH4OH. Eluted phosphopeptides were acidified immediately with 50 µL 4% formic acid, frozen, and dried by speed vac.

Analysis of phosphopeptide and input fractions by nLC-ESI-MS/MS

Samples were analyzed by nanoflow liquid chromatography (nLC) on an EASY-nLC 1200 UHPLC, followed by electrospray ionization (ESI) with an EASY-Spray source, and tandem mass spectrometry (MS/MS) with a Q Exactive Plus Hybrid Quadrupole-Orbitrap (all instrumentation from Thermo Fisher Scientific). Prior to injection, sample fractions were resuspended in 22 µL of 0.1% FA. Phosphopeptide fractions were analyzed with technical triplicate runs, injecting 6.5 µL (phosphopeptide). Input fractions were analyzed with at least duplicate runs, injecting 1 µL (approximately 1 µg) for the first and adjusting the volume for the second to target an MS1 base peak intensity of 1 × 1010. For each injection, the sample was first trapped on an Acclaim PepMap 100 C18 trapping column (3 µg particle size, 75 µm × 20 mm) with solvent A (0.1% FA) at a variable flow rate dictated by a maximum pressure of 500 bar, after which the analytical separation was performed over a 105 min gradient (flow rate of 300 nL/min) of 5% to 40% solvent B (90% ACN, 0.1% FA) using an Acclaim Pep-Map RSLC C18 analytical column (2 µg particle size, 75µg × 500 mm column (Thermo Fisher Scientific) with a column temperature of 55 °C. MS1 (precursor ions) was performed at 70,000 resolution with an AGC target of 3 × 106 ions and a maximum injection time (IT) of 60 ms. MS2 spectra (product ions) were collected by data-dependent acquisition of the top 10 most abundant precursor ions with a charge greater than 1 per MS1 scan, with dynamic exclusion enabled for a window of 30 s. Precursor ions were filtered with a 0.7 m/z isolation window and fragmented with a normalized collision energy of 30. MS2 scans were performed at 35,000 resolution, with an AGC target of 1 × 105 ions and a maximum IT of 60 ms.

Data analysis

Data availability

All Thermo .raw files were archived on Proteome Xchange [6] (accession: PXD030674) using the jPOST Repository (accession: JPST001437) [7]. The entire analyzed dataset (doi: 10.17632/495z7s99zs.2), and the custom R scripts (doi: 10.5281/zenodo.6226588) used to produce it are available online. See below for detailed descriptions of data acquisition and processing, as well as the methods behind additional bioinformatics analyses highlighted in Wilson et. al. [1] Figure 6 (KSEA) and Tables S1-2 (top regulated phosphosites/proteins).

Raw data analysis with proteome discoverer

Thermo .raw files were processed with Proteome Discoverer 2.5 (PD2.5) from Thermo Fisher Scientific. A complete mouse proteome .fasta file containing 55,485 reviewed (Swiss-Prot) and unreviewed (TrEMBL) protein sequences, downloaded from UniProt on 2/3/2021, was imported into PD2.5. Study Factors were designated as “Fraction” (Input or Phosho) and “Genotype” (KO, WT, or pool). Separate Processing Workflows were created for searching the Input and Phospho .raw files, as PTM-containing peptides are scored differently by search algorithms than unmodified peptides and should ideally be considered independently for FDR estimation when working with PTM-enriched and unenriched samples [8]. MS1 precursor mass recalibration was performed after a pre-search using the Spectrum Files RC node. Subsequent searches with both Sequest HT and MS Amanda 2.0 were conducted at 10 ppm (recalibrated) and 0.02 Da mass tolerances for precursor and fragment ions, respectively. Modifications considered included oxidation on methionine (variable +15.995 Da), carbamidomethylation of cysteine (fixed +57.021 Da), and TMT6plex addition to lysine peptide N-termini (fixed +229.163 Da). Note, it is typical to use the default “TMT6plex” for searching, as the nominal mass is the same as for TMT11plex. However, the quantitative method was specified as TMT11plex, including isotopic impurity corrections for each channel. Searches allowed for up to 2 missed K/R cleavages by trypsin (full specificity). All .raw files from phosphopeptide analyses also considered modification of serine, threonine, and tyrosine residues with a phosphoryl group (variable + 79.966 Da). Percolator [9] was used to filter peptide spectral matches (PSMs) from each search algorithm to an estimated FDR of less than 1%, and ptmRS [10] was used in “PhosphoRS Mode” to determine the probabilities of localizing phosphorylation to specific S/T/Y residues. As part of a common Consensus Workflow, the Peptide Isoform Grouper node was used to collapse PSMs down to unique peptides, while maintaining FDR < 1% at the peptide level and requiring a site localization probability of 90% or greater for differentiating phosphopeptide positional isoforms. The Protein FDR Validator node was used to group peptides into proteins using the rules of strict parsimony and maintain FDR < 1% at the protein level. The Reporter Ion Quantifier node was used to filter PSMs for quantitative quality, setting an upper threshold of 0.5 precursor co-isolation and an average reporter ion S/N > 2.5. For PSMs meeting these criteria, quantitative values were summed together at the peptide and protein levels. Peptides shared between protein groups and peptides with oxidized methionine were excluded from protein quantitation calculations. The quantitative values for each data type (Phospho, Input) were kept separate based on the “Fraction” Study Factor designation noted above. As another key to our workflow for PTM quantitation, none of the normalization options in PD2.5 were selected, as this is performed by our custom R script as described below.

Statistical analysis using a custom R script

Upon completion of the search in PD2.5, The “Protein” and “Peptide Isoforms” tabs were exported from the .pdResult file as tab-delimited .txt files. Subsequent analysis used custom R scripts (doi: 10.5281/zenodo.6226588), modified from previous code [11] for use in R (version 4.0.5) in conjunction with Rstudio [12]. First, the phosphopeptide and protein abundance data were normalized to correct for any subtle (and unintended) deviations in equal amounts of peptide on each TMT channel in the input fraction, as described previously [13], [14]. To calculate channel-specific loading normalization factors, reporter intensities for each peptide in the input fraction were summed together for each TMT channel and each channel's sum was divided by the average of all channels’ sums. Reporter intensities for peptide isoforms from the phosphopeptide fraction and proteins from the input fraction were divided by the loading control normalization factors for each TMT channel, respectively (Fig. 2). All loading control-normalized quantitation values were converted to log2 space. The limma (3.46.0) R package [15] was applied to create linear models for differential relative abundance of proteins and phosphopeptides, with moderated statistics computed from the eBayes function. Benjamini-Hochberg correction for testing of multiple hypotheses was applied, considering Padjusted ≤ 0.1 (10% FDR) as statistically significant. Protein level quantification was performed exclusively on Master Proteins—UniProt entries representing a group of parsimonious proteins containing common peptides identified at 1% qualitative FDR. Relative changes in phosphopeptide abundance were computed before (abundance) and after (relative occupancy) normalization for changes in the respective “Master Protein”. Relative Occupancy normalization was achieved by subtracting the Master Protein log2 fold changes (relative to the sample mean) from phosphopeptide abundance log2 fold changes (relative to the sample mean) for each sample (Fig. 3). The ggplot2 R package was used for making all graphs.

Kinase-substrate enrichment analysis (KSEA)

Using the previously published KSEAapp R package [3], the function KSEA.Complete, was altered to not perform log2 transformation on fold changes from the phosphopeptide data, as the KO vs. WT values were already log2 transformed from the analysis described above (section 2.2.3.). The kinase-substrate database used to predict kinase activities was updated from phosphositeplus [16] for compatibility with mouse samples [17]. First, the file “Kinase_Substrate_Dataset.gz”, which contains “experimentally determined substrates, sequences, cognate kinases, and metadata curated from the literature” was downloaded (phosphosite.org/staticDownloads). Kinase-substrate pairs were subsequently filtered to include only those with mouse substrates, regardless of the species of a given kinase. For example, some results populating the phosphositeplus [16] database may have used a purified human kinase in a reaction with mouse tissue homogenate, but it is critical to only consider sites of phosphorylation on mouse proteins, as the surrounding sequences are often different between species. The default horizontal barplot that displays the kinase with its respected z.score was recreated using ggplot2, but instead of having to alter the substrate cutoff argument in the function, the plot was designed to display the number of identified kinase peptide substrates next to the respective bar in the barplot. See Wilson et al. [1]. for KSEA results and kinase activity validation. Note, the altered KSEA code is not included in the present communication, but can be easily reproduced by applying the changes described above to the original R package [3].

Ethics Statement

The mouse experiments followed a protocol (A064-16-03) approved by the Duke University School of Medicine's Institutional Animal Care and Use Committee and NIH guidelines. Animals were grouped housed in a temperature (22 °C) controlled room, with successive 12 hour periods of light and darkness and continual access to water and food. Experimental groups consisting of male animals were chosen randomly.

CRediT Author Statement

Scott Lyons: Methodology, Software, Investigation, Data curation, Visualization, Writing – original draft; Rebecca Wilson: Investigation, Writing – review & editing; Deborah Muoio: Funding acquisition, Project Administration, Writing – review & editing; Paul Grimsrud: Conceptualization, Methodology, Investigation, Funding acquisition, Writing – original draft.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Subject	Omics: Proteomics
Specific subject area	Proteomics and Phosphoproteomics of skeletal muscle
Type of data	Raw Mass Spectrometry DataAnalyzed Proteomics DataCustom R ScriptsFigures
How the data were acquired	nLC-MS/MS analysis on Q Exactive Plus Hybrid Quadrupole-Orbitrap mass spectrometer (Thermo Scientific)
Data format	Raw DataAnalyzed DataData Analysis Code
Description of data collection	STIM1fl/fl mice were treated with tamoxifen (TMX), resulting in STIM1 knockout (KO) in skeletal muscle (STIM1−/−). After 8 weeks, STIM1−/− animals and STIM1fl/fl controls were subjected to a 16 h fast. Gastrocnemius muscle samples were collected and proteins were isolated and digested with trypsin. Peptides from five animals per condition (STIM1−/− vs. STIM1fl/fl), along with a pooled mixture of each sample, were labeled with TMT11plex reagents. Samples were mixed, fractionated by high pH reversed-phase (HPRP) chromatography, and subjected to phosphopeptide enrichment using immobilized metal affinity chromatography (IMAC) after retaining a small portion of the input material from each fraction. All fractions were resuspended in 0.1% formic acid, injected onto a precolumn using a EASY-nLC 1200 UHPLC, separated over a 120-minute gradient of increasing concentrations of acetonitrile on an analytical column, and subjected to data-dependent acquisition (DDA) on a Q Exactive Plus Orbitrap. The raw MS/MS data were analyzed with Proteome Discover 2.5 (PD), searching against the mouse UniprotKB database and calculating quantitative values for peptides and proteins. Output of PD was then further processed using a custom R script for normalizing protein and post-translational modification (PTM) quantitation for any unintended differences in sample input across TMT channels. Relative changes in phosphopeptide abundance between genotypes were further normalized by any changes in the corresponding protein to calculate phosphopeptide relative occupancy changes with loss of STIM1.
Data source location	Duke Molecular Physiology InstituteDuke University School of Medicine300 North Duke St.Durham, NC 27701
Data accessibility	Raw DataRepository Name: Proteome XchangeData identification number: PXD030674Direct link to the dataset: http://proteomecentral.proteomexchange.org/cgi/GetDataset?ID=PXD030674Repository Name: jPOSTData identification number: JPST001437Direct link to the dataset: https://repository.jpostdb.org/entry/JPST001437Analyzed DataRepository Name: Mendeley DataData identification number: doi: 10.17632/495z7s99zs.2Direct link to the dataset: https://data.mendeley.com/datasets/495z7s99zs/2
	Data Analysis Code:Repository Name: GitHubData identification number: doi: 10.5281/zenodo.6226588Direct link to the dataset: https://github.com/dmpio/stim1-data-in-brief
Related research article	R. J. Wilson, S. P. Lyons, T. R. Koves, V. G. Bryson, H. Zhang, T. Li, S. B. Crown, J. D. Ding, P. A. Grimsrud, P. B. Rosenberg, D. M. Muoio, Disruption of STIM1-mediated Ca(2+) sensing and energy metabolism in adult skeletal muscle compromises exercise tolerance, proteostasis and lean mass, Mol Metab (2021) 101429. https://doi.org/10.1016/j.molmet.2021.101429.