Altered Tyrosine Phosphorylation of Cardiac Proteins Prompts Contractile Dysfunction in Hypertrophic Cardiomyopathy


 Altered Serine/Threonine phosphorylation of the cardiac proteome is an established hallmark of heart failure (HF). However, the contribution of tyrosine phosphorylation to the pathogenesis of these diseases remains unclear. The cardiac proteome was explored by global mapping to discover and quantify site-specific tyrosine phosphorylation in two cardiac hypertrophic models; cardiac overexpression of ErbB2 (TgErbB2) and cardiac expression of α-Myosin heavy chain R403Q (R403Q-αMyHCTg) compared to control hearts. Phosphoproteomic changes found in R403Q-αMyHC Tg mice indicated EGFR1, Focal Adhesion, VEGF, ErbB signaling, and Chemokine signaling pathways activity were likely to be activated. On the other hand, TgErbB2 mice findings displayed significant overrepresentation of Right Ventricular Cardiomyopathy, Hypertrophic Cardiomyopathy (HCM), and dilated cardiomyopathy (DCM) KEGG Pathways. In silico kinase-substrate enrichment analysis (KSEA) highlighted a marked downregulation of canonical MAPK Pathway Activity downstream of k-Ras in TgErbB2 mice and activation of EGFR, PP2 inhibition of c-Src, and Hepatocyte growth factor stimulation. In vivo ErbB2 inhibition by AG-825 decreased cardiac fibrosis, cardiomyocyte disarray, and rescued contractile function on TgErbB2 mice. These results suggest that altered tyrosine phosphorylation may play a regulatory role in cardiac hypertrophic models, suggesting that tyrosine kinase inhibitors could be used therapeutically in Hypertrophic Cardiomyopathy.

and discover potential targets for therapy [2][3][4] . For instance, PTMs of cardiac Troponin I (cTnI), a sarcomere protein that is centrally involved in myocardial contractility regulation, have been extensively studied, particularly for the functional role of PKA-dependent phosphorylation 5 . Notably, phosphorylation of cTnI-S22/23 -one of the most relevant regulatory sites of cTnI -is known to be down-regulated in human heart failure (HF) 6 and leads to contractile dysfunction 7 .
Tyrosine phosphorylation is essential for cardiac structural development and myofibril organization during embryogenesis 8,9 . For example, several tyrosine phosphatases have been linked to heart disease and even proposed as a therapeutic target for some conditions. Moreover, mutations at the level of tyrosine-protein phosphatase non-receptor type 11 (PTPN11) can lead to HCM or dilated cardiomyopathy (DCM) 10,11 . In the heart, PTPN11 likely plays a role in systolic dysfunction produced by pressure-overload 12 . Another tyrosine phosphatase is associated with cardiac pathophysiology, specifically acid phosphatase 1 (ACP1):-the deletion of ACP1 results in a protection against stress-induced cardiomyopathy 13 . Our group first showed that cTnI-Y26 phosphorylation is readily detected in healthy human hearts and downregulated in human HF and DCM 14 . However, how these alterations contribute to the onset and progression of cardiac disease remains poorly understood. A better grasp of additional site-specific changes of individual pTyr sites would help in this direction, substantially.
In the present study, a label-free and tandem mass tagging (TMT) quantitative global phosphotyrosine proteomics approach was applied to determine which sarcomeric sites have altered amounts of Tyr phosphorylation at specific sites in two unrelated models of HCM. The first model is secondary to the overexpression of the tyrosine kinase receptor ErbB2 15,16 ; the second recapitulates features of the human disease, more specifically, an R403Q mutation of the myofilament protein myosin heavy chain, distinctive of familial hypertrophic cardiomyopathy 17 . The rationale for this choice was to determine whether triggering HCM through different mechanisms elicits similar, pTyrrelated pathways/regulatory sites within the heart. In doing so, the manipulation of these Tyrosine phosphorylation sites can offer new opportunities for therapeutic targeting in different forms of human HCM

Choice of the Models
The approach combined examination of the global proteome and myofilament targeted tyrosine phosphorylation proteomics approach in two HCM mice models to identify the potential impact of altered Tyr phosphorylation on underlying disease mechanisms. One model has HCM caused by cardiac-specific overexpression of a tyrosine kinase receptor (ErbB2); the second is a classic model of a point mutation in a significant sarcomeric protein (R403Q-aMyHC mice). TgErbB2 mice initially develop a striking, concentric cardiac hypertrophy, which evolves to diffuse fibrosis and myocyte disarray 16 with HCM.
This line of mice also has abnormal calcium handling, and are prone to arrhythmias and develop Hypertrophic Obstructive Cardiomyopathy 2 . Similarly, R403Q-aMyHC mice reproduce human familial HCM by progressing from mild hypertrophy, mild fibrosis to frank myocyte disarray, HF, and arrhythmias 17 .

Immunoblotting Reveals Increased Cardiac Tyrosine Phosphorylation in HCM models
First, immunoblot analysis of the myocardium of TgErbB2 and Ntg mice for tyrosine phosphorylation was applied. In the whole heart homogenates of ErbB2 overexpressing mice, there was a global increase of Tyrosine phosphorylation compared to nontransgenic (Ntg) littermates (Fig.1A). These findings were corroborated with an immunoblot using pan-tyrosine phosphorylated antibodies detected by a fluorescent secondary and normalizing the tyrosine phosphorylation signal to Troponin I levels (Fig.1B).
Although c-Src (proto-oncogene tyrosine-protein kinase Src) is downstream of ErbB2 signaling, we did not anticipate its increased expression in TgErbB2 mice (Fig1.C, D) when normalized to GAPDH signals. The c-Src activity was enhanced in TgErbB2 mice compared to Ntg, as revealed by enhanced Tyr416 phosphorylation (Fig.1E). Altogether these data suggest that tyrosine phosphorylation is a broadly distributed post-translational modification (PTM) in the myocardium. Therefore, alterations in cardiac tyrosine phosphorylation may play a regulatory role in the disease progression of non-ischemic cardiomyopathies, such as familial cardiomyopathy (HCM) and dilated cardiomyopathy (DCM), as part of the pathophysiologic response to the underlying disease.

Tyrosine phosphorylation in myofilaments and cross-talk with serine/threonine kinases/phosphatases
A label-free phospho-tyrosine proteomics approach was undertaken to compare the cardiac tyrosine phosphorylation profile of Non-transgenic controls (Ntg), TgErbB2, and R403Q-aMyHC mice (Fig.2). It was hypothesized that unbiased global pTyr profiling would yield salient information about which specific tyrosine sites are phosphorylated in well-known essential cardiac-specific proteins and tyrosine kinases, thus providing clues about which tyrosine kinases could be mediating cardiac tyrosine phosphorylation. A total of 1,800 peptide spectra matched (PSM) were collected from the heart ventricle whole proteome, and 50% were identified with high confidence resulting in 431 unique tyrosinephosphorylated peptides corresponding to 239 proteins. Tables with all identified phosphoproteins and phosphopeptides are provided in Supplementary Table 1-7. An evaluation of the MS data quality, intensity, and distribution post-median sweep normalization are shown in Supplementary Fig.S1. Confidence of Phosphorylated site localization was evaluated for annotation, and a score of more than 49% was required.
This data set indicates that a comparable yield of peptide identification was achieved, reproducibility on the enrichment of pTyr peptides, and high-quality MS/MS data using similar approaches to others [18][19][20] .
The label-free phospho-tyrosine proteomics approach core findings identified 239 different tyrosine-phosphorylated proteins, 169 tyrosine-phosphorylated proteins that overlap among the control mouse hearts (Ntg), and the myocardium from two cardiac hypertrophic models (TgErbB2 and R403Q-aMyHC). At the peptide level, 431 pTyr sites were detected with an overlap of 271 tyrosine phosphorylated residues (Fig.2E, F). Also, 18 novel pTyr sites were detected with an overlap of 13 pTyrosine sites for the three mouse groups (Fig.2G). These are novel sites not reported before, nor found in the extensively curated database PhosSitePlus (Cell Signaling Technology).
Next, the focus was on the cardiac-specific proteins from the myofilament apparatus. Nine major myofilament proteins were tyrosine phosphorylated. They harbor multiple tyrosine phosphorylation sites, many of which are novel (See Supplementary Fig.S2); the best example is Titin with 36 pTyr sites (See Supplementary Fig.S3). Also, seventeen tyrosine kinases had detectable tyrosine phosphorylation sites, and thus they could potentially be involved in regulating cardiac myofilament tyrosine phosphorylation.
In contrast, only two tyrosine phosphatases (Ptpn11 and Ptpra) demonstrated detectable tyrosine phosphorylation sites and, therefore, are likely to play a role in regulating cardiac tyrosine phosphorylation levels in the sarcomere. In addition to the number of tyrosine kinases found with pTyr sites, fifteen serine/threonine kinases had detectable tyrosine phosphorylation sites; however, only one serine/threonine phosphatase (Ppp1r12b) demonstrated tyrosine phosphorylation sites. These data show that cardiac myofilament tyrosine phosphorylation is broadly distributed and that there may be a cross-talk between tyrosine kinases/phosphatases with serine/threonine kinases/phosphatases.

ErbB2 Cardiac Overexpression and R403Q-aMyHC Point Mutation Remodeled the Cardiac Tyrosine Phosphorylation Proteome
Next, studies were undertaken to define the impact of ErbB2 cardiac overexpression and R403Q mutation of the cardiac myosin heavy chain on the specific tyrosine phosphorylation changes. Data were then analyzed to determine the signaling transduction pathways associated with two unrelated types of HCM. Spectral normalization and statistical methods were performed, as previously described 21,22 .
Then, principal component analysis (PCA), k-means and hierarchical clustering were applied to determine if the pTyr profiles were similar within the groups. This approach determined that, indeed, the three groups segregate distinctively by the principal component 1 (PC1= 36.6%), whereas no group-wise distinction was evident in the second component (PC2=16.4%), see Fig.3A. Of note, one technical/biological replicate of R403Q-aMyHC was removed because it had more than 50% of absent pTyr peptides compared to the other two replicates. Hence, it was considered a technical failure.
Heatmaps of hierarchically clustered up-or down-regulated pTyr peptides helped visualize the technical reproducibility and the specificity of how the genotype of cardiac hypertrophy largely influenced the phospho-tyrosine proteome remodeling (Fig.3B, p<0.05 by ANOVA). Tyrosine phosphorylated peptides are color-coded according to their extracted chromatogram MS1 signal intensities; yellow is upregulated, and blue is downregulated. This approach is a relative measure for peptide abundance. Additionally, k-means, which is an unsupervised machine learning algorithm, discovered the underlying pattern, which is distinct to every group, see the clusters in Fig.3C   A pathway enrichment analysis was performed on significant genes/sites (Log2 FC>1 and p-value < 0.05) using a gene-protein level with MSigDB. It was found that TgErbB2 mice's most significant pathways by KEGG are right ventricular cardiomyopathy (ARVC), see In comparison, the PTMsigDB revealed that pathway analysis at the phosphorylated sitespecific level did not reach statistical significance for the EGFR1 pathway. In the case of TgErbB2 mice, only three site-specific phosphosites matched the EGFR1 pathway; they are highlighted in cyan color at the Volcano plot ( Fig.3D) and PTMsigDB table (Fig.3E).
Intriguingly, however, in R403Q-aMyHC mice, the same pathway was statistically significant; see the 16 site-specific changes that matched the EGFR1 pathway (p=0.001) represented in the volcano plot as cyan data points (Fig.3F) and PTMsigDB table (Fig.   3G). EGFR1 is a receptor tyrosine kinase, and the following most significant change is a signature secondary to Erythropoietin (PSP_EPO) stimulation (p= 5.6 -4 ); see three purple and one blue data points on the volcano plot (Fig.3F), the blue dot is to show an overlapping site, Stat5a Y694p, of EGFR1 pathway and PSP_EPO (Erythropoietin stimulation). These data show the utility of phosphorylation site-specific databases to narrow the search of a biologically relevant pathway on phospho-proteomics data sets, particularly in understudied phosphorylation events such as pTyr that have small data sets to date.

Tyrosine Phosphoproteome Dysregulation
A tandem mass tagging (TMT) quantitative labeling proteomics was used to gain more specific insight into the myofilament tyrosine phosphorylation changes. The hypothesis is that myofilament enrichment would enhance site-identification in myofilament proteins, specifically those with low abundance phospho-Tyrosine modifications, that may have been missed using the global approach. To that end, Ntg and TgErbB2 mice were evaluated using freshly isolated myofilaments, as described previously 29 with minor modifications, from three hearts per genotype and evaluated the tyrosine phosphoproteome ( Fig.4A, B), as described previously for the Heart Failure phosphoproteome 20 . Briefly, in this 6-plex TMT experiment, the peptides identified with single spectra were removed, which lead to 24,727 peptide spectra matched (PSM). After median sweep normalization, there were 1,116 peptides that corresponded to 1,092 proteins. For spectral intensity distribution, before and after normalization, see peptides was quantified because they had corresponding peptides from the full proteome's expression data. For further statistical data analysis, a workflow similar to that used for Label-Free proteomics was adopted. The only difference was that specific phosphorylation levels were calculated matching pTyr peptides intensity levels to their respective protein intensity levels on the flow-through fraction of peptides used for full proteome determination and unsupervised statistical analysis obviated k-means because it is designed for experimental groups of three or more.
The first step of the analysis uses unsupervised principal component analysis (PCA) and hierarchical clustering. PCA analysis showed the correlation between members in the same group and the ellipse represent the confidence interval of 95% (See

Kinase-Substrate Enrichment Analysis (KSEA) Implicated Downregulation of MAPK in TgErbB2 and Upregulation of EGFR in R403Q-aMyHCTg
Kinase-Substrate Enrichment Analysis (KSEA) was used to characterize genotypeinduced signaling changes by estimating the kinase's relative activity in TgErbB2 and R403Q-aMyHC mice using the global pTyr label-free and TMT quantitative sarcomere pTyrosine proteomics data, using their respective Ntg groups as a reference. KSEA 30 is a method that infers the kinases' differential activity based on the differential phosphorylation of its substrates and computes scores that reflect the directional change in each kinase's signaling. This method assumes that the differential activity of kinases is correlated with phosphorylation changes in its substrates. A positive score corresponds to a kinase with phosphorylated substrates in Tg mice relative to Ntg control. Likewise, a negative score is a hypophosphorylation in Tg relative to their Ntg control. The kinasesubstrate interaction data was downloaded from the PhosphoSitePlus (PSP) 31 website.
Next, the KSEA method was applied to all the pTyr sites identified in these experiments. Since the pathway enrichment analysis of pTyr TMT data did not yield further insights distinct from label-free pTyr data, a protein-protein interaction (PPI) network approach was employed. Pathway analysis can miss protein groups mainly because pathways algorithms are predefined and rigid. PPI might better capture signaling responses in these models that are not classically detected by pathway analysis. To do so, MoBaS Analysis 22,32 was employed, i.e., to identify densely-connected subnetworks that are related and might exhibit differential phosphorylation in Tg models. Several modules were identified and focus on the top two most statistically significant PPI modules from the label-free data set; for purposes of interpretation, they were designated as Module 1 and Module 2. The substrates of Modules 1 and 2 are illustrated in Fig.6B and Fig.6D -αMyHC n=352, Fig.2F). This evidence suggests that the differences observed can be attributed to genotype or phenotype differences between these two mouse models. Interestingly, we unearthed the same HCM/DCM tyrosine-phosphorylated peptides in TgErbB2 mice and R403Q-αMyHC mice. However, only on R403Q-αMyHC mice, the changes reached statistical significance for HCM/DCM KEGG pathways' involvement. Some of the pathways highlighted by the KSEA in R403Q-αMyHC Tg mice include the Focal Adhesion, VEGF, ErbB signaling, Chemokine Signaling, and Angiopoietin receptor pathways.
Initially, it was reasoned that overexpression of a receptor tyrosine kinase in TgErbB2 mice hearts would result in a more significant alteration of tyrosine phosphoproteome. However, the data suggest that a single sarcomere point mutation, as in the case of R403Q-αMyHC mice, is enough to alter the cardiac tyrosine phosphoproteome significantly. Another surprising finding was that the ErbB signaling pathway elements were not significantly altered in TgErbB2 mice, but rather in R403Q-αMyHC mice. Indeed, MSigDB analysis in the TgErbB2 mice yielded only two targets (CRK and SHC1) associated with the ERBB Signaling KEGG pathway. The TgErbB2 mouse model is important because it shows pronounced cardiac hypertrophy at an early age, sarcomere dysfunction, and abnormal calcium handling 16 , thus connecting the tyrosine kinase pathway and tyrosine phosphorylation to cardiac hypertrophy. While ErbB2 was previously reported to be overexpressed ~40-fold in TgErbB2 hearts 12 , the tandem mass tagging quantitative proteomics data, which allowed normalization of the peptide phosphorylation signal to total protein levels of ErbB2, revealed that the phosphorylation level of ErbB2-Y1006 peptide is reduced in this model (Fig.5G). The low stoichiometry of phosphorylated:non-phosphorylated ErbB2 may account for the apparent lack of alteration of the ErbB signaling pathway.
The mechanisms by which a pathogenic mutation in the myosin heavy chain may affect tyrosine phosphorylation regulation is unknown and warrants further investigation.
On the other hand, pharmacological inhibition of ErbB2 (HER2/neu) with the inhibitor lapatinib in breast cancer patients treated with doxorubicin increases the risk of developing heart failure compared to patients treated with doxorubicin alone 35 . This finding suggests that maintaining homeostasis of tyrosine phosphorylation may be important for regulating cardiomyocyte function and homeostasis.

EGFR1 pathway is central to HCM
The EGFR1 (ErbB1) Pathway seems to play a central role in cardiac hypertrophy in both Tg models. However, at the site-specific level, the molecular signature (PTMsigDB) of R403Q-αMyHC had 16 specific peptides altered consistent with the activation of the EGFR1 pathway. The TgErbB2 mice tyrosine-phosphorylated peptides did not reach statistical significance, with only three specific peptides for PTMsigDB (Fig.3E, G). EGFR pharmacological inhibition, using AG-1478, protects against Angiotensin II-induced cardiac hypertrophy in vitro and in vivo 36 . The concentric hypertrophy associated with ErbB2 cardiac-specific overexpression can be reversed by early administration of Lapatinib, which inhibits EGFR receptor Tyrosine kinase 12 in addition to ErbB2. In this study, the ErbB2 receptor's pharmacological blockage by the small molecule inhibitor AG-825 in TgErbB2 adult mice preserved cardiac function without affecting cardiac hypertrophy. In addition, the specific blockage of ErbB2 by AG-825 reduced myocyte disarray and cardiac fibrosis. These findings indicate some overlap in cardiac hypertrophy signaling mediators downstream of EGFR1 and ErbB2, and the sarcomeric mutation R403Q-αMyHC appears to amplify the EGFR Signaling Pathway preferentially.

Tyrosine phosphorylation is present in multiple proteins associated with the regulation of cardiac function and disease
A wide variety of other functionally relevant targets were identified in this study as tyrosine phosphorylated, from sarcomere proteins to Z-disk and desmosome components, focal adhesion, and adherence junction proteins, membrane receptors, kinases, and phosphatases. For instance, tyrosine phosphorylation in several Z-disk associated proteins was noted (For a full list, see Supplementary tables 2-7). Z-disk proteins are crucial for muscle contraction and mechanical stress, growth, and metabolic signaling 37 .
Also, the alteration in the phosphorylation levels of the desmosome key component Plakophilin-2 protein (Pkp2) peptide Pkp2-Y119 was noted in both mouse models, upregulation in R403Q-aMyHC Tg and down-regulation in TgErbB2. The functional effect of up-or down-regulation of Pkp2-Y119 phosphorylation is not known. However, Pkp2 homozygous deletion disrupts heart architecture and is lethal in the embryo 38 . In the heart, Pkp2 is required for the assembly of the desmosome and the PKC activity 39 .
Autosomal dominant mutations in this gene are responsible for 25% to 50% of arrhythmogenic cardiomyopathy (ARVC). Interestingly, TgErbB2 mice have increased susceptibility to arrhythmias and myofibrillar disarray 16 , similar to patients with myosin mutations and HCM. Pkp-Y119 phosphorylation changes could impact the phenotype.

Bioinformatic Identification of Potentially Activated Kinases that are Drug Targets
The The MAPK pathway is involved in the adaptive and maladaptive response that could lead to heart hypertrophy (for review 40 ). The RAS-RAF-MEK-ERK signaling pathway is an attractive target for therapeutic intervention in oncology, and several selective RAF and MEK small molecule inhibitors have been tested in clinical trials 41 . The present study found increased phosphorylation at MAPK1-Y185, which activates the kinase, in R403Q-αMyHC, and MAPK1 activity directly relates to heart hypertrophy 42 . c-Src figures as one of the potential regulators for the phosphotyrosine proteome changes observed in R403Q-αMyHC because it phosphorylates the EGFR receptor upstream and phosphorylates Stat, among other targets, downstream. c-Src affects the response to mechanical cardiomyocyte stretching by triggering a cascade of intracellular signaling in cardiomyocytes towards a hypertrophic response 43,44 . Also, c-Src phosphorylates PXN-Y118 45 , and this site had a 25-fold increase in phosphorylation in R403Q-αMyHC mice.
Interestingly, c-Src mediates the activation of MAPK1 and MAPK3 in response to pressure overload 43 . Notably, a previous study has shown that pressure overloadinduced cardiac hypertrophy is exacerbated in R403-αMyHC Tg mice 46 , suggesting that the mechanism by which R403-αMyHC mutation produces heart hypertrophy might be by sensitizing cells to pressure-overload induced signaling via c-Src.
R403Q-αMyHC mice displayed enhanced JAK2, STAT5A, and STAT5B phosphorylation in the activation sites (Y570, Y694, and Y699, respectively), which indicates not only activation of MAPK signaling but also activation of JAK-STAT signaling. Jak-Stat Signaling: IL-6 pathway is activated in response to the IL-6 family of cytokines (IL-6, cardiotrophin 1, and leukemia inhibitor factor), cross-talks with the EGFR pathway, and is involved in cardiac hypertrophy. This pathway has cardioprotective effects, but chronic activation may lead to heart hypertrophy (for review 40

Limitations & Studies in Perspective
The studies were limited to two transgenic mouse models that develop cardiac hypertrophy; other models of sarcomere mutations or pressure overload such as aortic banding could be explored and compared. Some relevant tyrosine-phosphorylated sites might be missed due to their low stoichiometry and technical challenges. This problem is particularly true for membrane-bound proteins, which are difficult to evaluate in phosphoproteomic studies.

Conclusions
This study reveals that altered patterns of Tyrosine phosphorylation are striking in two separate models of hypertrophic cardiomyopathy. Moreover, despite some shared sites, etiologically different forms of HCM may harbor specific tyrosine phosphorylation signatures. This evidence, combined with inhibition of a specific receptor tyrosine kinase using tyrphostin AG-825, can reverse cardiomyocyte disarray, thus preserving function in a TgErB2 mouse HCM model, could rationalize approaches to manipulate the Tyr phosphoproteome as a therapeutic approach for HCM in general. Furthermore, these studies indicate that tyrosine kinase inhibitors, now used broadly in cancer therapies, may change cardiac function by directly modifying the heart's tyrosine kinase profile.

Western blot
Whole heart lysates from Heat Stabilized Tissue were resuspended in 1% SDS buffer, Signaling Technology) was mixed with peptide solution and incubated on a rotator at 4°C for 2 h. After incubation, the beads were washed with IP buffer and water for two and three times, respectively. The phosphotyrosine peptides were eluted using 0.1% TFA. The eluted peptide samples were desalted using C18 STAGE tips, vacuum dried, and kept at −80°C before LC-MS analysis. From flow-through, a quantification of the full proteome was made to serve as a reference for pTyr proteome findings.

LC-MS/MS Analysis: The phosphotyrosine peptides were analyzed on Orbitrap Fusion
Lumos Tribrid (Thermo Scientific, San Jose, CA, USA) coupled with Easy-nanoLC 1200 nanoflow liquid chromatography system (Thermo Scientific). The peptides from each fraction were reconstituted in 10% formic acid and loaded on Acclaim PepMap100 Nano Trap Column (100 μm × 2 cm) (Thermo Scientific) packed with 5 μm C18 particles at a flow rate of 5 μl per minute. Peptides were resolved at 250 nl/min flow rate using a linear gradient of 10 to 35% solvent B (0.1% formic acid in 95% acetonitrile) over 95 minutes on the EASY-Spray column (50cm x 75µm ID, PepMap RSLC C18, and 2µm C18 particles) (Thermo Scientific) and it was fitted on EASY-Spray ion source that was Two missed cleavages were allowed and 'match between runs' was enabled. Peptides and proteins were filtered at a 1% false discovery rate. Proteins with a q-value lower than 0.05 will be considered as differentially expressed with statistical significance. In a complementary way, R software was used.
A two-sided t-test and the moderated p-values and q-values calculation was performed using an algorithm developed by Herbrich et al 50 .

Systems Biology Analysis
Statistical analysis of the label-free data set was performed using Perseus software 51 .
Peptides with less than 50% of detection were filtered and the missing values of the intensities of the reminding peptides were filled with k-nearest neighbors. Then, the log2 intensities were normalized by mean subtraction subsequently ANOVA and false discovery rate calculations. The K-mean clustering with Spotfire, 4 clusters were made using correlation as similarity measurement and data centroid based search. A heatmap was used to visualize the unsupervised hierarchical clustering of proteins with p-value <0.05 by ANOVA test.
The web-based software Ingenuity Pathway Analysis (IPA) was used for both data sets for pathway analysis. The results obtained were uploaded to IPA software. The Uniprot Accession ID was selected as an identifier for each site, and the cut of values for each intensity was a p-value < 0.05 and an absolute fold change of 2. Interaction networks were generated, including endogenous chemicals with 70 molecules per network and 25 networks per analysis. All node types and data sources were selected, considering all experimentally observed and predicted (high and moderate confidence) in all species.

Pathway Enrichment Analysis:
The hypergeometric p-value was used to identify the processes that are significantly enriched in the proteins and phosphosites that are differentially phosphorylated between case and control samples. For this purpose, MsigDB 24 was used to analyze the protein level pathways and PTMsigDB 25 to analyze the processes at the site level. The population/background for which the enrichment is calculated against is restricted to all the proteins in which their sites are quantified in the phosphoproteomics experiment, rather than all universally known proteins/phosphosites.

Kinase-Substrate Enrichment Analysis (KSEA): Kinase Substrate Enrichment
Analysis seeks to identify kinases whose targets exhibit significantly altered phosphorylation levels in a given condition. KSEA scores each kinase k with a set of substrates S as follows:  Whole Heart Lysates n=5 Ntg mice