This article is mentioned in:

Introduction

Primary liver cancer is the sixth most prevalent type of cancer worldwide and the second leading cause of cancer-associated mortality (1). The incidence of primary liver cancer has increased globally, especially in developing regions. Of the 782,000 new cases diagnosed in 2012, China accounted for 50% (1). Mortality due to primary liver cancer has also increased, causing 745,000 deaths worldwide in 2012 (2). The majority of diagnosed primary liver cancer is classified as hepatocellular carcinoma (HCC), which accounts for 85–90% of all primary liver malignancies (3). Although the treatments for HCC have greatly improved patient survival, the overall 5-year survival rate remains poor (5%) due to the difficulty in early diagnosis (4). Currently, the only curative therapy for HCC is surgical resection, and the survival rate can reach 70% if the tumor diameter is <2 cm (4). However, even after hepatic resection, the recurrence rate of HCC is 70% at 5 years, which limits the long-term survival of patients with HCC (5). The most common methods to diagnose HCC are imaging tests and biopsy; however, these methods have a relatively high false positive or negative rate, respectively (6). α-fetoprotein is also used as a biomarker for HCC, but its sensitivity and specificity are limited (7). Due to this lack of effective diagnostic methods, only 10–20% of patients with HCC are diagnosed at the early stage (8). Therefore, identifying reliable predictive markers involved in tumorigenesis and understanding the molecular mechanism of HCC will contribute to early diagnosis and delay tumor progression.

Aberrant gene and protein expression accounts for various human diseases, such as cancer, bacterial infections and neurodegenerative disease (9–11). In addition to genome and proteome profiles, post-translational modifications (PTMs) are important in various cellular processes, including regulation of chromatin remodeling, cell cycle and cellular signaling events (12,13) associated with human health. Phosphorylation at serine, threonine and tyrosine residues is the most widespread mechanism for regulating protein activity, with ~50% of human proteins undergoing phosphorylation during their life cycle (14). Protein phosphorylation functions as a regulator of protein function, cellular localization and biological processes. Phosphorylation and dephosphorylation of proteins are involved in key biological mechanisms that either switch on or off enzyme activity (15). Aberrant protein phosphorylation is associated with numerous diseases, including tumors such as HCC (16). Therefore, it is important to develop global and quantitative methods for elucidating phosphorylation alteration in diseases.

Currently, label-free and stable isotope labeling approaches combined with liquid chromatography tandem-mass spectrometry (LC-MS/MS) are two mainstream quantification methods used to quantify cellular proteins and peptides (17). Phosphoproteomics research may have a potential relevance in HCC treatment and diagnostics, since it could be used to identify both potential biomarkers and drug targets for cancer. Elevated phosphorylation levels of 4E-binding protein 1 on Thr46 were reported as a potential biomarker for early recurrence and metastasis of HCC (18). Phosphorylation of branched-chain α-ketoacid dehydrogenase kinase on Ser31 can affect the activation of the ERK1/2 signaling pathway mediated by aminopeptidase N, which may be helpful to identify novel biomarkers and therapeutic targets for HCC (19). Proviral integration site of murine leukemia virus 3 promotes HCC migration by regulating the phosphorylation of multiple Rho GTPases that modulate RhoA (20). CeO2 nanoparticles (NPs) were reported as a potential therapeutic agent for HCC (21). Phosphoproteomic analysis revealed that CeO2 NP treatment affected the phosphorylation of proteins mainly associated with cell adhesion and RNA splicing (21). Additionally, sorafenib treatment can affect ferroptosis by altering protein phosphorylation in HCC (22).

Motif analysis of phosphorylated proteins can help identify potential kinases. Every phosphorylation site can be the substrate of one or more specific kinases (23). Casein kinase 2 (CK2), CDK1/2, protein kinase A (PKA), PKC, MAPK and ERK1/2 are some of the important kinases associated with cancer (24–37). CK2 kinases are a family of serine/threonine kinases that are involved in the regulation of diverse biological processes, such as cell proliferation, apoptosis and oncogenesis (24–27). CDKs are a family of serine/threonine kinases that are pivotal in cell cycle progression and transcription (28). Interference with CDK1 and CDK2 can affect chromosomal stability, as well as S and G2/M phase control, thereby leading to oncogenic events (29). Both PKA and PKC are serine/threonine kinases belonging to the family of AGC kinases, which are critical regulators of numerous cellular processes, including glucose metabolism, lipid metabolism, stress responses and apoptosis (30–32). MAPKs are proline-directed serine/threonine kinases, which are components of highly conserved signal transduction pathways (33). MAPKs are involved in diverse cellular processes, such as cell proliferation, differentiation, migration and apoptosis (33). There are four major groups of MAPK signaling pathways that have been characterized in mammalian cells, including ERK1/2, ERK5, Jun N-terminal kinase 1/2/3 and the p38 isoforms α/β/γ(ERK6)/δ (34,35). ERK1/2 kinases are a family of serine/threonine kinases that are the main effectors of the Ras-Raf-MEK-ERK signal transduction cascade, which is involved in a number of processes, including cell cycle progression, proliferation and oncogenesis (36,37).

In the present study, the differences in the phosphoproteome between HCC and normal human liver tissues were profiled using LC-MS/MS in combination with tandem mass tag (TMT). Phosphosites with differential abundance were identified, and motif analysis was applied to identify potential kinases. Functional enrichment analysis was performed to reveal the biological functions and cellular localization of phosphorylated proteins. Pathways and network analyses were applied to illustrate important changes in the components of regulated pathways in HCC. The present findings may provide a global map of the phosphoproteome of HCC, which may expand the current knowledge on the differences in protein phosphorylation between tumor and normal tissues of the liver, thus providing fundamental information for future studies.

Materials and methods

Protein sample preparation

Patients with HCC (n=4) were enrolled in the present study, and their tumor and adjacent non-tumor tissue samples were collected between May 2017 and July 2018 at the First Affiliated Hospital of Zhengzhou University (Zhengzhou, China). The distance between the adjacent non-tumor and tumor tissues was 1–5 cm. All patients were males and had hepatitis B virus infection. The mean age of the patients was 43 years old (SD, 2.8; age range, 41–46 years). The pathological grading of all 4 cases was moderate differentiated HCC according to the World Health Organization HCC grading classification (38). The diagnosis of HCC was based on the histological results of liver surgical resections. All patients provided written informed consent, and the study was approved by the Ethics Committee of Zhengzhou University (Zhengzhou, China).

To successfully extract proteins from liver tissues, the samples were first ground using liquid nitrogen and subsequently transferred to 5-ml centrifuge tubes. Subsequently, lysis buffer (8 M urea, 1% Protease Inhibitor Cocktail and 1% phosphatase inhibitor cocktail; MedChemExpress) was added at a ratio of 4:1 to the tissue. Subsequently, the lysate was sonicated for 3 sec with a 5 sec interval for three times using a high intensity Sonotrode set (Ningbo Scientz Biotechnology Co., Ltd.) at 4°C and then centrifuged at 12,000 × g for 10 min at 4°C. The supernatant was collected and the protein concentration was determined using a BCA protein assay kit (Beyotime Institute of Biotechnology).

Protein digestion and peptide labeling

Disulfide bonds of proteins were reduced using dithiothreitol (Sigma-Aldrich; Merck KGaA) at a final concentration of 5 mM for 30 min at 56°C. Proteins were alkylated using 11 mM iodoacetamide (Sigma-Aldrich; Merck KGaA) for 15 min at room temperature in the dark. Protein solutions were then diluted with 100 mM triethylammonium bicarbonate buffer (TEAB; Sigma-Aldrich; Merck KGaA) to an urea concentration <2 M. Trypsin was added at a 1:50 trypsin-to-protein mass ratio, and the sample was incubated overnight at 37°C. Subsequently, a second aliquot of trypsin was added at a 1:100 trypsin-to-protein mass ratio for 4 h at 37°C.

The digested peptides were desalted using a Strata X C18 SPE column (Phenomenex, Inc.) and were lyophilized in a vacuum centrifuge at 300 × g for 6 h at room temperature. Lyophilized peptides were reconstituted in 0.5 M TEAB buffer and were labeled using a TMT kit according to the manufacturer's protocol (cat. no. 90406; Thermo Fisher Scientific, Inc.). TMT reagents were dissolved in acetonitrile and then added to the sample. The labeling reaction was conducted for 2 h at room temperature. Subsequently, the labeled peptides were pooled and desalted, followed by lyophilization in a vacuum centrifuge, as aforementioned.

Phosphopeptide enrichment

The phosphopeptides were fractionated via high-pH reverse-phase high-performance liquid chromatography (HPLC; Agilent 1260 HPLC; Agilent Technologies, Inc.) using a Thermo BetaSil C18 column (5-µm particles; 10-mm ID; 250-mm length; Thermo Fisher Scientific, Inc.). The column temperature was 35°C. The detection wavelength was 214 nm. The sample quantity was 1 ml. The mobile phase A was water and acetonitrile (98:2, v/v). The mobile phase B was water and acetonitrile (2:98, v/v). The elution gradient was as follows: 0–5 min (5% mobile phase B); 5–10 min (5-8% mobile phase B); 10–67 min (8-32% mobile phase B); 67–69 min (32-95% mobile phase B); 69–80 min (95% mobile phase B); 80–85 min (95-5% mobile phase B). The flow rate was 1 ml/min. The fractionated peptides were first suspended in loading buffer (50% acetonitrile and 6% trifluoroacetic acid) and then incubated with pre-treated immobilized metal ion affinity chromatography (IMAC) microspheres (J&K Scientific Ltd.). Phosphopeptide-bound IMAC microspheres were washed once with the loading buffer, followed by three washes with the washing buffer (30% acetonitrile and 0.1% trifluoroacetic acid). Subsequently, the phosphopeptides were eluted using 10% ammonium hydroxide, desalted with C18 ZipTips (EMD Millipore) and dried by vacuum centrifugation for LC-MS/MS analysis.

LC-MS/MS analysis

The enriched phosphopeptides were analyzed using an EASY-nLC 1000 UPLC system coupled in line with an Orbitrap Fusion™ mass spectrometer (Thermo Fisher Scientific, Inc.) via a nanospray ionization source in positive mode. Each sample was loaded onto a reversed-phase analytical column (75-µm inner diameter; 15-cm length) packed in-house in solvent A (0.1% formic acid) and resolved at a 400-nl/min flow rate using the following gradient: 6–23% solvent B (0.1% formic acid in 98% acetonitrile) for 26 min, 23–35% solvent B for 8 min and 35–80% solvent B for 3 min, then holding at 80% solvent B for the last 3 min. The electrospray voltage was set to 2.0 kV. The nitrogen gas temperature was 320°C and the nebuliser pressure was 15 psi. MS was operated in the data-dependent acquisition mode. The full scan MS1 spectra range was 350–1,800 m/z. The resolution in the Orbitrap was 60,000. The MS2 scan was measured at a resolution of 30,000 at an m/z of 100. The 20 most intense precursor ions were selected and fragmented by higher-energy collision-induced dissociation with a normalized collision energy at 35%. Automatic gain control was set at 5×104 to determine the target value. The signal threshold was set to 5,000 ions/sec and the maximum ion injection time was set to 200 msec. Dynamic exclusion was set at 15 sec to avoid the repeated detection of the same fragment ion peaks. The mass spectrometry proteomics data were deposited in the ProteomeXchange Consortium via the Proteomics Identification Database (39) (dataset no. PXD013934). For the quantitative analysis, statistical significance was determined using unpaired Student's t-tests. P<0.05 and a differentially expressed ratio >1.5 or <0.67 was considered to indicate a statistically significant difference.

Database search

The LC-MS/MS data were analyzed using Maxquant v.1.5.2.8 (40). The data were searched against the SwissProt Human database (https://www.uniprot.org/taxonomy/?query=Homo+sapiens+%28Human%29+&sort=score) (20,130 entries) concatenated with its sequence-reversed database. Contaminants frequently observed were also added. Trypsin/P was the proteolytic enzyme, and two missed cleavages were allowed. Carbamidomethylation of cysteine was defined as a fixed modification, and oxidation of methionine, acetylation of N-terminal protein and phosphorylation of serine, threonine and tyrosine were specified as variable modifications. The minimum peptide length was set to 7 amino acid residues, and the maximum number of modifications per peptide was set to 5. Quantification was performed using TMT 10plex (Thermo Fisher Scientific, Inc.). The mass tolerance for precursor ions was set to 20 ppm in the first search and to 5 ppm in the main search, while for fragment ions, the mass tolerance was set to 0.02 Da. Regarding the level of proteins and peptide spectrum matches, the maximum false discovery rate was adjusted to 1%.

Bioinformatic analysis

The UniProt-Gene Ontology annotation (GOA) database (http://www.ebi.ac.uk/GOA/) was used to detect the significantly enriched Gene Ontology (GO) terms. Since some identified phosphoproteins were not annotated on the UniProt-GOA database, the InterProScan software (version 5.47-82.0; http://interproscan-docs.readthedocs.io/en/latest/Introduction.html) based on the protein sequence alignment method was used. Subsequently, the identified phosphoproteins were automatically assigned to three categories of biological processes, cellular components and molecular functions. Wolfpsort software (version 0.1; http://wolfpsort.hgc.jp/) was used to predict subcellular localization. The InterPro database (https://www.ebi.ac.uk/interpro/) was used to analyze the functional domains of differentially expressed proteins. Pathway enrichment analysis was conducted using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https://www.kegg.jp/). A two-tailed Fisher's exact test was used to filter the enrichment result to a corrected P<0.05. The protein-protein interaction (PPI) network was obtained using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING; http://www.string-db.org/) database. To determine the significance of the differentially expressed phosphoproteins against all the identified phosphoproteins, a two-tailed Fisher's exact test was used in each category, and a corrected P<0.05 was considered to indicate a statistically significant difference. Motifs of the phosphorylated sites for serine, threonine and tyrosine residues were analyzed using Motif-X (41). Peptide sequences were centered on each phosphorylation site and extended to 10 amino acids upstream and downstream of the site. The minimum occurrence was set to 20, and the significance was set to 0.000001. The new phosphorylation sites found in the present study were identified by comparing the dataset in the current study with the phosphoSitesPlus database (https://www.phosphosite.org/).

Western blotting

Tumor and adjacent non-tumor tissue samples collected from 6 different patients with HCC were used to further verify the bioinformatic results using western blotting. All patients were males and their mean age was 53 years (SD, 7.5; age range, 46–67 years). The tissues were also collected between May 2017 and July 2018 at the First Affiliated Hospital of Zhengzhou University. All patients provided written informed consent. Proteins were extracted from tissues on ice using RIPA buffer (Beyotime Institute of Biotechnology) supplemented with a protease inhibitor cocktail and phosphatase inhibitors for 30 min at 4°C. Protein concentration was determined using the BCA method. A total of 20 µg/lane protein was separated via 10% SDS-PAGE and transferred onto polyvinylidene difluoride membranes. The membranes were blocked at room temperature for 2 h with 5% BSA (Beyotime Institute of Biotechnology), followed by incubation overnight at 4°C with primary monoclonal antibodies against minichromosome maintenance complex component 2 (MCM2; 1:1,000; cat. no. 10513-1-AP; ProteinTech Group, Inc.), phosphorylated (p-)MCM2 (1:10,000; cat. no. ab109270; Abcam), vimentin (1:1,000; cat. no. 10366-1-AP; ProteinTech Group, Inc.), p-vimentin (1:1,000; cat. no. 13614; Cell Signaling Technology, Inc.) and GAPDH (1:10,000; cat. no. 60004-1-Ig; ProteinTech Group, Inc.), with the latter used as the loading control. After washing 3 times with TBS-Tween (0.1% Tween-20) for 5 min at 37°C, the membranes were incubated with HRP-conjugated Affinipure goat anti-mouse or anti-rabbit secondary antibodies for 2 h at room temperature (1:10,000; cat. nos. SA00001-1 or SA00001-2, respectively; ProteinTech Group, Inc.). Proteins were visualized using enhanced chemiluminescence detection reagents (Beyotime Institute of Biotechnology). The software used for densitometry was Quantity One 4.6.6 (Bio-Rad Laboratories, Inc.).

Statistical analysis

GraphPad Prism 5.0 (GraphPad Software, Inc.) was used for statistical analysis. The results of the experiments were expressed as the mean ± SD (n=6). Paired Student's t-test analysis was used to determine the statistically significant differences. P<0.05 was considered to indicate a statistically significant difference.

Results

Phosphoproteome profiles of HCC and normal tissues

To identify the phosphoproteomic changes that result from HCC, four pairs of tumor and adjacent non-tumor tissues were resected from patients diagnosed with HCC. These tissues were frozen immediately after resection and stored in liquid nitrogen until lysis to maintain potentially unstable phosphorylation status. The proteins from HCC and normal tissues were extracted and then digested using trypsin. The tryptic peptides were labeled with TMT isobaric mass tags for multiple and precise quantification of phosphorylation levels, and then mixed together in equal quantities, followed by fractionation of the peptides using high-pH HPLC. The phosphopeptides were subsequently enriched by IMAC, and analyzed via LC-MS/MS. The distribution of mass errors was ~0, and all were <10 ppm (Fig. S1A). The majority of the identified peptides were between 8 and 20 amino acids in length, which is in accordance with the pattern of trypsin-digested peptides (Fig. S1B). These results indicated that the sample preparation method met the standards of LC-MS/MS for quantitative proteomics and that the MS data were accurate. The LC-MS/MS quantitative method then determined whether the peptides present in the samples were phosphorylated (Fig. 1), with MCM2, ADH4 and TP53BP1 used as examples.

Figure 1. Tandem-mass spectrometry spectra for MCM2, ADH4 and TP53BP1. The relative intensities of the tandem mass tag reporter ions in the different tissue samples (A1, C1, D1 and F1 represent tumor tissues, while A2, C2, D2 and F2 represent normal tissues) are shown on the upper right side MCM2, minichromosome maintenance complex component 2; ADH4, alcohol dehydrogenase 4; TP53BP1, TP53-binding protein 1.

The phosphoproteomic analysis of HCC and normal liver tissues resulted in the identification of 6,129 total unique phosphorylation sites from 2,687 proteins, which were predicted via Maxquant. A total of 4,780 phosphorylation sites from 2,209 proteins were quantified and were considered for further analysis (Table SI). Among these, single, double, triple and quadruple or higher phosphorylated peptides represented 55, 19, 11 and 15% of the total phosphopeptides, respectively (Fig. 2A). The observed phosphorylation sites were divided into four classes according to their phosphorylation site localization probability score calculated using Maxquant: 6,129 phosphorylation sites were classified as class I with a specific residue with high confidence (≥0.75); 1,008 were class II sites (0.5-0.75); 989 were class III sites (0.25-0.5); and 24 were class IV sites (<0.25). In total, 75% of all phosphorylation sites were classified as class I sites and were considered to be accurately quantified and localized, and were therefore further analyzed (Fig. 2B). The non-redundant phosphorylation sites in class I comprised 5,446 phosphorylated serine residues (89%), 655 phosphorylated threonine residues (11%) and 28 phosphorylated tyrosine residues (<1%) (Fig. 2B). Comparing the dataset in the current study with the phosphoSitesPlus database, 317 phosphoproteins and 138 phosphorylation sites were identified for the first time in the present study (Fig. 2C). For the quantitative analysis, statistical significance was determined using unpaired Student's t-tests. P<0.05 and a differentially expressed ratio >1.5 or <0.67 was considered to indicate a statistically significant difference. The difference between the HCC and normal tissue groups of each phosphosite was plotted with the corresponding P-value, as shown in Fig. 2D. A total of 533 differentially expressed phosphosites from 376 proteins were identified in HCC, specifically 74 phosphosites that were upregulated and 459 phosphosites that were downregulated compared with in normal tissues (Table SII).

Figure 2. Identification and characterization of HCC phosphoproteome. (A) Distribution of phosphorylation sites on each phosphopeptide. (B) Distribution of phosphorylation sites localization probabilities, and the proportion of phosphorylated amino acid residues on phosphopeptides of class I sites. (C) Overlap of phosphoproteins and phosphorylation sites identified in the present study with those in the phosphoSitesPlus database. (D) Volcano plot of significantly dysregulated phosphosites in HCC. S, serine; T, threonine; Y, tyrosine; HCC, hepatocellular carcinoma.

Motif analysis

Considering that tyrosine phosphorylation comprised <1% of all the localized phosphorylation sites observed, only phospho-serine motifs and phospho-threonine motifs were analyzed using Motif-X. The motif enrichment heat map of all the upregulated and downregulated phosphoproteins is shown in Fig. S2. In total, 53 phospho-serine motifs and 8 phospho-threonine motifs were identified (Table SIII). The majority of phospho-serine motifs were categorized into three classes according to the literature (42): Acidic, proline-directed and basic motifs. Proline-directed and acidic motifs accounted for 15 and 22 of the 53 identified phospho-serine motifs, respectively, while 12 of the 53 phospho-serine motifs belonged to the basic motifs class (Table SIII). The proline-directed motifs were predicted to be associated with MAPK, CDK, ERK1/2, PKA and PKC (Fig. 3A). The acidic motifs were predicted to be associated with CK2 (Fig. 3B), while the basic motifs were predicted to reflect activation of PKA and PKC (Fig. 3C).

Figure 3. Motif analysis of hepatocellular carcinoma phosphoproteome. The phosphorylated residue is located at the central position, within a 13-mer phosphorylated peptide sequence. Motifs were classified according to the chemical properties of the peptide sequences as (A) proline-directed, (B) acidic and (C) basic. ‘N’ represents the number of times that a motif was observed. PKA/C, protein kinase A/C; CK2, casein kinase 2.

GO annotation classification

All the identified proteins were subjected to GO annotation in order to assess differences in biological processes, molecular functions and cellular components to explain the biological effects of phosphoproteins. The phosphoproteins identified via LC-MS/MS were involved in various biological processes (Fig. 4A), including ‘cellular process’ (14%), ‘single-organism process’ (12%), ‘biological regulation’ (12%) and ‘metabolic process’ (9%). The identified phosphoproteome in the differentially expressed phosphoproteins was also classified according to molecular functions (Fig. 4B), revealing that most phosphoproteins were classified into the categories of ‘binding’ (56%), ‘catalytic activity’ (20%), ‘structural molecule activity’ (7%) and ‘molecular function regulator’ (7%). The cellular components of phosphoproteins were determined as ‘cell’ (24%), followed by ‘organelle’ (22%), membrane-enclosed lumen (13%), macromolecular complex (12%) and ‘membrane’ (11%) (Fig. 4C). Subcellular localizations of differentially expressed phosphoproteins were predicted and classified using wolfpsort software. The sub-cellular protein annotation for the differentially expressed phosphoproteins revealed that the majority of phosphoproteins were localized to the nucleus (53%), followed by cytoplasm (33%), mitochondria (9%) and plasma membrane (3%) (Fig. 4D).

Figure 4. Gene Ontology analysis of phosphoproteins. The differentially regulated phosphorylation sites (HCC versus normal tissues) were classified according to (A) biological processes, (B) molecular functions, (C) cellular components and (D) subcellular localization.

Functional enrichment analysis

Functional enrichment analysis revealed that the target phosphoproteins were enriched in numerous processes. GO enrichment analysis of molecular functions revealed that the phosphorylated proteins were significantly enriched in processes involved in binding activity and structural molecule activity. GO enrichment analysis based on cellular components illustrated that the phosphorylated proteins were enriched in various cellular components, such as ‘anchoring junction’, ‘extracellular organelle’ and ‘extracellular exosome’. GO enrichment analysis of biological processes demonstrated that the phosphorylated proteins were mainly enriched in biosynthetic and metabolic processes (Fig. 5A and B). KEGG enrichment analysis revealed that the upregulated targets were enriched in eight pathways (Fig. 5C), while the downregulated targets were enriched in 11 pathways (Fig. 5D). The upregulated phosphorylated proteins were most significantly enriched in the processes of ‘RNA transport’, ‘DNA replication’ and ‘tyrosine metabolism’, while the downregulated phosphorylated proteins were most significantly enriched in the processes of ‘arrhythmogenic right ventricular cardiomyopathy’, ‘ABC transporters’ and ‘glycolysis/gluconeogenesis’. These pathways were involved in genetic information processing, metabolism, environmental information processing, cellular processes and human diseases. Additionally, a domain enrichment analysis was performed. The identified proteins were annotated using the InterProScan software and the InterPro domain database. There were 9 domains that matched the upregulated proteins, such as ‘alcohol dehydrogenase, N-terminal’, ‘Gas2-related domain’ and ‘NAD(P)-binding domain’ (Fig. 5E). Additionally, there were 25 domains that matched the downregulated proteins, such as ‘ABC transporter type 1, transmembrane domain’, ‘calponin homology domain’ and ‘riboflavin synthase-like beta-barrel’ (Fig. 5F).

Figure 5. Functional enrichment analysis of phosphoproteins. Gene Ontology enrichment analysis of (A) upregulated and (B) downregulated phosphoproteins in terms of molecular function, cellular component and biological process. Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis of (C) upregulated and (D) downregulated phosphoproteins. Protein domain enrichment analysis of (E) upregulated and (F) downregulated phosphoproteins.

Prediction of PPI network

To further understand the cellular processes regulated via phosphorylation in HCC, a PPI network of the phosphoproteins was established (Fig. 6A). All the upregulated and downregulated phosphoproteins were uploaded to STRING, presenting a global view of the networks of phosphoproteins in HCC. It is well known that protein phosphorylation is closely associated with cell signal transduction. Therefore, STRING was used to construct a PPI network to illustrate the association between the identified phosphoproteins and signal transduction pathways in HCC (Fig. 6B). The PPI network was composed of 173 phosphoproteins participating in the hypoxia-inducible factor-1 (HIF-1), mTOR and MAPK signaling pathways, as well as in cell cycle and apoptosis, which may be involved in networks associated with cancer functions, directly or indirectly. There were 4 upregulated and 22 downregulated phosphoproteins, while LMNB2 included both upregulated and downregulated phosphorylation sites (Fig. 6B). The network indicated phosphoprotein dysregulation in HCC, thus providing information on the initiation and progression of HCC, and revealing potential targets for therapeutic intervention.

Figure 6. PPI network. (A) PPI network for the upregulated and downregulated phosphoproteins in hepatocellular carcinoma. (B) PPI network of phosphoproteins associated with the HIF-1, mTOR and MAPK signaling pathways, cell cycle and apoptosis. PPI, protein-protein interaction; HIF-1, hypoxia-inducible factor-1.

Validation of phosphorylated proteins via western blotting

To verify the protein phosphorylation data obtained via LC-MS/MS, western blotting was performed. As shown in Fig. 7, MCM2 exhibited an increased phosphorylation level in HCC, while vimentin exhibited a decreased phosphorylation level in HCC compared with in normal tissues. Therefore, western blotting confirmed the reliability of the LC-MS/MS data.

Figure 7. Western blot validation. Differential phosphorylation levels of vimentin and MCM2 in HCC and normal tissues were verified using western blot analysis. *P<0.05. HCC, hepatocellular carcinoma; MCM2, minichromosome maintenance complex component 2; p, phosphorylated.

Discussion

In the present study, a high-resolution mass spectrometry-based quantitative phosphoproteomics approach was conducted to analyze the phosphoproteome dysregulation in HCC. The IMAC method was used to enrich the phosphorylated peptides, which were then analyzed via LC-MS/MS, so that the identification of phosphorylated peptides was not affected by the abundance of the protein. A total of 6,129 unique phosphorylation sites were identified. Among these phosphorylation sites, 4,780 were quantified, thus allowing to calculate the site-specific differences of phosphoprotein expression between HCC and normal tissues. The phosphorylation sites with differential abundance may help to reveal the possible roles that phosphoproteins serve in cancer development and progression. The distribution of phosphorylation sites indicated that those with high occupancy tended to be associated with serine, while those with low occupancy were likely to be associated with threonine and tyrosine, which was in accordance with previous studies (43,44). Focusing on phosphorylated proteins that were 1.5-fold upregulated and 0.67-fold downregulated, 533 dysregulated phosphorylation sites were identified. Since protein kinases phosphorylate their substrates at specific motifs, motif analysis was applied to clarify the presence of activated protein kinases.

The functions of the phosphorylation sites with differential abundance were classified based on GO annotation. The functional annotation of the upregulated and downregulated proteins was mainly based on biological processes, molecular functions and cellular components. According to the literature (45), phosphorylation is often involved in the temporary regulation of protein function, which can alter protein conformation, target proteins for destruction, promote or interrupt protein interactions, and activate or deactivate enzymatic activity. In the present study, the largest group of phosphorylated proteins consisted of proteins with binding activity, accounting for 56% of all the identified phosphorylated proteins, which indicated that phosphorylation may serve critical roles in regulating PPIs. The second most abundant protein group was associated with catalytic activity (20%). In addition, >146 phosphoproteins were identified to be associated with signal transduction activity, ranking fourth among all categories. Since multiple phosphoproteins may serve a critical role in signal transduction pathways, alternations in phosphorylation patterns may explain the differences in oncogenic properties between normal and tumor tissue. The sub-cellular localization of proteins is important for their functions. Previous studies revealed that the nuclear organization of cancer cells is modified (46–49). During cancer development, several mechanisms associated with the control of gene expression are modified. The concurrent regulation of phosphorylation events responding to modifications in the nuclear environment of cancer cells suggests that phosphorylated proteins may be associated with the transcriptional regulation of tumor-associated genes (50). Sub-cellular localizations of HCC phosphoproteins implied a high probability that the mechanisms of hepatic oncogenesis are associated with the phosphoproteins identified in the cell nucleus. Mitochondria, which serve critical roles in cell metabolism, survival and apoptosis (23), may also be important in hepatic oncogenesis, as 29 dysregulated phosphoproteins were associated with this organelle. As the volume of mitochondria is small, the proportion of phosphoproteins in mitochondria was relatively high. Furthermore, mitochondrial proteins can either act as protein kinases and protein tyrosine phosphatases, or as kinase substrates (51). Therefore, phosphorylation serves an important role in the biological processes of mitochondria. The GO annotation results revealed that the phosphoproteins identified in the present study were associated with several aspects of HCC development and progression, including metabolism, cell proliferation, migration, metastasis and transcriptional control of gene expression.

GO enrichment analysis further demonstrated that biological processes, molecular functions and cellular components may be regulated by protein phosphorylation. Phosphorylation appears to be a key regulatory mechanism in HCC. KEGG pathway enrichment analysis was used to deepen the understanding of molecular mechanisms and biological pathways associated with HCC. KEGG-based functional enrichment analysis revealed that the upregulated phosphorylated proteins were most significantly enriched in the process of ‘RNA transport’, ‘DNA replication’ and ‘tyrosine metabolism’, while the downregulated phosphorylated proteins were most significantly enriched in the processes of ‘arrhythmogenic right ventricular cardiomyopathy’, ‘ABC transporters’ and ‘glycolysis/gluconeogenesis’. It is well known that major features of tumors are due to errors in RNA transport and DNA replication, which can lead to cancer and metabolic abnormalities. Abnormal tumor metabolism mainly includes glucose, lipid and amino acid metabolism (52–55). Additionally, KEGG enrichment analysis of the phosphorylated proteins revealed pathways associated with arrhythmogenic right ventricular cardiomyopathy. This may be explained by the elevated levels of aldosterone and neutrophil gelatinase-associated lipocalin that are induced by cancer cachexia (56).

Since protein phosphorylation can regulate small protein motifs associated with protein interactions, a domain enrichment analysis was conducted in present study. Upregulated phosphorylated proteins were enriched in 9 protein domains, including ‘Aldolase-type TIM barrel domain’, ‘Gas2-related domain’ and ‘NAD(P)-binding domain’. Phosphorylation of proteins containing Aldolase-type TIM barrel domains may serve an important role in glucose metabolism and oxidative stress (57,58). The Gas2-related domain is critical in cell cycle, cell abscission after division and cytoskeletal function in tissues (59). Proteins with an NAD(P)-binding domain are important coenzymes involved in catalyzing redox or non-redox reactions (60). The results suggested that phosphorylation has important effects on a variety of processes, including stress response and biosynthesis.

Pathway analysis is commonly used to identify known biological pathways or biological processes that are potentially involved in differentially expressed proteins associated with certain diseases (61). The information derived from this analysis can improve the understanding of biological processes behind diseases. As pathway analysis is based on known data, it cannot be used to identify new biological pathways or biological processes. PPIs occur widely in living organisms and serve important roles in regulating cellular processes (62). To further understand the effects of interactions between phosphoproteins on the regulation of HCC development and progression, a PPI network of the upregulated and downregulated proteins was established using STRING. The interactions among these phosphoproteins are likely to contribute to their cooperation and coordination in HCC. As apoptotic disorders and uncontrolled cell proliferation contribute to tumor malignancy, phosphoproteins involved in the cell cycle and apoptosis, as well as phosphoproteins participating in associated signaling pathways, such as the MAPK, mTOR and HIF-1 signaling pathways, were uploaded onto STRING to predict a PPI network. The results suggested that the PPIs among these phosphoproteins were involved in multiple signaling pathways, such as the MAPK, mTOR and HIF-1 signaling pathways, in HCC. As protein phosphorylation is important for the activities of signal effectors, signaling regulation and cellular processes (63), further analysis of the identified phosphoproteins may help to identify important regulatory factors.

The present study has some limitations. Firstly, the sample size was small and therefore the current findings should be validated in future large scale studies. Additionally, there is a lack of TNM staging information, which may be used to further understand the association between protein phosphorylation and HCC progression and prognosis. Despite these limitations, the results of the present study may contribute to improve the current knowledge on the phosphoproteome of HCC and provide scientific evidence that phosphorylation of proteins may be closely associated with the development and progression of HCC.

In conclusion, the present study used a TMT approach combined with LC-MS/MS to characterize the global phosphoproteome profile alteration in HCC tissues. The phosphosites identified both in HCC and normal liver tissues were compared, and 4,780 phosphosites were identified and quantified. Dysregulation of phosphoproteins was associated with various biological processes and molecular functions, which may be important in HCC. Identifying differentially expressed phosphoproteins and analyzing their functions, associated pathways and protein networks may contribute to improve the understanding of the development and progression of HCC.

Supplementary Material

Supporting Data

Supporting Data

Supporting Data

Supporting Data

Acknowledgements

Not applicable.

Funding

The present study was supported by the National Science and Technology Major Project of China (grant no. 2018ZX10302205), the Natural Science Foundation of Henan Province (grant no. 182300410361) and the Major Project of Science and Technology in Henan Province (grant no. 161100311400).

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Authors' contributions

JYZ and JTZ contributed to the fundraising. JYZ, JTZ and YL contributed to the experimental design. JTZ and YL drafted the manuscript. YL, QZ, FX, YZ and HC performed the experiments. KW, WW and WH analyzed the data. All authors read and approved the final version of the manuscript.

Ethics approval and consent to participate

All patients provided written informed consent, and the study was approved by the Ethics Committee of Zhengzhou University (Zhengzhou, China).

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References