Integrated in silico MS-based phosphoproteomics and network enrichment analysis of RASopathy proteins

Background RASopathies are a group of syndromes showing clinical overlap caused by mutations in genes affecting the RAS-MAPK pathway. Consequent disruption on cellular signaling leads and is driven by phosphoproteome remodeling. However, we still lack a comprehensive picture of the different key players and altered downstream effectors. Methods An in silico interactome of RASopathy proteins was generated using pathway enrichment analysis/STRING tool, including identification of main hub proteins. We also integrated phosphoproteomic and immunoblotting studies using previous published information on RASopathy proteins and their neighbors in the context of RASopathy syndromes. Data from Phosphosite database (www.phosphosite.org) was collected in order to obtain the potential phosphosites subjected to regulation in the 27 causative RASopathy proteins. We compiled a dataset of dysregulated phosphosites in RASopathies, searched for commonalities between syndromes in harmonized data, and analyzed the role of phosphorylation in the syndromes by the identification of key players between the causative RASopathy proteins and the associated interactome. Results In this study, we provide a curated data set of 27 causative RASopathy genes, identify up to 511 protein–protein associations using pathway enrichment analysis/STRING tool, and identify 12 nodes as main hub proteins. We found that a large group of proteins contain tyrosine residues and their biological processes include but are not limited to the nervous system. Harmonizing published RASopathy phosphoproteomic and immunoblotting studies we identified a total of 147 phosphosites with increased phosphorylation, whereas 47 have reduced phosphorylation. The PKB signaling pathway is the most represented among the dysregulated phosphoproteins within the RASopathy proteins and their neighbors, followed by phosphoproteins implicated in the regulation of cell proliferation and the MAPK pathway. Conclusions This work illustrates the complex network underlying the RASopathies and the potential of phosphoproteomics for dissecting the molecular mechanisms in these syndromes. A combined study of associated genes, their interactome and phosphorylation events in RASopathies, elucidates key players and mechanisms to direct future research, diagnosis and therapeutic windows. Supplementary Information The online version contains supplementary material available at 10.1186/s13023-021-01934-x.


Background
RASopathies are a group of phenotypically overlapping syndromes caused by germline mutations that encode components of the RAS/MAPK signaling pathway, affecting growth and development [1,2]. The RAS/ MAPK signaling pathway is a chain of proteins in the cell that communicates a signal from a receptor on the surface of the cell to the DNA, and has major impact in human health [2][3][4]. These disorders include neurofibromatosis type 1 (NF1), Legius syndrome (LS), Noonan syndrome (NS), neurofibromatosis-Noonan syndrome (NFNS), Noonan syndrome-like (NSL), Noonan syndrome with multiple lentigines (NSML), formerly known as LEOPARD syndrome, Noonan syndrome-like with loose anagen hair (NSLSH) also known as Mazzanti syndrome, Costello syndrome (CS), cardiofaciocutaneous (CFC) syndrome, capillary malformation-arteriovenous malformation syndrome (CM-AVM), intellectual disability associated with autism spectrum disorder and juvenile myelomonocytic leukemia (JMML) [3][4][5]. To date, mutations in 27 genes have been proven as cause of these RASopathies (Fig. 1).
Phosphoproteomics is a valuable approach to understand diseases linked to phosphorylation events, and it has a potential value for the clinics and for a personalized medicine [83,84]. Classical immunoblot assays have been used to study phosphorylated amino acids in RASopathy proteins [85][86][87]. Whereas with the advent of more sensitive and accurate liquid chromatography-mass spectrometry (LC-MS) instrumentation, the number of large-scale mass spectrometry-based phosphoproteomic studies has swelled over the past decade [88], initiating its application to RASopathies [85,[89][90][91][92]. These techniques should allow to pinpoint biologically relevant phosphorylation events of utmost importance in the mechanisms of RASopathies [5]. However, only a few studies focused on their dysregulated phosphoproteins and their biological implications. Approaches mimicked mutations found in human samples associated to NS, NSML, and NF1 [85][86][87][89][90][91][92][93][94]. In NF1 studies, NF1null cells derived from malignant peripheral nerve sheath tumors (MPNST) showed RAS cascade hyperactivation typical of NF1. The interplay of H-Ras, N-Ras, and K-Ras produced a counter-inhibition of these classical Ras proteins, revealing an intense phosphoproteome remodeling [89]. Other studies were also conducted focusing on proteins such as CRMP-2 [93] and dynein IC2-C [85]. RAF1 [86], BRAF [94], PZR and SIRPa [87] are related to both NS and NSML, and changes in their phosphorylation were also observed. High-throughput screening strategies based on mass spectrometry have also been used to study NS and NSML [90][91][92], revealing dozens of dysregulated phosphoproteins.
The current state of knowledge about the RASopathy interactome is mainly based on an integrated network presented at genome, interactome, and phenome levels [1]; Twelve causative genes and clinical symptoms were collected from OMIM and NCBI GeneReviews databases for 6 syndromes: NS, NSML, NF1, CFC, Legius and Costello syndrome. In particular, they created an interactome network based on interactions between 12 proteins (PTPN11, SOS1, RAF1, BRAF, KRAS, NRAS, HRAS, MAP2K1, MAP2K2, SPRED1, NF1, and RIT1) of the RASopathies and another one based on their 10 first neighbors. However, we still lack a comprehensive view about the protein-protein associations involved in all RASopathy proteins and their neighbors. Mutations in some of these genes may drive RASopathies, cancer, or both in the same patient, but that is something that needs to be empirically tested. Phosphoproteomics data in RASopathies is still scarce and its potential to decipher crucial signaling pathways involved in this family of disorders needs to be considered. The application of LC-MS analysis, relying on high-sensitive nanoflow setups and phosphoprotein enrichment, is encouraged to unveil underlying molecular mechanisms affected in these syndromes, and suggest therapeutic targets for clinical implementation.
In this work, we aim to provide an up to date panel of proteins underlying RASopathies in order to identify protein-protein associations, and assign them to their cognate syndromes. The interactome was further analyzed in silico by binary interactions between the RASopathy proteins, a search for main hub proteins in the interactome and their classification based on GO terms for biological processes. We also carried out an analysis of phosphorylation level changes in the proteins from the RASopathy interactome integrating previous phosphoproteomic studies in the context of RASopathy syndromes. Finally, we analyzed the phosphosite abundance in causative RASopathy proteins in silico.

RASopathy protein-protein associations
In order to identify protein-protein associations within the RASopathy proteins we used STRING database to download a subgraph composed of the initially 32 RASopathy proteins selected in this study. Among the panel of 32 RASopathy proteins, our results show that 6 proteins belong to the RAS family, 15 are directly associated with RAS, and 11 proteins do not seem to be associated with RAS ( Fig. 2A). RAS includes the classical RAS (H-RAS, K-RAS and N-RAS) and the three RRAS (RRAS, RRAS2 and RRAS3/MRAS). The 11 proteins that do not directly associate with RAS are ANKRD11, MEF2C, SHOX and SRCAP (aforementioned as non-classical RASopathy proteins), KAT6B, CBL, MAP3K8 and A2ML1 (which cause NSL), LZTR1 and RIT1 (responsible for NS), and PPP1CB (which causes NSLSH). The proteins that associate with all RAS include but are not limited to BRAF, RAF1, as well as several guanine nucleotide exchange A B  Table S1) factors (GEFs) and GTPase-activating proteins (GAPs) that regulate RAS activity ( Fig. 2A). We then generated an interactome for the 32 RASopathy proteins including their direct interacting partners. This interactome contains a total of 765 proteins. Some of them were found associated to more than one RASopathy protein and therefore they appeared duplicated. When eliminating the duplicated proteins, we ended up with 511 unique proteins based on unique UniProt IDs (Additional file 1). On average, the overlap of the interactome among the 32 RASopathy proteins is 56%. According to our results, RRAS2 interactome fully overlaps with other RASopathy protein's interactome, followed by SHOC2 (90%), MEK2 (87%) and RAS (87-53%) among other proteins (Fig. 2B). Interestingly, FGFR3 did show direct interaction with RAS, was previously shown to exert and impact on RAS-MAPK signaling pathway [73][74][75] and had an interactome overlap of 24% with other RASopathy protein neighbors (Fig. 2B). On the other hand, ANKRD11, MEF2C, SRCAP, SHOX and KAT6B does not interact with RAS/MAPK and they show a low overlap (Fig. 2B). Therefore, they were not considered for further analyses in this study. However, A2ML1 was not excluded because according to our results A2ML1 interactome has a protein-protein association overlap with the rest of the RASopathy interactomes above 25%, which is even higher than the one obtained for other RASopathy proteins. Also, A2ML1 may interacts with the RAS-MAPK signaling pathway; A2ML1 is known to bind to lipoprotein receptor-related protein 1 activating of the Ras/MAPK pathway through its association with SHC domain proteins and CBL during recruitment to the plasma membrane [95,96]. Considering the above data, from the initial 32 RASopathy proteins we continue studying 27 proteins, excluding ANKRD11, KAT6B, MEF2C, SHOX and SRCAP from further analyses.

The interactome of RASopathy proteins by syndrome
The interactome of 27 RASopathy proteins including their direct interacting partners, yield a total of 687 proteins (432 unique entries based on UniProt IDs). These proteins were assigned to their cognate syndromes in order to identify their occurrence by syndrome, purely reflecting the number of proteins associated with the syndrome (Fig. 3A). Our results show that the largest overlap correspond to NS with a total of 196 proteins (41.2%) shared with other RASopathy interactomes. On the other hand, FGFR3-related syndromes had the smallest overlap with only 4 proteins (0.8%) (Fig. 3A).
The analysis of all possible binary interactions between the interactome of each of the RASopathy proteins shows that 124 out of 351 binary interactions do not have in common a single protein (Additional file 2). Our results also show that 185 binary interactions share between 1 and 5 proteins, 38 between 6 and 10, and only 4 binary interactions share more than 10 proteins in their interactomes. This last set of binary interactions was found between HRAS and NRAS with 15 proteins, followed by BRAF and RAF1 with 14 proteins, HRAS and KRAS with 11 proteins, and HRAS and RRAS2 with 11 proteins (Additional file 2).
A search for main hub proteins in the interactome resulted in 12 nodes. The protein corresponding to the node with most associations with other genes within the RASopathy interactome is AKT1 with 199 associations reported. The second node with more interactions, 186 edges, corresponds to HRAS. PIK3R1 is the third with 124 known interactions. SRC yielded 118 interactions, SOS reported 109 interactions, GRB2 has 104, and CDC42 has 100 interactions. Other highly connected nodes found on the interactome network were NRAS with 94 interactions, SHC1 with 87, PPP2CA reporting 77, BDNF with 76 and JNK1 with 76 interactions.
An enrichment analysis based on GO terms for biological processes was done using the 432 unique proteins (Fig. 3B). Interestingly the largest group of proteins correspond to proteins with phosphorylated tyrosine residues. Followed in abundance by PKB, RAS and MAPK signaling proteins, phosphatase activity and PIP3 biosynthesis and signaling proteins (Fig. 3B). We also found that several biological processes directly related with the nervous system may be altered, such as axonogenesis, peripheral nervous system development and neuron apoptosis.

Analysis of dysregulated phosphoproteins in RASopathies
Functional annotation clustering of the interactome of the 27 RASopathy proteins and their interacting partners shows that 370 out of 432 proteins are annotated as phosphoproteins (Additional file 3). In order to study dysregulated phosphoproteins in RASopathies, protein phosphorylation level changes associated to RASopathies were compiled (Additional file 4) from nine phosphoproteomic (Fig. 4A) and twenty-seven immunoblotting studies ( Table 1). The analysis of phosphorylation protein level changes within the RASopathy interactome comprises a total of 37 phosphoproteins (Additional file 5). In total and by syndrome, phosphoproteomic studies show that 8 phosphosites from 3 proteins are dysregulated in NF1, 53 phosphosites from 21 proteins in NS, 75 phosphosites from 31 proteins in NSML, 1 phosphosite in NSL, 2 phosphosites from 3 proteins in NSLH, 3 phosphosites from 2 proteins in LS, 4 phosphosites from 2 proteins in CS, 9 phosphosites from 3 proteins in CFC and 2 phosphosites from 2 proteins in JMML (crossing Additional file 1 and Additional file 4). Immunoblotting results show that dysregulation on several phosphoproteins account for more than one RASopathy, including ERK1/2 (Thr202/185 and Tyr204/187), AKT1 (Thr308 and Ser473), MEK1 (Thr292) and RAF1 (Ser259), which are the most co-occurring studied proteins along the syndromes (Table 1). For most RASopathies, phosphorylation of downstream effectors was upregulated compared to the control, including but not limited to AKT, MEK1 and ERK1/ERK2.
The most studied syndromes are NS and NSML, thanks to the availability of phosphoproteomics data. In NS the average fold change in upregulation is 2.61 and in downregulation 0.29, the highest fold change corresponds to PZR:Y241 (7.72) and the lowest to PTPN11:Y63 (0.02) with nearly a complete depletion of phosphorylation. In NSML the average fold change for upregulation is 2.29 and for downregulation 0.80, the highest value also corresponds to PZR:Y263 (6.7) whereas the lowest corresponds to PDGFR (0.02). These values demonstrate  Clustering analysis of the dysregulated phosphosites identified in phosphoproteomics and within the RASopathy interactome in either NS or NSML. An asterisk is added to zero values representing reported non-altered phosphorylation levels, to distinguish them from those representing just lack of data. In case of duplicated quantitative values from different studies, ' symbol was added [91]. B GO Biological processes overrepresentation test results using the 33 unique UniProt IDs from human dysregulated phosphosites reported in the RASopathy interactome on PANTHER  Figure 4 shows a comparison of quantitative values available for NS and NSML, that fulfil a threshold of 1.5 in both downregulation and upregulation events. Clustering analysis of the dysregulated phosphosites identified in phosphoproteomics and within the RASopathy interactome in either NS or NSML (Fig. 4A)

Analysis of phosphosite occurrence in RASopathy proteins
Studies of the phosphoproteome in RASopathies are scarce, making difficult to evaluate how phosphorylation affects proteins. In order to have an overview of phosphorylation sites in RASopathy proteins, we used Phosphosite database to consider all documented phosphosites rather than only those with studied dysregulation. We analyzed the phosphosites for the 27 RASopathy proteins and their cognate syndromes. Information relies on low-throughput analysis (a total of 1041, 5.6%) and high-throughput analysis (17,656, 94.4%). Using this information, our results show that all RASopathies are associated to at least one phosphoprotein (Fig. 5). NS, with 13 proteins, contains the highest number of RASopathy proteins susceptible of phosphorylation, followed by NFNS (5 proteins) and JMML (5 proteins). Interestingly, all RASopathy proteins have residues potentially subjected to regulation by phosphorylation. BRAF, CBL, NF1, PTPN11, RAF1, SOS1 and SOS2 contain at least 20 phosphosites each, and they account for 326 of the

Discussion
In this study we provide a global comprehensive picture of phosphoproteome remodeling in RASopathies, including but not limited to the most up-to-date information regarding the genes that cause the RASopathies. 5 (ANKRD11, MEF2C, SRCAP, SHOX and KAT6B) out of the initial 32 proteins were excluded from further analysis based on two factors: no direct association with RAS/MAPK signaling pathway and having an interactome with a protein-protein association overlap with the rest of the RASopathy´s interactomes below 16%. A2ML1 was not excluded because according to our results A2ML1 interactome has a protein-protein association overlap with the rest of the RASopathy´s interactomes above 25%, which is even higher than the one obtained for other RASopathy proteins. FGFR3 directly interacts with RAS, and our results suggest that 24% of its interactome is shared with other RASopathy proteins. This 24% is higher than the one observed for other RASopathy proteins, including but not limited to PPP1CB, and similar to NF1. FGFR3 was considered responsible for mosaic RASopathies [73], but has not been included in any RASopathy protein panel yet in clinical studies. FGFR3 has been found to influence the RAS/MAPK signaling pathway in many studies. FGFR3 mutations are frequent in superficial urothelial cell carcinoma (UCC) as oncogenic activation of FGFR3 is predicted to result in stimulation of the MAPK pathway [97]. MIP-1α (CCL3) is a downstream target of FGFR3 and RAS-MAPK signaling in multiple myeloma [98]. In hematopoietic cells, FGFR3 activates RSK2 to mediate activation of the MEK/ ERK pathway [99]. Mutations in FGFR3 and PIK3CA, singly or combined with RAS and AKT1, are associated with AKT but not with MAPK pathway activation in urothelial bladder cancer [100]. However, in another study it was found that enhanced activation of FGFR3 is linked to Ras and MAPK activation; in particular the authors described a novel FGFR3/Ras mediated mechanism for acquired-resistance to B-RAF inhibition [101]. RAS and FGFR3 mutations in urothelial carcinoma are mutually exclusive and non-overlapping events which reflect activation of oncogenic pathways through different elements. Also, papillary cancers typically exhibit activation of the MAPK pathway, as a consequence of oncogenic mutations in FGFR3 or RAS genes [102,103]. Forced expression of FGFR3 mutants in NIH-3T3 cells resulted in cellular transformation and mitogen-activated protein kinase (MAPK) activation, resembling the transfection effects observed with activated HRAS [104,105].
Activated FGFR3 seemed to be linked to RAS through adaptor proteins (that is, growth factor receptor-bound protein 2 (GRB2)-son of seven less (SOS) complexes) that are common to the RTK activation pathway [102]. These strong evidences of FGFR3 implications in the RAS/MAPK signaling pathway along with our in silico results of the interactome, let us suggest that FGFR3 disorders (thanatophoric dysplasia, achondroplasia, and hypochondroplasia) may be included among the RASopathies. On the other hand, KAT6B has no direct RAS interaction, poor direct association with RAS/MAPK signaling pathway with just one study so far [52] and its interactome is not associated at all with any other RASopathy protein. A translocation breakpoint 10q22.3 in a clinically diagnosed NS individual identified disruption of the KAT6B gene and hyperactivated MAPK signaling in humans and mice [52]. In particular, they found that functional studies using a patient-derived lymphoblastoid cell line causing KAT6B haploinsufficiency demonstrated an increase of RAS/MAPK pathway activity. Therefore, the authors postulated that altered expression of multiple genes associated with RAS/MAPK pathway regulation may be responsible for the increase in pathway activation and the NS-like phenotype [52]. However, that correlation has not been confirmed with more research so far. In a more recent study of a de novo heterozygous variant within exon 16 of KAT6B that was detected in a 7-months-old Chinese female infant, the patient presented symptoms of short stature, global developmental delay, and clinical features consistent with blepharophimosis mental retardation syndromes (SBBYSS, also called Ohdo syndrome) [106]. From a clinical point of view, the KAT6B phenotypic spectrum is broad and naming of the KAT6B-related disorders has been problematic, and suggested that considering this whole group as 'KAT6B spectrum disorders' may be more helpful [107]. As expected, RAS proteins are the ones with the highest number of associations, although RASopathies related to SPRED1, MEK1/2, PTPN11, FGFR3 and SPRY1 may regulate RAS signaling differently, affecting primarily the classical RAS versus RRAS signaling. We believe that further analysis is needed to determine the biological roles of the master nodes that we have identified in the RASopathy interactome. The functional characterization of these 12 nodes (AKT1, HRAS, PIK3R1, SRC, SOS, GRB2, CDC42, NRAS, SHC1, PPP2CA, BDNF and JNK1) may uncover implications and common patterns on the RASopathies signaling pathways, and their suitability as druggable targets. AKT is involved in many processes, including metabolism, proliferation, cell survival, growth and angiogenesis [108][109][110][111]. HRAS mutations are known to cause Costello syndrome [8,11,59]. Phosphatidylinositol 3-kinase plays an important role in the metabolic actions of insulin, and a mutation in this gene has been associated with insulin resistance [112]. SRC proto-oncogene may play a role in the regulation of embryonic development and cell growth [113]. GRB2 provides a critical link between cell surface growth factor receptors and the Ras signaling pathway [114,115]. Cell division cycle 42 (CDC42) regulates signaling pathways that control diverse cellular functions including cell morphology, migration, endocytosis and cell cycle progression [116][117][118]. Interestingly, some of these nodes, such as SHC1:(SHC1:Y427), AKT1 (S473, T308) and JNK1 (T183, Y185) are also subjected to dysregulated phosphorylation according to our phosphoproteomic analysis.
Many of the signaling cascades associated to the RAS/ MAPK pathway imply phosphorylation and dephosphorylation of proteins [5], and our results suggest that phosphoproteome remodeling in RASopathies is strongly relevant. The study of phosphosites is based on residues present in human and non-human model organisms that are extrapolated to the human sequences, so that dysregulated phosphosites may vary in human models and must be confirmed. A high proportion of the interactome (around 72%) are phosphoproteins, when compared to a 10% accepted value of phosphoproteins in the cytoplasm [89]. Further, all RASopathy proteins contain annotated phosphosites. Also, the gene ontology analysis of the RASopathy interactome and their dysregulated phosphoproteins shows mainly phosphorylation-related processes highlighting the PKB and MAPK pathways. Previous phosphoproteomics studies in RASopathies have revealed dysregulated phosphoproteins. Our results highlight several phosphoproteins dysregulated in NS, NSML and NF1 syndromes. The impact of a particular NS or NSML associated variant on downstream targets will depend first on whether the mutated RASopathy gene is linked particularly to NS (LZTR1, MRAS, NRAS, RASA2, RIT1, RRAS2, SOS1, SOS2) or both NS and NSML (BRAF, MEK1, PTPN11, RAF1). Second, mutations within the same causative gene linked to both NS and NSML may or may not have a different impact on downstream targets. Phosphoproteomic studies suggest that a higher proportion of phosphosites is upregulated (147 versus 47 downregulated), whereas immunoblotting suggest that most studied downstream effectors are upregulated, including but not limited to AKT1, BRAF, ERK, FAK, JNK, MEK and RAF. However, in some studies MEK and RAF have been also found downregulated in both syndromes. The relative variance in symptoms between RASopathies caused by activating signaling components in NS and NSML may suggest that there are different mechanisms at play. Systematic investigation of mutation strength across other pathway components rather than RAS, ERK or AKT may yield further insights into disease etiology. Subsequent studies on candidates highlighted by phosphoproteomics, as PZR protein in NS and NSML, have been already proved useful to explain crucial pathophysiological events, such as hypertrophic cardiomyopathy and cardiac fibrosis in NSML [119]. Phosphorylation of the RASopathy proteins RAF1 [86], PTPN11 [91,92], and BRAF [94] is altered in NS and NSML due to mutations in their corresponding genes, whereas RASA1 altered phosphorylation is associated to mutant SHP2 models [92] and not to its own mutation. The most upregulated phosphoproteins from the RASopathy interactome in NS and NSML, except for CRKL, have already been proposed as key players in these syndromes, namely PZR [90] and Fer kinase [91]. Similarly, downregulated phosphorylation levels of DPLY2 in NF1 have been well studied in neuronal differentiation [93]. Remarkably, contradictory values for identical phosphosites between different studies is an intriguing reason to further evaluate.
Noteworthy, most phosphosites in RASopathy proteins are reported on high-throughput mass spectrometry analysis. Protocols rely on different enrichment methods for phosphoproteins (once trypsinization, S-alkylation and reduction are performed). Nano-LC-MS/MS with reversed phase columns efficiently detects phosphopeptides and the modified amino acid, although occasionally this information is limited. Strategies of isobaric labelling (TMT, iTRAQ, SILAC) offer relative quantification of phosphorylated forms or, at least, qualitative information. For instance, the iTRAQ quantifications on SHP2 mutants [90][91][92] and immunoprecipitated BRAF phosphomapping [94] for NS and NSML, or SILAC strategies in NF1-null cells [89]. On the other hand, lowthroughput analyses include immunoblots with antibodies designed towards specific phosphoresidues in purified proteins as documented for RAF-1 [86] or CRMP-2 [93]. Commercially-available kits for western-blotting of phosphorylated proteins as ProQ-Diamond dying have also been used [93], although information of the phosphoresidue must be obtained somehow else. This is the case of the two undefined phosphosites in PZR and SIRPa in our dataset. High-throughput mass spectrometry analysis broadens, with its sensitivity limitations in the pmol range, to the examination of the whole phosphoproteome, whereas immunoblotting is limited to a few chosen proteins. Nonetheless, immunoblot analAysis offer specific and reliable information that may complement the monitorization of hundreds of proteins identified by MS strategies.

Conclusions
Herein, we provide for further studies with a comprehensive up-to-date analysis of 27 RASopathy genes involved in the RAS-MAPK pathway and their interactome, which comprises at least 432 different proteins including 12 nodes, highlighting AKT, HRAS and PIK3R1 as the top three, and involved in several biological processes directly related with the nervous system including but not limited to axonogenesis, peripheral nervous system development and neuron apoptosis.
Current state-of-the-art allows obtaining rich information based on the hypothesis that phosphoproteome remodeling is a key event in RASopathies. Accordingly, we present an integrative and comprehensive analysis of their interactome with regards to dysregulated phosphoproteins. 370 out of the 432 proteins are annotated as phosphoproteins, from which we identify a set of 37 phosphoproteins dysregulated in NF1 (3 proteins), NS (33 proteins), NSML (19 proteins), NSL (1 protein), NSLH (3 proteins), LS (2 proteins), CS (2 proteins), CFC (3 proteins) and JMML (2 proteins), some of them found in more than one syndrome. For most RASopathies, phosphorylation of downstream effectors was upregulated compared to the control, including but not limited to AKT, MEK1 and ERK1/ERK2. All RASopathies are associated to at least one phosphoprotein, in which the NS leads the number of identified phosphoproteins with 13 proteins, followed by NFNS (5 proteins) and JMML (5 proteins). Although just a few phosphoproteomic assays have been done in the onset of RASopathies, we encourage to outline future strategies to unveil the molecular events underlying the RASopathies using mass spectrometry approaches. This is an exciting field to widen molecular knowledge in rare diseases and even aim for translation into clinics.

Generation of a RASopathy proteins interactome
STRING (https:// string-db. org/) is a database of known and predicted protein-protein interactions. The interactions include direct (physical) and indirect (functional) associations; they stem from computational prediction, from knowledge transfer between organisms, and from interactions aggregated from other (primary) databases. In order to analyze the RASopathy proteins using STRING, we first compiled an updated panel of RASopathy The UniProt IDs for all the above 32 proteins were submitted to STRING database under the multiple proteins search. The following basic settings were used to generate a portable network graphic: confidence for the meaning of network edges (line thickness indicates the strength of data support); active interaction sources are based on experiments and databases only; and high confidence (0.700) as the minimum required interactome score. Edges represent protein-protein associations. For the generation of the RASopathy proteins interactome ( Fig. 2A) in STRING "none/query proteins only" and "none" values were used for the 1st and 2nd shell, respectively. These associations are meant to be specific and meaningful (i.e. proteins jointly contribute to a shared function), although this does not necessarily mean they are physically binding each other.

Generation of a RASopathy proteins and neighbors interactome
An interactome was generated including the above 32 RASopathy proteins and their neighbors. We first downloaded all protein-protein associations for each individual RASopathy protein from STRING (using export your current network protein annotations link) and complemented with those annotated in UniProt database (under the Interaction section, which provides information on interaction(s) with other proteins or protein complexes). Data compilation yielded a total of 765 proteins (Additional file 1: Tab1), including 511 unique proteins based on UniProt IDs (Additional file 1: Tab2). Functional annotation of the 511 unique proteins was carried out using DAVID (https:// david. ncifc rf. gov/).

Analysis of the main hub proteins of the interactome network
In order to further analyze relationships within the whole interactome derived from the 32 RASopathy proteins, we used STRING to build a network of the interactions with the above list of 511 unique proteins (Additional file 1). The resulting network has a total of 8191 edges representing protein-protein associations. Then, in order to find highly connected nodes within the network, the number of edges for every single node was calculated. We focused our analysis on selected nodes that had a number of interactions greater than 4 times the average of interactions found in the network (mean = 18.4 interactions).

Analysis of dysregulated phosphoproteins present in the RASopathy interactome
Differential phosphorylation levels were analyzed in the RASopathy interactome. UniProt IDs from the interactome (Additional file 1) and from the dataset of dysregulated phosphosites in RASopathies (Additional file 4) were crossed. The identified proteins were selected for further analysis. Phosphorylation levels are normalized to the total abundance of each phosphoprotein, providing harmonization of the data between experimental models and organisms. Fold changes in the RASopathy versus control were used to represent dysregulation in each phosphosite, transformed to base 2 logarithm (log2[fold change]). A threshold of 1.5 both in upregulation and downregulation was applied to filter out noise in the data. Qualitative data was processed in parallel assigning upand down-regulation based on original data. R 3.63 [142] and RStudio 1.2.5042 [143] were used to generate a heatmap representing fold changes of each phosphosite (R package pheatmap [144]). Clustering analysis was done in NS and NSML by selecting manually from quantitative data dysregulated phosphoproteins with available information for only NS syndrome, phosphoproteins equally dysregulated in both syndromes and phosphoproteins dissimilarly dysregulated in both syndromes.
The analysis of dysregulated phosphoproteins present in the RASopathy interactome presented some limitations. Different models, even when they mimic the same syndrome, may show a different response in terms of phosphorylation changes that may be driven in part, by the different experimental conditions. Therefore, experimental validation is recommended for those less studied phosphoproteins. In terms of statistical validation, some original phosphoproteomic studies [89][90][91][92] were done for screening purposes. Most immunoblots data was statistically validated including several biological replicates. These aspects are included in Additional file 4.

Analysis of RASopathy proteins phosphosites
We investigated the presence of phosphosites in the 27 confirmed RASopathy proteins using available information in the Phosphosite [141] database. All phosphosites up to available referenced works demonstrating phosphorylation of amino acids in the human sequence were included, independently of the study on their dysregulation in RASopathies (Additional file 6). Classification of the references into high-throughput and low-throughput studies was annotated, but both equally considered. Total number of phosphosites per protein with correlation to the cognate RASopathy syndrome, were represented in a bubble plot using R 3.63 [142] and RStudio 1.2.5042 [143] (R package ggplot2).

Gene ontology analysis and biological processes
The 27 RASopathy proteins (Fig. 3B) and the 33 dysregulated phosphoproteins in the RASopathy interactome (Fig. 4B) were annotated for GO terms using PANTHER (http:// www. panth erdb. org/) classification system tool. A statistical over-representation test using the human PANTHER GO-Slim Biological Process database was used. All groups were selected attending to its statistical significance (FDR value < 0.01) and according to the most specific subclass on the GO hierarchy term they belong to.

Limitations of the study for phosphoproteomics data
In its conceptualization, our study compiles information from articles including different models and statistical validations assuming some limitations. MS phosphoproteomics studies were used for screening purposes [89][90][91][92], and statistics are limited since no p-value test has been performed. Qualitative data from immunoblots is validated by replication of the experiment and fullfill a p-value threshold of 0.05. Details for each study can be consulted in Additional file 4.Therefore, we can only suggest herein the impact of phosphorylation changes in RASopathies and key phosphoproteins from a point of view of screening and based in the co-occurrence and relevance considering crossed data.