Natural aging and ovariectomy induces parallel phosphoproteomic alterations in skeletal muscle of female mice

Mouse and Tibialis Anterior (TA) muscle morphometrics and torque analysis

Significant differences in body mass across the groups were measured (p < 0.001) (Figure 1B). Mean body masses for Sham, Ovx, YA, and OA mice were 35.8 ± 3.9, 38.4 ± 6.2, 22.3 ± 1.0, and 30.2 ± 5.4, respectively. Uterine mass was significantly different between Ovx and Sham mice (p = 0.016), with the mean uterine masses being 19.3 ± 5.1 and 135.7 ± 70.3, respectively (Figure 1C). In addition, cytology confirmed persistent diestrus in Ovx mice and estrous cycling in Sham mice, as expected. Sham and Ovx TA muscle mass were significantly different from OA (p = 0.001) (Supplementary Figure 1A), with mean muscle masses of 49.6 ± 4.6, 53.7 ± 3.1, 45 ± 1.2, and 37.1 ± 4.7 mg for Sham, Ovx, YA, and OA, respectively. TA muscle mass normalized to body mass was greater in YA mice than Sham, Ovx, and OA (p = 0.003) (Supplementary Figure 1B).

No significant differences among the four groups in absolute torque for pre-tetanic twitch torque, maximal isometric tetanic torque, or post-tetanic twitch torque were measured (0.61 ± 0.09, 2.54 ± 0.53, and 0.84 ± 0.12 mN∙m, respectively for all mice; p ≥ 0.158) (Supplementary Figure 1C–1E). However, when body mass was considered, significant differences in normalized torque (torque/body mass) were measured in pre-tetanic twitch torque (p = 0.034), tetanic torque (p = 0.044), and post-tetanic twitch torque (p = 0.003) (Supplementary Figure 1F–1H). These data demonstrate robust muscle contraction intended to trigger protein phosphorylation within each group, as also shown by the twitch torques after the tetanic contraction being greater than the pre-tetanic twitch torques. To note, this study was not designed or powered to detect contractile differences between estrogen-deficient and -replete groups as has been shown previously [15, 29, 35, 36].

Comparative skeletal muscle phosphoproteomes in Ovx/Sham and OA/YA mice

In the two models of estrogen-deficient mice, a total of 2,593 phosphopeptides and 3,507 phosphopeptides were identified in Ovx/Sham and OA/YA TA muscles, respectively (Figure 2A and 2B). Phosphorylation was most prevalent on serine (Ser, ~ 78%), followed by threonine (Thr, ~ 17%), and tyrosine (Tyr, ~ 4%), which is consistent with literature and was similar across both datasets (Figure 2C and 2D) [37], demonstrating reproducible sample preparation and peptide detection. Further analysis identified 222 and 408 significant and differentially regulated phosphopeptides in Ovx/Sham and OA/YA datasets, respectively (Figure 2E and 2F).

Figure 2. Characteristics of the Ovx/Sham and OA/YA phosphoproteomes. Proteome Discoverer (v2.4) was used for database search and identification of phosphopeptides and proteins. Venn diagram of identified phosphopeptides unique to each group and common to both groups in (A) Ovx/Sham and (B) OA/YA datasets. Prevalence of phosphorylation on amino acid residues serine (S), threonine (T), and tyrosine (Y) in (C) Ovx/Sham and (D) OA/YA datasets. Volcano plots of differentially regulated phosphopeptides (p < 0.05 and 1.4-fold change) in (E) Ovx/Sham and (F) OA/YA mice.

The datasets were processed to identify intersecting proteins and phosphosites, which were then deemed to be associated with estrogen deficiency. After filtering for robustness (i.e., relative abundance detected in at least 3 out of the 4 biological replicates in Ovx or Sham and 2 out of 3 biological replicates in OA or YA groups), 66 estrogen deficiency-associated proteins were identified (Supplementary Table 1). Of these 66 estrogen deficiency-associated proteins, four were significant and differentially phosphorylated (eukaryotic translation initiation factor 4E-binding protein 1, ankyrin repeat domain-containing 2, heat shock protein beta-6, and heterogeneous nuclear ribonucleoprotein U-like protein 2) in the OA/YA Adult dataset. Two (phosphoacetylglucosamine mutase and myc box-dependent interacting protein1) were differentially phosphorylated in the Ovx/Sham dataset.

A total of 21 estrogen deficiency-associated phosphosites were found in the two datasets, of which 4 phosphosites were differentially regulated in both, tumor protein D54 Ser-166, ATP synthase subunit alpha Ser-521, calpastatin (CAST) Ser-82, and H/ACA ribonucleoprotein complex subunit DKC1 (DKC1) Ser-481 (Table 1). CAST Ser-82 and DKC1 Ser-481 were the only two phosphosites that had the same directionally across both datasets, upregulation and downregulation, respectively in estrogen-deficient mice. Although we did identify estrogen deficiency-associated phosphosites on six sarcomeric proteins (myozenin-1 Ser-134, desmin Ser-28, nexilin Ser-559, junctophilin-1 Ser-501, troponin T Ser-2, and myosin regulatory light chain 2 Ser-16), only myozenin-1 Ser-134 was differentially regulated in the OA/YA dataset; whereas, myozenin-1 Ser-164 and troponin T Ser-2 were differentially regulated in the Ovx/Sham datasets. These data show that Ovx and OA mice compared to their respective control, Sham or YA, had some similar phosphoproteomic alterations in contracted muscles induced by estrogen deficiency.

GO term overrepresentation analysis

To gain insight on the biological significance of estrogen deficiency in the skeletal muscle phosphoproteome after contraction in Ovx and OA mice, a comparative GO term overrepresentation analysis was performed. The top 10 most overrepresented GO terms for cellular compartment, molecular function, and biological process are shown in Figure 3A–3C, respectively. There was significant overrepresentation of GO terms (FDR < 0.05) across all GO domains with numerous similarities between both datasets. GO molecular function analysis had similar terms enriched in both datasets relating to actin/actin filament binding, structural constituent of the cytoskeleton and muscle, calmodulin binding, tropomyosin binding, and translation regulator activity. However, there was a preferential enrichment in phosphatase activity in the Ovx/Sham dataset; whereas, a preferential enrichment in translation activity in the OA/YA dataset (Figure 3A). All top 10 GO cellular component terms pertaining to the muscle fiber were enriched across both datasets, with reduced enrichment in the myosin complex in the OA/YA dataset. Likewise, all GO biological process terms except one, glycogen metabolic process, was enriched across both datasets. Overall, the top 10 GO term enrichments across many molecular functions, and essentially all cellular components and biological processes were consistent in the Ovx/Sham and OA/YA datasets suggesting that estrogen deficiency confer parallel changes in the phosphoproteome of contracted muscles.

Figure 3. Comparative GO term enrichment analysis between Ovx/Sham and OA/YA. Phosphopeptides were mapped back to their precursor protein and submitted for Gene Ontology (GO) overrepresentation analysis using the clusterProfiler package in R. The top 10 overrepresented GO terms in the dataset for (A) molecular functions, (B) cellular components, and (C) biological processes are listed. P-value was adjusted using Benjamini – Hochberg post-hoc analysis for multiple comparison. Significant GO terms were accepted at p.adjusted < 0.05.

KEGG and Reactome pathway overrepresentation analysis

To further explore biological pathways associated with estrogen deficiency in the datasets, overrepresentation analysis of KEGG and Reactome pathways was performed. Using the K-means clustering algorithm, mutual overlapping pathways were clustered together via the similarity between nodes, and the stronger the similarity, the shorter and thicker the connecting lines (Figure 4). Similar to the GO term analyses, among the top 10 pathways, all but one KEGG and three Reactome pathways were associated across both datasets. Only the “pathways of neurodegeneration – multiple diseases” from the KEGG analysis was unique to the OA/YA dataset. Out of the four apoptosis related pathways from the Reactome analysis, three were unique to the Ovx/Sham dataset suggesting that apoptosis was more of a factor for Ovx than it was for natural OA mice. In summary, many enriched KEGG and Reactome pathways were similar across both datasets suggesting that these pathways may be sensitive to estrogen levels and may contribute to muscle contractile dysfunction.

Figure 4. Comparative KEGG and Reactome pathway enrichment analysis between Ovx/Sham and OA/YA. Significant and differentially regulated phosphopeptides (adjusted p-value < 0.05 and |FC| ≥ 1.4) were mapped back to their precursor protein and the list was submitted to R for analysis of KEGG and Reactome pathways. The top 10 overrepresented (A) KEGG and (B) Reactome pathways were clustered using K-means clustering. The size of the circle represents the number of proteins associated to the pathway, and the connecting lines represents the strength of the similarity (i.e., shorter and thicker lines correspond to stronger similarity and inversely.

IPA’s predicted downstream effects and upstream regulator analytics

A comparative analysis between the two core analyses was then performed to identify parallel activation states via IPA’s Z-score algorithm across IPA-derived molecular and cellular functions, canonical pathways, and upstream regulators. For reference, IPA's Z-score indicates a predicted overall activation or inhibition of functions/pathways and activated or inhibited upstream regulators, where a negative Z-score signifies inhibition/inhibited and a positive Z-score signifies activation/activated [38].

IPA molecular and cellular functions analysis

Seventeen biological processes were identified across both datasets with similar directionality in activation and were broadly related to muscle quantity and function, filament formation and stabilization, transport of molecules, autophagy, and RNA expression and translation (Figure 5A). Four molecular and cellular functions (quantity of muscle, necrosis, force generation, and expression of RNA) were in activation states, with necrosis activation being significant (Z-score = 2.52) in the OA/YA dataset. The remaining molecular and cellular functions were in inhibition states, with a strong inhibition of polymerization of filaments (Z-score = 1.94) in the Ovx/Sham dataset. These parallel molecular and cellular functions identified across both datasets suggest a compromise not only in the force-generating capacity of skeletal muscle but also mechanisms involved in skeletal muscle maintenance and integrity when estrogen is deficient.

Figure 5. IPA’s predictive downstream effect and upstream regulator analytics between Ovx/Sham and OA/YA. Phosphopeptides identified in both the Ovx/Sham and OA/YA datasets were submitted to IPA for comparative analysis. (A) canonical pathways, (B) functions, and (C) kinases and phosphatases using IPA’s predictive activation Z-score to determine downstream and upstream effects of estrogen deficiency. Pathways, functions, and kinases identified in both datasets are represented. Significant Z-scores were accepted at |Z| ≥ 2. Red box denotes significant Z-scores.

Animals

Female C57BL/6J (6, 4, and 24 mo) mice were purchased from Jackson Laboratories (Bar Harbor, ME, USA). All mice were housed in groups of four to five with access to phytoestrogen-free food and water ad libitum. The room was maintained on a 14:10 h light/dark cycle.

Experimental design

The study design to investigate the skeletal muscle phosphoproteome after contraction in two model of estrogen deficiency, Ovx vs. Sham (6 mo; n = 4/group) and OA (24 mo) vs. YA (4 mo; n = 3/group), is summarized in Figure 1A. C57BL/6J female mice were randomly assigned to a sham or ovariectomy surgery. Sham and Ovx mice were 14 mo at the terminal experiment and Ovx mice were estrogen deficient for 32 wk. Vaginal cytology was performed 4 weeks post-surgery and uterine mass was measured at the terminal experiment. OA and YA mice were purchased from Jackson Lab and were acclimatized for two weeks before the terminal experiment. We chose 24 mo of age for the OA group, as C57BL/6J mice are ovarian-senescent as previously determined by cessation of estrous cycling between 13-16 mo [57] and in another study at <20 mo as determined by vaginal cytology and plasma E2 [58]. Age, evaluation of vaginal cytology, and measurement of uterine mass are proxies to indicate estrogen deficiency; a shortcoming of this study is that serum E2 was not measured.

An in vivo muscle contraction protocol of the anterior crural muscles was performed on each mouse in all four groups. The protocol (detailed below) primarily consisted of one maximal isometric tetanic contraction. A single tetanic contraction was chosen to mitigate potential complications of muscle fatigue that might occur with multiple maximal contractions. Immediately following the muscle contraction protocol, the contracted tibialis anterior (TA) muscles from anesthetized mice were dissected within 5 min, flash frozen in liquid nitrogen, and stored at −80ºC until sample preparation and phosphopeptide enrichment for nanoflow LC-MS/MS acquisition. The TA muscle was selected for study for several reasons including; (1) consistency with phosphoproteomic results from our previous work on non-contracted muscles [32], (2) specificity of stimulating contraction and measuring peak tetanic and twitch torques of the anterior crural muscles [59], and (3) the muscle is composed of type II fibers [60] rendering it more sensitive to contraction-induced contractile protein phosphorylation relative to a muscle composed of type 1 fibers [61].

Sham and ovariectomy surgeries

Anesthetized (1.75% isoflurane and 200 ml O₂ per min) mice received a subcutaneous injection of slow-release buprenorphine (1 mg/kg) immediately prior to surgeries. Bilateral ovariectomy via abdominal incisions were made to locate the ovaries which were excised in Ovx mice, and located but not removed in Sham mice. The abdominal muscle wall incisions were closed with 6–0 silk sutures and 7 mm wound clips closed the skin incision.

In vivo muscle contraction protocol

Mice were anesthetized (1.25% isoflurane and 125 O₂ per min) and positioned on a temperature-controlled platform as previously described [62]. Briefly, the left knee was immobilized and the left foot secured to an aluminum “shoe” that was attached to the shaft of an Aurora Scientific 300B servomotor. Sterilized platinum needle electrodes (FE212, Grass Technologies, Warwick, RI, USA) were inserted to stimulate the left common peroneal nerve, and voltage and electrode placement were optimized with 3–5 twitch contractions (0.1 ms pulse). Following optimization, an isometric tetanic contraction of the anterior crural muscles (the TA being the primary muscle in that group) was elicited (1000 Hz for 1000 ms with 0.1 ms pulses), with a twitch contraction performed 10 s prior and two post-tetanic twitch contractions at 2 s and 30 s after the tetanic contraction [36, 59].

Protein extraction, digestion, and phosphopeptide enrichment

Frozen dissected TA muscles were prepared as previously described [32]. Briefly, TA muscles were pulverized into powder with a cryo-grinder (liquid nitrogen cooled mortar and pestle), lysed (10 μl lysis buffer per mg of tissue) in protein lysis buffer (7 M urea, 2 M thiourea, 0.4 M Tris pH7.5, 20% acetonitrile, 4 mM TCEP) with 1X HALT Protease and Phosphatase Inhibitor Cocktail (Thermo Fisher Scientific, Rockford, IL, USA), and sonicated for 5 s using a probe sonicator (Branson Digital Sonifier, Emerson, St. Louis, MO, USA) set at 30% amplitude, all done on ice. After sonication, a 160 μl aliquot of each lysate was placed in the Barocycler^® NEP2320 (Pressure Biosciences, South Easton, MA, USA) at 37ºC, with pressure cycles set at 35,000 psi for 20 s, then 0 psi for 10 s for 60 cycles for further protein homogenization. Once pressure cycling was complete, samples were transferred to a new 1.5 ml Eppendorf protein LoBind tube and 200 mM chloroacetamide stock solution was added for a final concentration of 8 mM to alkylate proteins (1:24 dilution) and incubated for 15 min at room temperature. Samples were spun down at 15,000 × g for 10 min at 18°C. Aliquots of supernatant were used to determine protein concentration using the Bradford assay. For trypsin digestion, 500 μg total protein was digested with 12.5 ug sequencing grade modified trypsin (Promega, Madison, WI, USA) and incubated at 37°C in a warm air incubator overnight (~ 16 h). Samples were acidified to 0.2% trifluoroacetic acid (TFA) to a final volume of 1 ml and extracted using 1 cc Oasis HLB Solid Phase Extraction cartridges (Waters Corporation, Milford, MA, USA) for clean-up. Briefly, 1cc HLB cartridges were equilibrated by passing 1 mL of 80% acetonitrile, 0.1% TFA over the cartridge followed by passing 1 mL of 0.1% TFA. The sample was passed over the cartridge followed by washing the sample with 1 mL 1% acetonitrile, 0.1% TFA over the cartridge. Peptides were eluted with 0.5 mL 50% acetonitrile, 0.1% trifluoroacetic acid and vacuum dried to remove acetonitrile. Lysates were stored at −80°C until phosphopeptide enrichment. Phosphopeptide enrichment was performed with the High-Select^™TiO₂ Phosphopeptide Enrichment Kit (Thermo Fisher Scientific, Rockford, IL, USA). Eluted peptides were dehydrated using a speed-vac.

Nanoflow LC-MS/MS

Approximately 600 ng of peptide mixture, contained in a 1 μl aliquot of 98:2, water:acetonitrile, 0.1% formic acid, were analyzed in data dependent acquisition mode by liquid chromatography (LC)-nano ESI-mass spectrometry (MS) with an Ultimate 3000 Dionex RSLC nano LC system online with an Orbitrap Fusion Tribrid mass spectrometer (Thermo Fisher Scientific, Rockford, IL, USA). Peptides were separated during a linear gradient with the following profile: 5–22% solvent B in 70 min, 22–35% solvent B over 35 min, and 90% solvent B held for 10 min. Solvent A was water with 0.1% formic acid and solvent B was 80% acetonitrile with 0.1% formic acid. The 50 cm column was packed in house using ReproSil-Pur 120 C18-AQ 1.9 μm (Dr. Maisch, Ammerbuch-Entringen, Germany) in a PicoTip 75 μm inner diameter (New Objective, Littleton, MA, USA). Data acquisition was acquired with the following MS parameters: ESI voltage 2.1 kV, ion transfer tube 275°C; Easy-IC internal calibration; Orbitrap MS1 scan 120k resolution in profile mode from 380–1580 m/z with 100 msec injection time; 100% (4 × 10E5) automatic gain control (AGC); MS2 was triggered on precursors with 2–6 charges above 5 × 10E4 counts; MIPS (monoisotopic peak determination) was set to peptide; MS2 settings were: 1.6 Da quadrupole isolation window, higher energy collisional dissociation activation at 35% collision energy, Orbitrap detection with 60K resolution at 200 m/z, first mass fixed at 110 m/z, 150 msec max injection time, 100% (5 × 10E4) AGC and 40 sec dynamic exclusion duration with +/− 10 ppm mass tolerance.

Phosphoproteomics database search, phosphoprotein and phosphopeptide quantification

The raw MS files were processed by Proteome Discoverer v2.4 (Thermo Fisher Scientific, Rockford, IL, USA). MS/MS spectra were searched against the UniProtKB/Swiss-Prot mus musculus database (55,474 entries, UniProt UP000000589, downloaded November 2019) with the Sequest HT search engine embedded in Proteome Discoverer v2.4. Parameters were set as follows: MS1 tolerance of 15 ppm, MS/MS mass tolerance of 0.05 Da, trypsin (full) digestion with a maximum of two missed cleavages, minimum peptide length of 6 and maximum of 144 amino acids. Cysteine carbamidomethylation (57.02 Da) was set as a fixed modification, and methionine oxidation (15.99 Da), asparagine and glutamine deamidation (0.98 Da), acetylation of the N-terminus (42.01 Da), and phosphorylation of tyrosine, serine, and threonine (79.97 Da) were set as dynamic modifications. A false discovery rate (FDR) of 1% was set for peptide-to-spectrum matches using the Percolator algorithm (v3.02.1) and for protein assignment. Phospho-localization scoring was performed with the IMP-ptmRS v2.0 node and only phosphopeptides with a localization score >0.8 were used for quantification. Unique and razor peptides were used for quantification. Precursor abundance quantification was based on area and normalized by total peptide amount. All peptides were used for normalization for protein quantification; however, only phosphorylated peptides were used for pairwise ratios and protein roll-up.

Label-free quantitation (LFQ) of phosphopeptides were performed with normalized abundances using the Proteome Discoverer LFQ algorithms. The protein ratio was calculated as the geometric median of the phosphopeptide group ratios, and the phosphopeptide group ratios were calculated as the geometric median of all combinations of phosphopeptide ratios from all the biological replicates in the study. Phosphopeptides that had a high 1% FDR confidence were used for further analysis. We applied a maximum p-value filter of 0.05 and a minimum relative fold change of phosphoprotein and phosphopeptide expression at 1.4 to the Proteome Discoverer quantification results. Phosphoproteins and phosphopeptides were considered significant and differentially regulated if they had an adjusted p-value < 0.05 and were defined as downregulated if they had a fold change ≤ −1.4 or upregulated if they had a fold change ≥ 1.4.

Kyoto Encyclopedia of Genes and Genomes (KEGG), Reactome, and Gene Ontology (GO) enrichment analysis

Significant and differentially regulated phosphopeptides were mapped back to their precursor protein and the list of proteins was used for overrepresentation analysis using the clusterprofiler and ReactomePA packages in R v4.1.1 [63–65]. Comparative enrichment analysis of the Ovx/Sham and OA/YA datasets were performed for GO terms in all three domains: cellular component, molecular function, and biological process. KEGG and Reactome pathway analyses were performed to identify overrepresented pathways. Cluster analysis using the K-means algorithm was performed on the top 10 KEGG and Reactome pathways. Overrepresented pathways and GO annotation terms were considered significant if they had a Benjamini – Hochberg adjusted p-value < 0.05.

Ingenuity pathway analysis (IPA)

IPA (Qiagen, Redwood City, CA, USA) was used to perform core analyses on the on the phosphopeptides from Ovx/Sham and OA/YA datasets. A comparative analysis across both core analyses was performed in IPA using the activation Z-score algorithm to predict the activation states of pathways, functions, and upstream regulators from both datasets. The Z-score measures how closely the observed expression pattern of the molecules from the datasets compare to the expected expression pattern based on the literature for a particular annotation. Molecules from the dataset that met the cutoffs, |FC| ≥ 1.4 and adjusted p-value < 0.05 were considered for analysis.

Statistical analyses

Relative phosphoprotein quantification was analyzed in Proteome Discoverer v2.4 (Thermo Fisher Scientific, Rockford, IL, USA) using a background t-test. Fisher’s exact test was used in IPA to calculate p-values for the association or overlap between the identified molecules in the dataset and a given pathway/process/function. Benjamini-Hochberg post-hoc analysis was used to correct for multiple comparisons. The predicted activation state in the pathways/functions/upstream regulators in IPA were measured as a Z-score. Significant activation or inhibition was accepted at |Z| ≥ 2. Significance was accepted at α < 0.05 level.

Supporting information

Supplementary Materials and Methods and Figures can be found online.

The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via PRIDE [66] partner repository with the dataset identifier PXD035171 (http://www.ebi.ac.uk/pride/archive/projects/PXD035107).