Prediction of Biomarkers for Hepatocellular Carcinoma Based on Proteomics and Phosphoproteomics

Introduction

Primary liver cancer ranks among the most prevalent malignant tumors in globally, standing as the second leading cause of cancer-related deaths. HCC is the main type of primary liver cancer, accounting for 75%–85% of cases globally.¹ As a widespread disease, HCC exhibits high incidence and mortality rates. Unfortunately, patients usually have no obvious symptoms in the early stages, typically being diagnosed at advanced stages. Consequently, identifying biomarkers and potential therapeutic targets for early diagnosis and intervention becomes paramount in unraveling the disease’s internal biological connections and molecular characteristics.

In recent years, with the maturation of proteomics technology, research on tumor markers has increasingly shifted focus from genomics to proteomics. Research has shown that the production of differentially expressed proteins is either derived from abnormal products of tumor cells themselves or substances produced by the host during the anti-tumor process.² Research pertaining to HCC has revealed that proteins can influence disease states through various pathological and molecular mechanisms, including metabolism, ferroptosis, cell apoptosis, inflammatory response, and immune regulation. This effect might be attributed to changes in protein activity following modifications such as phosphorylation, nitration, acetylation, and ubiquitination.³

Protein phosphorylation modification is the most common and crucial post-translational modification, playing a role in almost all life activities. Catalyzed by protein kinases, this modification entails the addition of a phosphate group to the amino acid side chains. Conversely, the dephosphorylation process is facilitated by phosphatases. Serine (Ser), threonine (Thr), and tyrosine (Tyr) are the most frequently phosphorylated amino acid residues.⁴ Serine phosphorylation primarily functions to activate conformational changes in proteins. Once initiated by phosphorylation, the primary signal typically leads to sequential phosphorylation of downstream proteins, forming a phosphorylation cascade. For instance, nuclear factor-κB (NF-κB) has emerged as a critical player in HCC. NF-κB-inducing kinase (NIK) phosphorylation and activates inhibitor of NF-κB kinase subunit alpha (IKKα), and the hepatic NIK/IKKα cascade suppresses liver regeneration.⁵ Additionally, studies have found that the phosphorylation of NF-κB p65 is significantly upregulated in inflammation-related HCC patients, and β-arrestin1 directly bounds p65 to promote the phosphorylation of NF-κB p65 at ser536, resulted in cell malignant proliferation through GSK3β/mTOR signalling. This underscores the significance of phosphorylation-mediated signaling regulation in cancer progression.⁶

To identify candidate biomarkers for the transition from LCT or ANT to HCT, we applied quantitative proteomics and TMT-labeled quantitative phosphoproteomics, systematically prioritized core differential proteins across both omics layers, and used KEGG to group core proteins into lipid, amino-acid, and carbohydrate metabolism modules. Integrative analysis of protein abundance and dysregulated phosphorylation establishes a robust framework for elucidating stage-specific molecular features and the modification networks that drive malignant transformation in HCC.

Materials and Methods

Patients

A total of 40 subjects were recruited from the Dalian Public Health Clinical Center between January 2016 and May 2021. All samples in this study were obtained with written informed consent. The research protocol was approved by the Ethics Committee of the Dalian Public Health Clinical Center (Approval ID: 2021–027-002). The collected tissue samples include HCT (n=20), ANT (>5cm) (n=20), and LCT (n=20). Liver cirrhosis (LC) and HCC were diagnosed according to the criteria of the Asian Pacific Association for the Study of the Liver. The patient’s basic information such as gender and age are listed in Table 1.

Table 1 The Characteristics of Study Population

Proteomic Mass Spectrometry Analysis

The following abbreviations are used: DTT, dithiothreitol; IAA, iodoacetamide; FA, formic acid; ACN, acetonitrile; DDA, data-dependent acquisition; DIA, data-independent acquisition.

Sample Preparation

Proteomic analysis was performed on HCT (n=20), ANT (n=20) and LCT (n=20) samples. Equal quantity of proteins for each sample were incorporated into SDT buffer (4% SDS, 100 mM DTT, 150 mM Tris-HCl pH 8.0). Protein concentrations were determined using the BCA Protein Assay Kit (Beyotime, P0012), and samples were then enzymatically hydrolyzed with trypsin. An equal aliquot from each sample in this experiment was pooled into one sample for construct a spectral library.⁷ DTT was added to achieve a final concentration of 20 mM, and the mixture was incubated at 37°C for 1 hour. Up cooling to room temperature, IAA (Bio-Rad, 163–2109) was added to a final concentration of 25 mM, and the solution was incubated in the dark for 30 minutes. An appropriate volume of NH4HCO3 buffer (50 mM) was added to dilute the UA concentration to less than 1.5 M. The protein suspensions were digested with 2 μg of trypsin in 40 μL of NH4HCO3 buffer overnight at 37°C, and the resulting peptides were collected as a filtrate. The peptides from each sample were desalted using C18 Cartridges (Sigma, 66872-U), concentrated by vacuum centrifugation and reconstituted in 0.1% FA. Peptide content was estimated by measuring the UV light absorbance at 280 nm. Take 100μg of high-abundance and low-abundance peptide segments after separation, using Thermo Scientific Pierce High pH Reversed-Phase Peptide Fractionation Kit (Thermo Scientific, 84868) for grading, and collect all components. After freeze-drying each peptide segment, it was dissolved in 10 μL of 0.1% FA, and the peptide concentration was determined using OD280.

LC-MS/MS Analysis and Data Processing

For DIA and DDA experiments, iRT (Biognosys) calibration peptides were spiked into 200 ng of each sample. All samples underwent chromatographic separation using the Evosep one system liquid chromatography (Denmark), followed by DDA mass spectrometry and DIA mass spectrometry using a timsTOF mass spectrometer (Bruker).⁸ Chromatographic separation was performed on an Evosep One system using the 30 samples per day (30SPD) method. The mobile phases were 0.1% FA in water (A) and 0.1% FA in ACN (B). DDA detection mode: positive ion, Analysis time: 44 min. We set the accumulation and ramp time was 100 ms each and recorded mass spectra in the range of 100–1700 m/z under positive electrospray mode. Additional parameters were as follows: dynamic exclusion: 24.0 s, Ion source Voltage: 1500 V, temperature: 180°C, dry gas: 3 L/min, TIMS range: 0.75–1.35, PASEF cycles: 8. DIA detection mode: positive ion. The mass spectrometer collected ion mobility MS spectra over a mass range of 100–1700 m/z, We defined up to eight windows for single 100 ms TIMS scans according to the m/z-ion mobility plane. During PASEF MSMS scanning, the collision energy was ramped linearly as a function of the mobility from 20eV at 1/K0=0.6 Vs/cm2 to 59 eV at 1/K0=1.6 Vs/cm2. For DDA data, we searched the FASTA sequence database (https://www.uniprot.org/) using SpectronautTM 14.4.200727.47784 software (Biognosys) and constructed a spectral library. DIA data were analyzed with the same software (SpectronautTM 14.4.200727.47784) by searching against the above-constructed spectral library.

Phosphoproteomic Mass Spectrometry Analysis

The following abbreviations are used: TMT, tandem mass tag; TiO₂, titanium dioxide; DHB, 2,5-dihydroxybenzoic acid; FA, formic acid; ACN, acetonitrile; HCD, higher-energy collisional dissociation.

Sample Preparation

Phosphoproteomic analysis was performed on HCT (n=3), ANT (n=3), and LCT (n=3) samples. Samples were ground into powder using liquid nitrogen, and proteins were extracted using SDT (4% (w/v) SDS, 100 mM Tris/HCl pH7.6, 0.1 M DTT) lysis method and quantified by BCA Protein Assay Kit (Thermo Scientific, 23225). Each sample was digested with trypsin using the Filter aided proteome preparation (FASP) method.⁷ The peptides were desalinated with C18 Cartridge (Sigma, 66872-U), lyophilized and then redissolved with 40μL dissolution buffer. Peptide concentrations were determined using OD280. A 100 μg aliquot of peptides from each sample was labeled according to the instructions of the TMT labeling kit (Thermo Scientific, A44520). After labeling, TiO₂ enrichment was performed. The labeled, mixed peptide solution was freeze-dried under vacuum and reconstituted with 1 × DHB buffer (prepared by diluting 5 × DHB buffer (3% DHB, 80% ACN, 0.1% TFA)). TiO2 beads were added to the solution and shaken for 40 minutes, then centrifuged, and the supernatant removed. The beads were transferred into a plug tip and washed three times with washing buffer 1, followed by three washes with washing buffer 2. After elution with elution buffer, the collected phosphorylated peptide segments were vacuum-concentrated, dissolved in 10 μL of 0.1% FA, and 5 μL of the solution was used for mass spectrometry analysis.⁹

LC-MS/MS Analysis and Data Processing

Chromatographic separation was performed using the Easy nLC system (Thermo Finnigan) with nanoliter flow rate HPLC. Buffer: the mobile phases were 0.1% FA in water (A) and 0.1% FA in ACN (B). The column was balanced with 95% A liquid. After the sample was injected into the Trap Column (Thermo scientific EASY column, 100 um*2 cm, 5 um, C18), it passed through the analytical columns (Thermo scientific EASY column, 75 μm*10 cm, 3 μm, C18) for gradient separation at a flow rate of 250 nl/min. The liquid phase separation gradient was as follows: 0 minutes–50 minutes, the linear gradient of liquid B was from 0% to 35%; 50 minutes–58 minutes, the linear gradient of liquid B was from 35% to 100%; 58 minutes–In 60 minutes, the linear gradient of liquid B rose to 100% and maintained. Perform mass spectrometry analysis using the Q Exactive HF-X mass spectrometer (Thermo Finnigan). Analysis time: 120 min, detection mode: positive ion. Parameters for the primary mass spectrometer: scan range 300–1800 m/z, mass spectral resolution: 60,000 (m/z 200), AGC target: 3e6, Maximum IT: 10 ms, dynamic exclusion time: 40 s. After each full MS scan (full MS scan), 10 fragmentograms (MS2 scans) were collected with the following settings: isolation window: 2 m/z, mass spectrum resolution: 17,500 (m/z 200), MS2 Activation Type: HCD, Normalized collision energy: 30eV, Underfill ratio: 0.1%. Finally, the raw data from mass spectrometry analysis were identified and quantitatively analyzed using the Mascot2.2 and Proteome Discoverer1.4 software.

Bioinformatics Analysis

KEGG enrichment analysis and GO enrichment analysis were conducted on differentially expressed proteins using the clusterProfiler, Blast2GO, STRING (https://string-db.org/) and KEGG Automatic Annotation Server (KAAS) software tools. Bubble charts were generated using the ggplot2 package in the R (version 4.5.2). Kinases of differentially phosphorylated proteins were predicted using the PhosphoSitePlus (PSP) database. Only predictions with a Site Percentile < 5 or the highest relative rank were retained for subsequent analysis.^10,11

Comprehensive Scoring Analysis

Composite scoring analysis: PPI data were downloaded from the STRING database and filtered to retain high-confidence interactions with combined score ≥ 0.7. An undirected, weighted network was constructed in R using the igraph package (v2.2.1), with edge weights set to the STRING combined score normalized to 0–1. Four network centrality metrics were calculated: betweenness centrality (B), closeness centrality (C), degree centrality (D), and eigenvector centrality (E). A composite network score for each node was computed as: NetworkScore(v)=10 × (0.30 × D + 0.40 × B + 0.20 × C + 0.10 × E). Multi-dimensional scoring framework: Using tissue-specific proteome and phosphoproteome data, a multi-tier scoring system was established: Protein expression score and phosphorylation validation score were derived from fold change (FC) and statistical significance (P-value). Network topology was represented by the composite NetworkScore. Biological validation score was assigned according to phosphorylation detection status and significance tiers. Final integration: The four component scores were integrated into a final composite score using the following weights: protein expression 45%, network topology 30%, phospho validation 15%, and biological validation 10%.

Statistical Analysis

All statistical analyses were performed using SPSS 22.0. Student t test was used for analyzing the data conforming to normal distribution, while non-parametric test was used for analyzing the data not conforming to normal distribution. In all cases, P<0.05 was considered significant and P<0.001 was considered highly significant. The prefix “p-” was used to indicate that the data originated from phosphoproteomics analysis.

Results

Identification of Differentially Expressed Proteins in Proteomic and Phosphoproteomic at Different Groups

In this study, HCT, ANT, and LCT tissues were analyzed by relative quantitative proteomics and TMT labeled quantitative phosphoproteomics to reveal important molecular features across different dimensions of liver diseases (Figure 1A). A total of 68,315 peptides and 7,198 proteins were identified in proteomics (Figure 1B). Differentially expressed proteins (DEPs) were screened based on an FC greater than 1.5 times (up-regulation greater than 1.5-fold or down-regulation less than 0.667-fold) and P<0.05, which effectively separated each group (Figure 1C). In HCT vs. ANT, there were 2,156 DEPs, with 1,595 up-regulated, and 561 down-regulated (Figure 1D); in HCT vs. LCT, there were 1,982 DEPs, with 1,392 up-regulated, and 590 down-regulated (Figure 1E); in ANT vs. LCT, there were 100 DEPs, with 19 up-regulated, and 81 down-regulated (Figure 1F). A total of 2,077 peptides were identified in phosphoproteomics, distributed across 1,237 proteins, encompassing 2,425 phosphorylation sites (Figure 1G). The proportions of phosphoserine (pSer), phosphothreonine (pThr), and phosphotyrosine (pTyr) were 87.34%, 11.75%, and 0.91%, respectively (Figure 1H). Additionally, proteins with only one phosphorylation site comprised the largest category, at 85.65% (Figure 1I). Differentially expressed phosphopeptides were screened with an FC greater than 1.2 times (up-regulation greater than 1.2-fold or down-regulation less than 0.83-fold) and P<0.05, effectively distinguishing each group (Figure 1J). In p-HCT vs. p-ANT, there were 235 differentially expressed phosphopeptides, with 127 up-regulated, and 108 down-regulated (Figure 1K); in p-HCT vs. p-LCT, there were 237 differentially expressed phosphopeptides, with 114 up-regulated, and 123 down-regulated (Figure 1L); in p-ANT vs. p-LCT, there were 82 differentially expressed phosphopeptides, with 25 up-regulated, and 57 down-regulated (Figure 1M). Thirteen DEPs and eight differentially phosphorylated peptides were shared among the three groups. The expression of differential proteins IL4I1, BCAT1, FAU, VNN1, FKBP8, POGZ, ARHGAP12, TIMM17B, and the expression of differential peptides aLGsPTk (NANS), tQLWASEPGTPPLPTsLPSQNPILk (SRRT), sGsSSPDSEITELk (CTPS1), lRLsPSPTSQR (LMNA), lRLsPsPTSQR (LMNA) were the highest in HCT, the expression of differential proteins TGFB1I1, CMC4, VSIG4, ENO3 and the differential phosphorylated peptides sPGADLLQVLTk (DAK), gGVDFALDcAGGsETmk (ADH4), sGtSEFLNk (HSP90B1) were the least expressed in HCT (Figure 1N and O).

Figure 1 Continued.

Figure 1 Summary of proteome and phosphoproteome data. (A) Workflow of proteomics and phosphoproteomics experiments. (B) Statistical results of proteome identification. (C) Cluster heatmap depicting DEPs in the proteome. (D) Volcano plot of HCT vs. ANT. (E) Volcano plot of HCT vs. LCT. (F) Volcano plot of ANT vs. LCT. In the volcano plots, blue indicates down-regulation and red indicates up-regulation. (G) Statistical results of phosphoproteome identification. (H) Distribution of identified phosphorylation sites on serine, threonine and tyrosine residues within the phosphoproteome. (I) Proportion of phosphorylation sites per protein within the phosphoproteome. (J) Cluster heatmap illustrating of differentially expressed phosphopeptides in the phosphoproteome. (K) Volcano plot of p-HCT vs. p-ANT. (L) Volcano plot of p-HCT vs. p-LCT. (M) Volcano plot of p-ANT vs. p-LCT. (N) Venn diagram of DEPs in the proteome. (O) Venn diagram of differentially phosphorylated peptides identified in the phosphoproteome.

Comparative Analysis of Differentially Expressed Proteins Between Proteomic and Phosphoproteomic

The alteration of expression in total protein is one of the primary reasons behind changes in phosphorylation modification levels. To gain insights into the shifts in expression levels of total proteins and phosphorylation modified proteins in the liver, the distribution of differentially expressed proteins (DEPs) and differentially expressed phosphorylated proteins (DEPPs) across various groups were compared. Thirty-two expressed proteins were overlapped between HCT vs. ANT and p-HCT vs. p-ANT, of which 20 up-regulated and 8 down-regulated expressions were consistent in both comparisons (Figure 2A). Similarly, forty expressed proteins were overlapped between HCT vs. LCT and p-HCT vs. p-LCT, of which 24 up-regulated and 11 down-regulated expressions were consistent in both comparisons (Figure 2B). Unlike the findings from the other two group comparisons, no overlapping DEPs between ANT vs. LCT and p-ANT vs. p-LCT (Figure 2C).

Figure 2 Overlap of the DEPs and DEPPs in Venn diagrams with different comparisons. (A) HCT vs. ANT. (B) HCT vs. LCT. (C) ANT vs. LCT.

GO Enrichment Analysis of DEPs and DEPPs

The top ten GO enrichment results demonstrated that the distributions of proteomics and phosphoproteomics differed across biological processes (BP), cellular components (CC), and molecular functions (MF). Proteomics predominantly centered on RNA and DNA processing and binding (Figure 3A–C), whereas phosphoproteomics mainly emphasized protein stability and the activity of structural molecules (Figure 3D–F). In addition, some functions were identified in both omics. In the comparison of HCT vs. ANT and p-HCT vs. p-ANT, significant differences were noted across several domains: In BP such as heterochromatin formation, telomere maintenance via telomerase, and glycolytic process; In CC such as nuclear matrix, chromosomal region, and nuclear membrane; In MF such as double-stranded DNA binding, transcription coactivator activity, and pyridoxal phosphate binding. In the comparison of HCT vs. LCT and p-HCT vs. p-LCT, significant differences were noted across several domains: In BP such as telomere maintenance via telomerase, cellular response to leukemia inhibitory factor, and glycolytic process; In CC such as nuclear matrix, nuclear envelope, and nuclear membrane; In MF such as structural molecule activity, minor groove of adenine-thymine-rich DNA binding, and arachidonate-CoA ligase activity. In the comparison of ANT vs. LCT and p-ANT vs. p-LCT, significant differences were noted across several domains: In BP such as response to food, wound healing, and cellular iron ion homeostasis; In CC such as spectrin complex; In MF such as structural constituent of cytoskeleton, ankyrin binding, and steroid hydroxylase activity (Table 2).

Table 2 GO Enrichment Analysis of Proteomic and Phosphoproteomic

Figure 3 Top 10 GO enrichment analysis of DEPs and DEPPs. (A) HCT vs. ANT. (B) HCT vs. LCT. (C) ANT vs. LCT. (D) p-HCT vs. p-ANT. (E) p-HCT vs. p-LCT. (F) p-ANT vs. p-LCT. The vertical axis represents biological processes (BP), cellular components (CC), and molecular functions (MF), respectively.

KEGG Enrichment Analysis of DEPs and DEPPs

In gain insights into the signaling pathways implicated in phosphorylation modified proteins and total proteins in the liver, KEGG enrichment analysis was conducted, revealing that the top ten enriched KEGG pathways varied between proteomics and phosphoproteomics. Proteomics was predominantly enriched for drug metabolism and immune-related pathways, with the chemical carcinogenesis-DNA adducts signaling pathway present across all three groups (Figure 4A–C); In phosphoproteomics, p-HCT vs. p-ANT and p-HCT vs. p-LCT groups were mostly enriched in glycolipid metabolism-related pathways, whereas p-ANT vs. p-LCT group was predominantly associated with apoptosis. Notably, the three groups did not have the same signaling pathway (Figure 4D–F). Furthermore, the same signaling pathways were identified in both omics analyses, involving some phosphorylated total proteins within these pathways. In HCT vs. ANT and p-HCT vs. p-ANT comparisons, the down-regulated protein fructose-bisphosphate aldolase B (ALDOB) was involved in the glycolysis/gluconogenesis pathway, the down-regulated protein PSAT1 was involved in the cysteine and methionine metabolism pathway, the up-regulated protein CTP synthase 1 (CTPS1) and the down-regulated protein phosphoserine aminotransferase (PSAT1) were involved in the biosynthetic of factors pathway, and the up-regulated protein voltage-dependent anion channel (VDAC1) was linked to the necroptosis pathway; In HCT vs. LCT and p-HCT vs. p-LCT comparisons, the down-regulated proteins ALDOB and fructose-bisphosphate aldolase C (ALDOC) were involved in the glycolysis/gluconeogenesis pathway, and the down-regulated protein solute carrier family 27 member 2 (SLC27A2) was involved in the peroxisome pathway. In contrast, the ANT vs. LCT comparison revealed enrichment of the chemical carcinogenesis-DNA adducts pathway in both omics, although the differences were not statistically significant, and no identical proteins or phosphorylated proteins were found to be involved in this pathway (Table 3).

Table 3 KEGG Enrichment Analysis of Proteomic and Phosphoproteomic

Figure 4 Top 10 KEGG enrichment analysis of DEPs and DEPPs. (A) HCT vs. ANT. (B) HCT vs. LCT. (C) ANT vs. LCT. (D) p-HCT vs. p-ANT. (E) p-HCT vs. p-LCT. (F) p-ANT vs. p-LCT.

Core Protein Scoring and Selection of DEPs and DEPPs

In the ANT vs. LCT comparison, the two omics layers showed no overlap in the sets of DEPs, and their significantly enriched KEGG pathways were non-overlapping. This indicates that between ANT and LCT, changes in protein abundance and regulation by phosphorylation do not exhibit coordinated patterns at the level of molecular sets or pathways. Therefore, to focus on molecular events with putative driving relevance to hepatocarcinogenesis, subsequent analyses concentrated on HCT vs. ANT and HCT vs. LCT. To systematically identify key driver proteins in HCT, we integrated differential protein expression, phosphorylation changes, and protein–protein interaction network topology to construct a multi-tier scoring framework, which produced the top 20 proteins for each comparison (Figures 5A and B) and their score distributions (Figures 5C and D). Proteins were ranked by the composite score; the top 10 were defined as core drivers, and ranks 11–20 as important regulators. Comparison of top proteins from HCT vs. ANT and HCT vs. LCT revealed 12 overlapping proteins-cytochrome P450 1A2 (CYP1A2), urocanate hydratase (UROC1), ALDOB, FAH, histone H1.5 (H1-5), ASS1, CYP2C8, glycogen synthase 2 (GYS2), alpha-aminoadipate aminotransferase (AADAT), formimidoyltransferase-cyclodeaminase (FTCD), CYP4A11, and DNA topoisomerase 3-beta-1 (TOP3B)-which are candidate shared drivers of HCC initiation and progression. KEGG enrichment analysis indicated that the top proteins in both comparisons were significantly enriched in metabolic pathways, biosynthesis of amino acids, retinol metabolism, and steroid hormone biosynthesis (Figures 5E and F). Notably, GAPDH displayed a distinctive regulatory pattern: its protein abundance was unchanged in HCT vs. LCT, whereas its phosphorylation was significantly downregulated (FC=0.59, P=0.032), suggesting that phosphorylation-dependent modulation of GAPDH function may play a dominant regulatory role during the transition from cirrhosis to HCC.

Figure 5 Scoring and selection of core proteins from DEPs and DEPPs in HCT vs. ANT and HCT vs. LCT. (A and B) Top 20 proteins ranked by composite score. (C and D) Score distribution plots: NetworkScore on the x-axis and Expression Score on the y-axis, each point represents one protein. (E and F) KEGG enrichment analysis of the top 20 proteins.

Kinase–Substrate Relationships Among Core Proteins

Changes in phosphorylation levels were primarily attributable to alterations in total protein abundance and upstream kinase activity, with kinases playing a key role in phosphorylation-mediated activation. For proteins with multiple phosphopeptides, the phosphopeptide exhibiting the largest FC was selected for analysis. Among the top 20 core proteins from each of the two comparisons, we identified eight core substrate proteins that contained specific phosphopeptides. To elucidate upstream regulatory mechanisms, kinase prediction was performed on these phosphopeptides using the PSP database. To ensure prediction reliability, stringent selection criteria were applied: only high-confidence predictions with Site Percentile < 5 or the top-ranked prediction were retained (Table 4). As shown in Table 4, this yielded a set of putative kinase–substrate regulatory pairs. Based on kinase family membership and annotated functions, the predicted kinases were categorized as metabolic regulators (Ribosomal S6 kinase 1 (RSK1), TAO kinase 3 (TAO3)), cell-cycle and proliferation regulators (Cyclin-dependent kinase 2 (CDK2), Aurora A (AURA), PDZ-binding kinase (PBK)), fibrosis/differentiation regulators (Bone morphogenetic protein receptor type II (BMPR2), Bone morphogenetic protein receptor type 1B (BMPR1B)), and apoptosis/stress regulators (Death-associated protein kinase 1 (DAPK1)). Notably, most predicted kinase–substrate relationships were not enriched within the same canonical KEGG pathways, suggesting extensive cross - pathway phosphorylation signal reprogramming during HCC progression. These predictions provide novel insights into the complex signaling network interactions in HCC and nominate several regulatory axes for experimental validation.

Table 4 Kinase Prediction Analysis of Core Proteins

Discussion

HCC remains a significant global public health concern due to its high mortality rate and limited treatment options.¹² The etiology of HCC is multifaceted, and identifying the critical targets involved in the transition from cirrhosis to HCC continues to pose a challenge. Genomics technologies are extensively employed to unravel disease mechanisms. Proteins are the products of gene expression, yet proteomics offers a more comprehensive view than genomics alone with its intricate composition and dynamic functionalities. Particularly, post-translational modifications highlight the significance of proteomics, serving as supplementary information to genomics. Protein phosphorylation, a key regulatory mechanism in tumorigenesis, represents a promising therapeutic target for drug intervention.¹³ However, elucidating the connections and discrepancies between multi-omics approaches in HCC requires further investigation. TMT-labeled phosphoproteomics is a well-established technique widely applied in various diseases. Studies employing TMT labeling technology have identified seven kinases, including GRP78, WSTF, PKN2, PRP4, LOK, NEK1, and AMPKA1, in non-functional pituitary neuroendocrine tumors (NF-PitNETs). The substrates of these key kinases may offer novel insights into discovering drug targets for NF-PitNETs.¹⁴ Currently, there is a scarcity of research utilizing TMT labeling for phosphoprotein analysis in liver tissues. In this study, we conducted a combined analysis of relative quantitative proteomics and TMT-labeled quantitative phosphoproteomics for HCT, ANT, and LCT, aiming to identify crucial biomarkers pertinent to the onset and progression of liver cancer.

It is well known that proteomics has been employed in analyzing protein expression profiles in HCC. However, protein activity is also regulated by post-translational modifications, which may not manifest as significant changes at the protein level but rather at the post-translational modification level. Consequently, integrating proteomics and phosphoproteomics analysis provides a more robust foundation for identifying tumor markers. Our integrated multi-omics analyses showed that, during the transition from LCT or ANT to HCT, there are core hub proteins and pathways exhibiting coordinated regulation of protein expression and phosphorylation-mechanisms that correspond to drivers of malignant transformation and the establishment of tumor adaptive programs. By contrast, ANT vs. LCT lacked such coordinated changes, indicating that the two follow fundamentally different pathological trajectories: ANT reflects a pro-tumor microenvironment, whereas LCT represents a state of persistent precancerous injury. This finding suggests that hepatocarcinogenesis and progression are not linear processes but rather arise from the convergence of distinct driver events, and therefore stage-specific investigation is warranted. Given the lack of coordinated patterns between ANT and LCT, subsequent analyses were focused on HCT vs. ANT and HCT vs. LCT. To systematically identify core driver proteins in HCT, we integrated protein expression, phosphorylation changes, and PPI relationships to construct a multi-tier scoring system and selected the top 20 proteins for each comparison; 12 core driver proteins were highly overlapping between the two sets. Importantly, these proteins were not randomly distributed but were significantly enriched in liver-related metabolic pathways, encompassing lipid, amino-acid, and carbohydrate metabolism.

At the lipid metabolism level, we identified cytochrome P450 family enzymes CYP4A11 and CYP2C8. CYP4A11, a fatty acid ω-hydroxylase and a rate-limiting enzyme that initiates fatty acid catabolism, has been shown to correlate strongly with plasma levels of lipid peroxidation products.¹⁵ In our study, CYP4A11 was downregulated in HCT vs. ANT (FC=0.19, P=6.83E-13) and HCT vs. LCT (FC=0.17, P=4.68E-13). This result is consistent with prior reports of reduced CYP4A11 expression in HCC and its positive association with favorable prognostic features such as smaller tumor size, lower histologic grade, and earlier pathological stage, suggesting that CYP4A11 expression may represent an important regulatory node in cancer progression.¹⁶ CYP2C8, another key cytochrome P450 enzyme, was also downregulated (HCT vs. ANT: FC=0.17, P=8.48E-15; HCT vs. LCT: FC=0.18, P=8.05E-14). Although CYP2C8 overexpression has been reported to suppress malignant HCC phenotypes, including proliferation, migration, invasion, and sorafenib resistance, via inhibition of the PI3K/Akt/p27Kip1 axis.¹⁷ CYP2C8 can also metabolically inactivate chemotherapeutics such as paclitaxel, and its inhibition may sensitize cancer cells to these agents.¹⁸ The dual roles of CYP2C8 in drug response likely reflect its functions in signaling regulation and drug metabolism, and support its potential utility for guiding drug selection and optimizing therapy based on expression level.

Amino-acid and carbohydrate metabolism are also central hubs of hepatic metabolism. Amino-acid metabolism-related proteins include ASS1 and FAH. ASS1, the rate-limiting enzyme of the urea cycle and arginine synthesis, was downregulated in HCT. Sanghwa Kim et al reported similar downregulation and showed that ASS1 overexpression suppresses tumor growth by activating the PERK/eIF2α/ATF4/CHOP axis and potentiates the efficacy of combined chemotherapy in vitro and in vivo, suggesting ASS1 upregulation as a potential therapeutic strategy.¹⁹ FAH, the terminal enzyme in tyrosine catabolism, was also downregulated; hereditary loss of FAH causes hereditary tyrosinemia type I, leading to accumulation of toxic intermediates (eg., succinylacetone, maleylacetone, fumarylacetoacetate) in liver and kidney tissues and conferring markedly elevated HCC risk.^20,21 The downregulation of FAH observed in our HCT samples may perturb tyrosine catabolism, promote toxic intermediate accumulation and metabolic rewiring, and thereby play a unique role in HCT initiation or maintenance. Notably, FAH was detected as phosphorylated in our phosphoproteome, suggesting PTM-dependent regulation. At the carbohydrate metabolism level, ALDOB downregulation is commonly associated with liver tumor metabolism.²² Importantly, our phosphoproteomic data indicate that ALDOB phosphopeptides were significantly less abundant in HCT than in LCT and ANT. Although phosphorylation of its isozyme ALDOA has been extensively studied, post-translational modification-particularly phosphorylation-of ALDOB is rarely reported, and its functional significance in liver physiology and pathology remains unclear.²³

Another core enzyme significantly enriched in the Glycolysis/Gluconeogenesis pathway with ALDOB is GAPDH. Our data show that, in p-HCT vs. p-LCT,GAPDH phosphorylation is significantly decreased while total GAPDH protein abundance is unchanged. This regulatory pattern suggests that, at the critical transition from cirrhosis to HCT, GAPDH function may be predominantly regulated by phosphorylation. Although GAPDH has been implicated across diseases, the regulatory mechanisms controlling GAPDH phosphorylation in HCT remain unreported. Recent studies indicate that GAPDH ubiquitination can directly lead to reduced phosphorylation and enzymatic activity, thereby inhibiting glycolysis, whereas GAPDH mutations can enhance glycolysis and accelerate tumorigenesis in cervical cancer.²⁴ These findings are highly consistent with our data and support the hypothesis that cells may actively suppress GAPDH activity through PTM networks during the cirrhosis to HCT transition.

Moreover, our kinase predictions identified several kinase–substrate pairs, including RSK1-ALDOB, DAPK1-MAP4, and AURA-GAPDH. However, these kinase–substrate relationships were not enriched within the same canonical KEGG pathways, suggesting cross-pathway phosphorylation signal reprogramming during HCT progression. Although direct modification evidence is lacking for many candidate kinases, their functional pathways partially intersect with substrate metabolic pathways. For example, RSK, a downstream kinase of the MAPK pathway, promotes glycolytic metabolism and tumor growth in BRAF-mutant melanoma by phosphorylating a key glycolytic enzyme.²⁵ AURA kinase, a key regulator of the G2/M cell-cycle transition, can regulate pyruvate kinase M2 (PKM2) activity during S phase to promote anabolic glycolytic branch metabolism.²⁶ These indirect functional associations provide potential support for the predicted kinase–substrate regulatory relationships and offer new mechanistic clues for the regulation of the core proteins identified here. Based on the outcomes of this study, we intend to expand the sample size in future research to further validate these findings.

Conclusion

By integrating relative quantitative proteomics and TMT-labeled quantitative phosphoproteomics, this study systematically and, to our knowledge, for the first time delineates differences in the coordinated regulation of protein abundance and phosphorylation between HCT, ANT, and LCT, clarifies the non-linear, multi-pathway convergence features of HCT progression, and identifies a set of key driver proteins enriched in core hepatic metabolic pathways. From the perspectives of metabolic reprogramming and post-translational modification regulation, these findings deepen mechanistic understanding of HCT pathogenesis and provide molecular leads for identifying stage-specific biomarkers and therapeutic targets. Notably, core proteins such as CYP4A11, CYP2C8, ASS1, FAH, ALDOB, and GAPDH-especially the aberrant phosphorylation of ALDOB and GAPDH-represent promising candidates for novel biomarkers and targeted therapies in HCT. However, additional samples are required to validate these findings, and further research is necessary to bridge the gap from basic science to clinical application.