Transcriptome data and somatic mutation matrices of 14 types of cancer were collected from the University of California Santa Cruz (UCSC) database (https://xena.ucsc.edu/) (8). Meanwhile, corresponding clinical data were also obtained ( Supplementary Table S3 ). The list of 500 aging-related messenger RNAs (mRNAs) was downloaded from the Aging Atlas database (9) (https://ngdc.cncb.ac.cn/aging/index) ( Supplementary Table S1 ). To validate the conservatism of this subtype, additional independent LIHC validation sets, GSE76427 (10), GSE14520 (11), and GSE54236 (12), containing 167, 43, and 161 samples, respectively, were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). In addition, we also downloaded single-cell dataset, GSE149614 (13), to investigate the function of aging-related genes in the tumor microenvironment. In order to verify the expression of CCL19 in CTNNB1 mutation set and CTNNB1 wild-type set, we also downloaded the transcriptome data and corresponding clinical information of GSE9829 (14). In addition, to verify the credibility of the prognostic model, we downloaded the transcriptome data and corresponding clinical information of GSE20140 (15) from the GEO database.

For the purpose of identifying recurrent driver genes, we involved MutSig2CV (16), a tool that needs to be performed in MATLAB. In this analysis, MutSig2CV performs three tests of significance for each gene based on abundance (mutation rate), cluster (physical co-location of mutations), and evolutionary conservation. P value <0.05 and q value <0.05 were used as the criteria for screening recurrently mutated genes.

To delineate the association between subtypes and immunity, we downloaded the infiltration estimation matrices of immune associated cells of patients with LIHC, which were calculated by CIBERSORT (17), XCELL (18), TIMER (19), and MCPCOUNTER (20) algorithms, respectively, from TIMER2.0 database.

PROGENy (21), a method that can infer pathway activity from gene expression by using core pathway-responsive genes, which were computed by leveraging a large compendium of publicly available perturbation experiments, was involved to calculate the activity of various pathways of LIHC patients, including the Androgen, Estrogen, epidermal growth factor receptor(EGFR), Hypoxia, Janus kinase-signal transducers and activators of transcription (JAK-STAT), mitogen-activated protein kinases (MAPK), nuclear factor-κB (NFκB), PI3K, protein 53 (p53), transforming growth factor β (TGFβ), Tumor necrosis factor α (TNFα), tumour necrosis factor related apoptosis-inducing ligand (TRAIL), Vascular endothelial growth factor (VEGF), and Wingless-Type (WNT) pathways.

In order to probe the biological pathways associated with aging-related genes in the model, gene set enrichment analysis (GSEA), a knowledge-based approach for interpreting genome-wide expression profiles, was performed in the high-expression and the low-expression groups compared with the median level of gene expression (22).

The Human Protein Atlas (23) (HPA; https://www.proteinatlas.org/) database was used to explore the protein level of EEF1E1 and HDAC2 in hepatocellular carcinoma and normal tissue samples. In addition, we collected protein expression labeling information (weak/moderate/high) from the HPA website of EEF1E1 and HDAC2 in immunohistochemical figures and calculated the significance between the number of normal tissue samples and hepatocellular carcinoma samples in protein expression level by using the chi-square test.

Seurat v.3.0.0 (24) was involved for data normalization, dimensionality reduction, and clustering using default parameters. In preprocessing, cells are filtered based on the criteria of expressing a minimum of 200 genes. Also, cells that had more than 20% expression on mitochondrial genes were also removed. Cell clusters were annotated using singleR (25) with classical cell signatures.

The total cellular proteins from each group were extracted using Radio Immunoprecipitation Assay (RIPA) lysis buffer with 1% phenylmethanesulfonyl fluoride (PMSF). Then, equal amounts (20 μg) of protein determined by Bicinchoninic Acid (BCA) protein assay kit (Thermo Fisher Scientific, Waltham, MA, USA) were separated using 10% sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS-PAGE) gels. The proteins were then transferred to polyvinylidene difluoride(PVDF) membranes (0.45 mm, Solarbio, Beijing, China). The membranes were blocked with 5% nonfat milk for 1 h at room temperature and then incubated with primary antibodies at 4°C for 12 h. The following antibodies were tested: Histone Deacetylase 2(HDAC2), glyceraldehyde-3-phosphate dehydrogenase(GAPDH) (1:1,000; Proteintech Group Inc.) rabbit polyclonal antibodies. The secondary antibodies were anti-mouse or anti-rabbit antibody and conjugated to horseradish peroxidase (HRP) (1:4,000; Proteintech Group Inc.). The secondary antibodies were used at a 1:4,000 dilution and were incubated for approximately 1 h at room temperature. The bands were visualized with Enhanced chemiluminescence (ECL) reagents (Thermo Fisher Scientific) and developed by Omega Lum G (Aplegen).

Total RNA was extracted by TRIzol Reagent (Invitrogen) from cells. Complementary DNA (cDNA) was obtained from total RNA with PrimeScript™ RT reagent kit (Takara Bio, Inc., Otsu, Japan). mRNA expression was assessed by real-time quantitative PCR (RT-qPCR), which was carried out in triplicate by a SYBR Premix Ex Taq™ kit (Takara Bio) and ABI 7900HT Real-Time PCR system (Applied Biosystems Life Technologies, Foster City, CA, USA). The primers used are shown in Supplementary Table S3 .

Secondary liver cancer samples were collected, embedded, and sectioned. Immunohistochemical staining was performed according to standard protocol, and the primary antibody was anti-HDAC2 (1:100, Proteintech). The antibody catalog in the Western blot and immunohistochemistry (IHC) analysis is 12922-3-AP (Proteintech Group Inc.). The staining scores were calculated using the immunoreactive score (IRS) system. The percentage of positive cells was scored as follows: no stained cells = 0; <10% staining = 1; 10%–50% staining = 2; 51%–80% staining = 3; and >80% staining = 4. The staining intensity was scored as follows: no color reaction = 0; mild reaction = 1; moderate reaction = 2; and intense reaction = 3. The IRS scores = (scores of staining intensity) × (scores of percentage of positive cells).

According to the aging-related genes that were downloaded from the Aging Atlas database (9), we divided LIHC patients into three groups by applying R package “CancerSubtypes” (26). The survival rate of each subtype was calculated by using the Kaplan–Meier (K-M) method, and the log-rank test was used to assess the difference in survival among the three subtypes with a criterion level of P < 0.05. Wilcoxon test and Kruskal–Wallis test were involved to estimate the subtype-specific differences and the mRNA expression differences related to other factors. Briefly, when clinical traits have two characteristics, Wilcoxon test was involved to examine the significance. Meanwhile, Kruskal–Wallis test was applied when the clinical features have more than two characteristics. “Limma” (27) package was used to calculate the marker genes of each subtype. Furthermore, hierarchical cluster analysis was performed to classify the validation set via the ward.D method. In addition, the difference between the proportion of clinical factors and subtypes was computed by the chi-square test. Search Tool for the Retrieval of Interaction Gene (STRING) database was applied to calculate the relationship between marker genes of each subtype, and Cytoscape was involved to describe the network of these genes. The relationships between the expression of mRNAs and immune infiltrating cells were calculated by Spearman’s correlation coefficient. Aging-related mRNAs significantly associated with survival were identified by using univariate Cox proportional hazards regression. Then, least absolute shrinkage and selection operator (LASSO) regression analysis and multivariate Cox regression was applied to establish the aging-related prognosis model based on aging-related genes. Multivariate Cox regression was also used to assess the independence of the model compared with other clinical factors such as Grade, T stage, and Stage. Also, the credibility and predictive value of our model were evaluated through time-related receiver operating characteristic (ROC) curve. Moreover, waterfall plot was created by using the maftools (28) package. All statistical analyses were performed using R-version 4.1.1.

Principal component analysis (PCA) was applied to detect the expression patterns of aging-related mRNAs in 14 types of cancer containing LIHC from The Cancer Genome Atlas(TCGA). Based on the first two principal components (PCs) of PCA, LIHC was characterized by a unique expression pattern of aging-related genes, which was completely segregated from other cancer types ( Figure 1A ). Thus, patients with LIHC may have unique expression patterns of aging-related genes. In order to characterize the heterogeneity in patients with LIHC, CancerSubtypes, a method of classing, was implied to identify the molecular subtypes of LIHC patients. As a result, LIHC patients were divided into three aging microenvironment (AME) subtypes (AME_I [I = {1,2,3}]) with completely different prognosis characteristics ( Figure 1B , Supplementary Table S1 ). AME2 patients had the highest survival rate, whereas patients in AME3 subtype had the worst outcomes. To uncover the characteristics of age in these subtypes, we defined patients older than 65 as senility and those younger than 55 as youth. Using the chi-square test, we found that patients in the AME1 subtype had the highest proportion of senility ( Figure 1C ).

For the purpose of characterizing the molecular basis of these three aging microenvironment subtypes in LIHC, we downloaded the infiltration estimation matrix of immune cells in patients with LIHC, which were calculated by TIMER and MCPCOUNTER and so on, from TIMER2.0 database (http://timer.cistrome.org/). Using different algorithms, we found a consistent result that in the AME1 subtype, almost all types of immune infiltrating cells such as CD8 T cells, CD4 T cells, and Natural Killer (NK) cells were significantly suppressed, whereas samples in the AME2 subtype showed high infiltration of immune cells ( Figure 1D ). Meanwhile, the AME3 group was characterized by intermediate feature. Thus, we characterized AME1 as “Immune Exclusion” feature, AME2 as “Immune Infiltration” feature, and AME3 as “Immune Intermediate” feature, which were exactly consistent with the results of previous studies (29). Interestingly, the immune score, calculated by XCELL, also significantly activated in the AME2 subtype ( Figure 1E ). Also, immune infiltration cells, related to the inhibition of tumor cells, including Macrophage, B cell, CD4 T cell, CD8 T cell, and NK cell, were significantly activated in the AME2 subtype ( Figure 1D ). Also, the purity of tumor cells was significantly decreased in AME2 ( Figure 1F ). Nevertheless, immune infiltration cells, related to the proportion of tumor cells, just like regulatory T cells (Tregs) and myeloid-derived suppressor cells (MDSCs), were noteworthy suppressed in the AME2 subtype ( Figures 1G, H ). This result was consistent with the highest survival rate of the AME2 subtype.

Somatic mutation, a molecular characteristic associated with aging, has also been proven to be associated with immune cells (30). By applying MutSig2CV, we identified 15 recurrently mutated genes (q < 0.05), 12 of which were found in at least 10 patient samples ( Figure 2A , Supplementary Table S1 ). The top putative candidate drivers from these patients with LIHC were Tumor Protein P53 (TP53), AT-Rich Interaction Domain 2 (ARID2), BRCA1-Associated Protein 1 (BAP1), Albumin (ALB), RB Transcriptional Corepressor 1 (RB1), Axis Inhibition Protein 1 (AXIN1), Kelch-Like ECH-Associated Protein 1 (KEAP1), Bromodomain-Containing Protein 7 (BRD7), Ribosomal Protein S6 Kinase A3 (RPS6KA3), Activin A Receptor Type 2A (ACVR2A) ( Figure 2A ). Specifically, the mutation proportion of CTNNB1 was the highest in the AME1 subtype, the mutation proportion of TP53 and BAP1 was the highest in the AME3 subtype, and the mutation proportion of ARID2 was the highest in the AME2 subtype ( Figure 2B ).

Immune checkpoints, as primary immune therapeutic targets, played important roles in the human immune system and have been widely applied in clinical treatment, which has achieved good clinical outcomes (29, 31). Therefore, we investigated the expression characteristics of immune checkpoints in these three subtypes. According to the expression of immune checkpoints, we find that almost all immune checkpoint expression levels were significantly higher in the AME2 subtype than those in the other two subtypes, containing Cytotoxic T Lymphocyte-Associated Protein 4 (CTLA4), T-Cell Immunoglobulin Mucin 3 (TIM3), Programmed Cell Death 1 (PD-1), Indoleamine 2,3-Dioxygenase 1 (IDO1), Programmed Cell Death 1 Ligand 2 (PD-L2), T-Cell Immunoreceptor with Ig and ITIM domains (TIGIT), and Programmed Cell Death 1 Ligand 1 (PD-L1) ( Figure 2C ). However, the expression of PD-1/PD-L1 was significantly restrained in the AME1 subtype, which was in accordance with previous research. On the contrary, the expression levels of PD1/PD-L1 and T cells in the AME2 subtype were significantly increased ( Figures 2D , E ). These results indicated that patients within the immune infiltration subtype may most likely respond to immune checkpoint inhibitors (ICIs), while those within the immune exclusion subtype might have innate resistance to anti-PD1/PD-L1 or similar therapies (32). In addition, we also explored the association between the mutation of TP53 and CTNNB1 and clinical factors including Stage, T, M, N, and status. Consequently, we found that there was no significant association between TP53 and CTNNB1 mutation and clinical factors ( Supplementary Figures S1A , S1B ). Moreover, we also demonstrated the amino acid mutation sites of TP53 and CTNNB1 ( Supplementary Figures S1C , S1D ). By applying chi-square test, we found that the mutation rate of TP53 and CTNNB1 had significant differences among the three subtypes (P < 0.05), among which TP53 mutation rate was dramatically increased in the AME3 subtype ( Figure 2F ), while the CTNNB1 mutation rate was significantly increased in the AME1 subtype ( Figure 2G ). In agreement with previous research, the effects of PD-1/PD-L1 ICI therapy in LIHC patients with CTNNB1 mutation were very poor (32). This suggests that CTNNB1 may be responsible for the characteristics of immune exclusion in patients with the AME1 subtype.

Many previous classification methods failed to classify independent datasets other than the training datasets. We used the aging-related gene signatures in a pipeline that could perform pan-cancer classification. We designed our pipeline using supervised clustering. This pipeline can be divided into five steps: 1) mRNA expression profile between the AME_j (j={1,2,3}) group and other subtypes AME_m (m={1,2,3} & m!=j) and AME_k (k={1,2,3} & (k!=j) & (k!=m)) was compared by using limma package. 2) |LogFC| >0 and False discovery rate (FDR) <0.05 were used as the criteria for screening differentially expressed aging-related mRNAs. 3) Duplication genes that were specifically upregulated in each subtype were removed from each group. Then, the top 50 aging-related genes in each subtype were selected as markers in each group according to their logFC. 4) Then, samples were divided into three subtypes through hierarchical clustering, among which “ward.D” was selected as the method.

We calculated each subtype’s markers using the pipeline shown above ( Supplementary Figures S2A - S2C , Supplementary Table S1 ). Moreover, we also described the correlation between these markers by applying Spearman correlation coefficient. As shown in Supplementary Figure S2D , these markers can be divided into three groups, which suggested that the markers of each subtype were interrelated. Moreover, we also investigated the protein–protein intersection network between these markers by using STRING software ( Supplementary Figure S2E ). Then, three independent validations of LIHC transcriptomic data (GEO: GSE76427, GSE14520, GSE54236) were analyzed to determine whether our subtypes were conserved. In accordance with the classification characteristics of our subtypes, the three groups of GEO samples were divided into three subtypes: AME1, AME2, and AME3, indicating that the expression signature based on the aging-related genes can be applied to independent LIHC samples, among which the immune signatures in the AME2 subtype were higher than that in the AME1 and AME3 subtypes, and other features such as the Epithelial–mesenchymal transition (EMT) signature were not significant in the three subtypes ( Figures 3A–C , Supplementary Table S1 ). In addition, by applying the pipeline, three additional types of tumors [Colon adenocarcinoma (COAD), Breast invasive carcinoma (BRCA), Lung adenocarcinoma (LUAD] in TCGA were classified into the three subtypes ( Figure 3D , Supplementary Table S1 ), which had significant correlation to overall survival (OS) (P < 0.05), indicating that this classification pipeline can be broadly applied at the pan-cancer level ( Figure 3E ).

We next compared our subtyping method with three other classification methods including TME subtype (33), Immune-Subtype (34), and iCluster (35) in TCGA-LIHC samples ( Figures 4A, C, E ; Supplementary Table S1 ). The proportion of C3 subtype is the largest in AME2, while the proportion of the C1 subtype is the largest in AME3. Accordingly, the C3 subtype has the highest survival rate and the C1 subtype has the lowest survival rate, which are consistent with the characteristics of the AME2 subtype and AME3 subtype, respectively ( Figures 4A, B ). When analyzed within the two subtypes separately, our AME classification model exhibited more significant correlation with OS than that in the Immune-Subtype method ( Figures 1B , 4B ). Consistent with the previous results, the proportion of IE (Immune-Enriched Fibrotic) subtype and IE/F (Immune-Enriched, Non-Fibrotic) subtype was the highest in AME2, which further confirmed that AME2 is characterized by immune infiltration. Meanwhile, the survival rates of the IE subtype and IE/F subtype were significantly higher than those of D (Depleted) subtype and F (Fibrotic) subtype ( Figure 4C ). When compared to our AME classification model, TME subtypes have poor survival segregation in LIHC ( Figures 1B , 4D ). Lastly, survival rates among iCluster subtypes were not statistically different ( Figure 4F ).

Overall, our classification system using the aging microenvironment-associated gene expression pattern across various cancer types divided the AME into the three subtypes that significantly correlated with OS in comparison with other LIHC classification approaches including Immune-Subtype, TME, and iCluster.

Previously, the hypothesis that LIHC “cold” tumors defined by Wnt/CTNNB1 mutation are those refractory to ICIs has been confirmed by Pinyol et al. (32). Immunotherapies in patients with “hot” tumors that were characterized by the presence of tumor-infiltrating lymphocytes (TILs) were more likely to be efficacious than in patients with “cold” tumors (36–38). About 30% of patients’ enrichment in CTNNB1 mutation showed exclusion of TILs (36, 37). Consequently, tumors within the LIHC Immune Exclusion subtype might represent those with innate resistance to anti-PD1/PD-L1 or similar therapies. Therefore, it is necessary for us to explore the causes of immune checkpoint resistance in patients with CTNNB1 mutation. To check whether the mutation of CTNNB1 can affect the expression of immune checkpoints, we divided LIHC patients into CTNNB1 mutated group and CTNNB1 wild-type group. Then, by applying Wilcoxon test, we find that a large number of immune checkpoints were significantly differentially expressed between these two groups ( Figure 5A ). The same phenomenon was also verified in the GSE9829 dataset ( Figure 5B ). As shown in Figure 2G , the proportion of patients with CTNNB1 mutation in the AME1 subtype was significantly higher than that of the AME2 subtype and AME3 subtype. Thus, the mRNA expression profile based on aging-related genes between the AME1 and other subtypes AME_O (O={2,3}) was compared by using limma package. |LogFC| >0.5 and FDR <0.05 were used as the criteria for screening differentially expressed aging-related mRNAs. Among the identified differentially expressed mRNAs, 87 downregulated mRNAs and 26 upregulated mRNAs in the AME1 subtype were selected as candidate markers ( Supplementary Figure S2A ). By applying Spearman’s correlation coefficient, we finally screened 11 aging-related genes that have a strong relationship with the activation of T or B cells (|Cor| >0.4 and P < 0.05), including ZAP70, PIK3CD, PRKCQ, PRKCB, TNFSF13B, BTK, SYK, TNFRSF13C, CAMK4, CCL19, and TNFRSF11A ( Figure 5C ). Interestingly, the expression of these genes was significantly restrained in the AME1 subtype. By using Wilcoxon test, PRKCB was excluded because it was not significantly differentially expressed between wild-type and mutated CTNNB1 AME1 samples.

By comparing the correlation between the 10 screened genes and immune checkpoints (PD-1/PD-L1) in the CTNNB1 mutated group and the wild type in AME1 patients ( Figures 5D, E ), we found that the correlation between these genes and immune checkpoints (PD-1/PD-L1) in the mutant group notably increased ( Supplementary Figures S3A , S3B ). Especially, the correlation of CCL19, TNFRSF11A, TNFSF13B, and SYK with immune checkpoints (PD-1/PD-L1) in the CTNNB1 mutant group was significantly higher than that in the CTNNB1 wild type (Cor_mut-Cor_wild >0.1 and P < 0.05). Interestingly, the correlation between CCL19 and PD-L1 was not significant in the CTNNB1 wild type (ρ = 0.0435, P = 0.607), whereas in the CTNNB1 mutant group, the correlation between CCL19 and PD-L1 was significant (ρ = 0.254, P = 0.00339) ( Figures 5F, G ). Then, we investigated the correlation between the expression of these aging-related genes and CD8+ T cells ( Supplementary Figure S3C ). Consistent with previous results, the correlation between CCL19 and CD8+ T cells was significantly higher in the mutant group than that in the CTNNB1 wild-type group (ρmut = 0.681, P < 0.001; ρwild = 0.458, P < 0.001) ( Figure 5H ). Moreover, the expression of CCL19 was significantly decreased in the CTNNB1 mutation group ( Figure 5I ). We also found the same tendency in the external validation dataset: GSE9829 ( Figure 5J ). In addition, qPCR analysis was employed to verify this tendency. We found that the expression of CCL19 in CTNNB1 wild-type hepatocellular carcinoma cells (HepG2 cells) is higher than that in CTNNB1 mutant cells (Huh6 and SNU398 cells) ( Figure 5K , Supplementary Table S2 ). It has been reported that CCL19, as a chemokine, can promote the release of cytokines by CD8+ T cells by binding to its receptor CCR7, thereby inhibiting the proliferation, migration, and invasion of tumor cells (39). These results suggested that CTNNB1 mutation may inhibit the activation of immune cells by inhibiting the expression of CCL19, leading to the “immune cold” phenotype. Moreover, the correlation between the expression of CCL19 and immune checkpoints (PD-1/PD-L1) in patients with CTNNB1 mutation was significantly increased, which may also explain why AME1 subtype patients are resistant to anti-PD1/PD-L1 or similar therapies.

To explore the potential reasons for the poor prognosis in patients with the AME3 subtype, we calculated the activity of 14 cancer-related signaling pathways, including the Androgen, Estrogen, EGFR, Hypoxia, JAK-STAT, MAPK, NFκB, PI3K, p53, TGFβ, TNFα, TRAIL, VEGF, and WNT pathways, which are available in the R package “PROGENy” (21) ( Supplementary Table S1 ). By Kruskal–Wallis test, we identified significant differential expression among the three AME subtypes in all 14 pathways, among which the MAPK, PI3K, and VEGF pathways were specifically activated in the AME3 subtype, whereas the p53 pathway was specifically suppressed ( Figures 6A–D ). The percentage of tumors that carry TP53 mutation as one of the suppressors of the p53 signaling pathway (40) was significantly higher in the AME1 subtype than those in the AME2 and AME3 subtypes ( Figure 2G ). As expected, the p53 pathway was significantly suppressed in the TP53 mutation group ( Figure 6E ). Meanwhile, mTOR, a downstream molecule of the p53 pathway (41), was significantly activated in the TP53 mutation group ( Figure 6F ). These results suggested that the TP53 mutation may inhibit the activity of the p53 pathway, which promotes the expression of mTOR and then contributes to the proliferation, migration, and invasion of tumor cells.

In addition, we found that the PI3K pathway, another activator of mTOR, was specifically activated in the AME3 subtype ( Figure 6B ). To explore the reason for the specific activation of the PI3K signaling pathway, we investigated the relationship between the INPP4B mutation and the PI3K signaling pathway (42). As shown in Figure 6G , the activity of the PI3K signaling pathway and the expression of mTOR were significantly enhanced ( Figure 6H ), suggesting that the INPP4B mutation may promote the activation of the PI3K signaling pathway, thereby activating the expression of mTOR. Furthermore, we also computed the differentially expressed genes in these pathways including the MAPK, PI3K, VEGF, and p53 pathways between AME3 and other subtypes ( Figure 6I ). Collectively, we propose that specific activation of mTOR is responsible for the poor prognosis of AME3 subtypes; importantly, we uncovered several regulatory axes, including TP53/p53/mTOR and INPP4B/PI3K/mTOR ( Figure 6J ), providing clues for the treatment of the AME3 subtype.

In order to investigate the prognostic effects of these aging-related mRNAs, we established aging-related genes by applying LASSO and multiple Cox regression analysis. Different from the traditional methods, we adopted multiple iterations (times = 500) to avoid the contingency caused by the random error of LASSO regression. This process can be divided into five steps: 1) Aging-related genes associated with prognosis were screened out by univariate Cox regression, and those with low variance were removed (Var <5); 2) The selected aging-related genes were dimensionally reduced by LASSO regression, repeated 500 times, and the total number of occurrences of each gene was counted; 3) Genes, total occurrences, number of more than 100 times were retained; 4) The model with the maximum AUC value was selected as the final prognostic model through multiple repetitions; 5) Calculate the prognostic signature for each sample.

Through this process, we finally selected 27 genes to establish a model corresponding to the maximum AUC (5-year AUC = 0.799) ( Figures 7A, B ), including AKT2, ATG101, BMI1, BMP6, CCNA2, CD40LG, COQ7, CXCL8, DLAT, EEF1E1, EIF4E, GCLM, HDAC1, HDAC2, HTRA2, IL7R, MAPT, MMP1, PIK3R1, PON1, PPARGC1A, SOCS2, TFDP1, TRAF3, ZAP70, UBB, and JUND. The results from the K-M analysis suggested that patients with high risk have lower OS rate than patients with low risk in LIHC patients in TCGA (P < 0.05) ( Figure 7C ). In order to verify the advantage of this model, we applied these 29 genes to construct multiple models including Grading Boosting (1-year AUC = 0.790; 3-year AUC = 0.814; 5-year AUC = 0.786), Logistic Regression (1-year AUC = 0.810; 3-year AUC = 0.809; 5-year AUC = 0.785), K-Nearest Neighbor (1-year AUC = 0.774; 3-year AUC = 0.789; 5-year AUC = 0.741), Bayesian (1-year AUC = 0.810; 3-year AUC = 0.808; 5-year AUC = 0.74), Artificial Neutral Network (1-year AUC = 0.674; 3-year AUC = 0.670; 5-year AUC = 0.606), and Decision Tree (1-year AUC = 0.499; 3-year AUC = 0.426; 5-year AUC = 0.44), and we found that our model has dominant credibility and predictive value ( Figure 7D ). Univariate Cox regression results indicated that our model was significantly correlated with prognosis, and multivariate Cox regression results suggested that our model could be used as an independent prognostic factor that was independent from other clinical factors, including Stage, T, N, and M ( Figure 7E ). In addition, a cross-platform validation from GSE20140 was involved to verify the accuracy of the model. A total of 20 genes were included in the list of 29 genes in the GSE20140 validation set, and the risk model was established through multivariate Cox regression. The ROC curve prompted that this model had dominant credibility and predictive value (1-year AUC = 0.848; 3-year AUC = 0.824; 5-year AUC = 0.840) ( Figure 7F ). By constructing K-M analysis, this model was proven to have a great prognostic credibility (P < 0.05) ( Figure 7G ).

We further investigated molecular events at single-cell resolution by comparing high- and low-risk tumors categorized by our aging gene-based model. We annotated seven major cell types for 10 primary liver tumor tissues from the GSE149614 dataset using biomarkers ACTA2 (Fibroblast cells), CD68 (Myeloid cells), ALB (Hepatocyte cells), CD3D (T cells), JCHAIN (Plasmablast cells), MS4A1 (B cells), and PECAM1 (Endothelial cells) ( Supplementary Figures S4A , S4B ). By applying the AUCell (43) algorithm, we calculated the score of each cell through these 27 gene features that comprised our model and determined HCC01T, HCC06T, HCC07T, and HCC09T samples as high-risk samples and HCC02T and HCC04T as low-risk samples ( Supplementary Figure S4C ). To characterize intercellular interactions in high- and low-risk groups, we inferred putative cell-to-cell interactions based on ligand–receptor signaling using CellChat (44). Interestingly, we observed globally enhanced intercellular interactions for the low-risk group ( Supplementary Figures S4D , S4E ). In addition, we found that CD46 signaling network was strengthened in the low-risk group, which suggested that it may play an important role in the process of inhibiting tumor progression ( Supplementary Figure S4F ). These results confirmed that aging-related genes play a crucial role in affecting the tumor microenvironment.

Overall, aging-related genes have good efficacy in predicting the prognosis of patients with LIHC, and they can serve as regulators in the tumor microenvironment.

In order to explore hub aging-related genes associated with liver cancer disease progression, we first screened aging-related genes that were significantly associated with survival in the model. The results from the K-M analysis suggested that CCNA2, EEF1E1, HDAC2, IL7R, and PON1 were significantly associated with OS, progression-free survival, and disease-free survival (P < 0.05, Figure 8A ). To define the correlation between five aging-related genes and clinical feature, Kruskal–Wallis test was involved. The result suggested that EEF1E1 and HDAC2 were associated with the progression of tumors where they show significant correlation with grade, stage, and T stage (P < 0.05; Figures 8B, C ), and the expression levels of EEF1E1 and HDAC2 were significantly higher in tumor tissue than in normal tissue ( Figure 8D ). In addition, using qPCR, we found that the mRNA expression of EEF1E1 and HDAC2 in hepatocellular carcinoma cells (HepG2, Huh6, and SNU398) are higher than those in normal hepatocytes (LO2) ( Figure 8H , Supplementary Table S2 ).

Using the HPA database, we found that the protein expression level of EEF1E1 and HDAC2 were higher in liver cancer while lower in normal liver ( Figures 8E–G ). By using the chi-square test, we found a significant difference in the protein expression level of normal tissue samples and hepatocellular carcinoma samples in the HPA database ( Figure 8G ). Using qPCR and Western blotting, we found that the mRNA and protein levels of HDAC2 in hepatocellular carcinoma cells (HepG2, Huh6, and SNU398) are higher than those in normal hepatocytes (LO2) ( Figures 8H, I ). We also found an mRNA level increase in EEF1A1 expression in LIHC cancer cells ( Figure 8H ). Interestingly, we also found that HDAC2 levels in secondary liver cancer tissues are higher than those in normal liver tissues, which further supports that HDAC2 plays an important role in the progression of liver cancer ( Figures 8J, K ).

To explore the biological function of these aging-related genes (EEF1E1 and HDAC2) in the progression of LIHC, we performed GSEA based on TCGA cohort. Enrichment result indicated that the PI3K/AKT/mTOR signaling pathway can be activated by EEF1E1 and HDAC2 (P < 0.05, Figure 8L ), which is consistent with previous results. The PI3K/AKT/mTOR pathway has been proven that it is aberrantly hyperactivated in many types of cancer that has a strong relationship with poor clinical prognosis. In the tumor microenvironment, the PI3K/AKT/mTOR pathway plays an important role in promoting the proliferation and metastasis of tumor while strongly inhibiting the antitumor immune response (45).

Interestingly, this result indicated that the PI3K/AKT/mTOR pathway may be a potential pathway that aging-related genes converge on to liver tumor progression.