Abstract
Background: Aging is a well-studied concept, but no studies have comprehensively analyzed the association between aging-related genes (AGs) and hepatocellular carcinoma (HCC) prognosis.
Methods: Gene candidates were selected from differentially expressed genes and prognostic genes in The Cancer Genome Atlas (TCGA) database. A gene risk score for overall survival prediction was established using the least absolute shrinkage and selection operator (LASSO) regression analysis, and this was validated using data from the International Cancer Genome Consortium (ICGC) database. Functional analysis was conducted using gene ontology enrichment, Kyoto Encyclopedia of Genes and Genomes analysis, gene set enrichment analysis, and immune microenvironment and tumor stemness analyses.
Results: Initially, 72 AGs from the TCGA database were screened as differentially expressed between normal and tumor tissues and as genes associated with HCC prognosis. Then, seven AGs (POLA1, CDK1, SOCS2, HDAC1, MAPT, RAE1, and EEF1E1) were identified using the LASSO regression analysis. The seven AGs were used to develop a risk score in the training set, and the risk was validated to have a significant prognostic value in the ICGC set (p < 0.05). Patients with high risk scores had lower tumor differentiation, higher stage, and worse prognosis (all p < 0.05). Multivariate Cox regression analyses also confirmed that the risk score was an independent prognostic factor for HCC in both the TCGA and ICGC sets (all p < 0.05). Further analysis showed that a high risk score was correlated with the downregulation of metabolism and tumor immunity.
Conclusion: The risk score predicts HCC prognosis and could thus be used as a biomarker not only for predicting HCC prognosis but also for deciding on treatment.
Introduction
The morbidity of hepatocellular carcinoma (HCC) is steadily decreasing worldwide (Fitzmaurice et al., ), but there are still ~900,000 newly diagnosed patients per year. According to the 2020 estimated worldwide age-standardized incidence rates, more than half of liver cancer patients are aging: the percentage of patients aged ≥60 years was 51.4% per 100,000, that of patients aged ≥65 years was 58.9% per 100,000, and that of patients aged ≥70 years was 65.5% per 100,000 (GLOBOCAN, ). In light of the continuous population growth and aging globally, there would be a greater number of aging HCC patients, which is concerning; however, the cut-off value for advanced age remains controversial (Kim et al., ; Oishi et al., ; Nozawa et al., ).
In addition, the association between aging and cancer is still unclear. In general, genome alterations such as DNA mutations accumulate with age (Vijg and Suh, ). However, recent studies have found that up to 50% of somatic mutations occur in early life (Moskalev et al., ). In addition, aged cells typically present with short telomeres, which indicate that cancers are less likely to occur during aging. However, latest evidence suggests that cancers could also arise from cells with short telomeres, regardless of the underlying mechanism of telomerase activation (Bernardes de Jesus and Blasco, ; de Magalhães, ). As an important immune and metabolic organ, the liver is also greatly affected by aging in terms of its function and structure, which involves regulation of the tumorigenicity and development of HCC via liver microenvironment modulation (Bernardes de Jesus and Blasco, ; Fan et al., ).
Interestingly, a latest retrospective study of 131 HCC patients (98 elderly patients and 33 children) reported that patients aged >85 years had smaller tumors compared with children, and that their survival rate was 92.4%, which was comparable with that of patients aged ≤ 85 years and was significantly higher than that of the children (Atyah et al., ). AKR1B10, which is involved in leading theories on aging, has also been identified as a crucial gene in the increase of carcinogenesis with age (Matkowskyj et al., ; Atyah et al., ). However, to our best knowledge, the role of the aging process in HCC remains unclear, and the correlation between aging genes and HCC prognosis has not yet been systematically evaluated. Thus, this study aimed to verify the association between aging-related genes (AGs) and HCC prognosis using data from the human aging genome resource (HAGR) (Tacutu et al., ), and then establish an aging-specific risk score to predict the prognosis of HCC patients using data from The Cancer Genome Atlas (TCGA) and the International Cancer Genome Consortium (ICGC).
Materials and Methods
Data Collection
A total of 307 human AGs were downloaded from HAGR (http://genomics.senescence.info/genes/, Supplementary Table 1) (Tacutu et al., ). The RNA sequencing (RNA-Seq) expression profile dataset of 371 HCC patients, which included data on clinicopathological characteristics and survival, was downloaded from TCGA (https://portal.gdc.cancer.gov/). Corresponding data from another 260 HCC patients were downloaded from the ICGC portal (https://dcc.icgc.org/projects/LIRI-JP). Considering that data from TCGA and ICGC are publicly available, the present study did not require approval of local ethics committees. However, this study was conducted according to the guidelines of TCGA and ICGC.
Construction and Valisdation of a Prognostic AG Signature
The gene expression profiles were standardized with the scale method using the “limma” R package. Prognostic AGs were screened using univariate Cox analysis of overall survival (OS). Differentially expressed genes (DEGs) between tumor tissues and adjacent non-tumorous tissues were distinguished using the “limma” R package, with a false discovery rate (FDR) < 0.05 in the TCGA set. Volcano plots and DEG heatmap analyses were depicted using the “ggplot2” and “pheatmap” R packages, respectively.
Construction of Protein–Protein Interaction Networks and Gene-Related Networks
A protein–protein interaction (PPI) network for overlapping prognostic DEGs was created using the STRING database (https://string-db.org/) and was visualized using Cytoscape (version 3.7.2). Meanwhile, based on TCGA gene expression data, an interaction gene-related network for overlapping prognostic DEGs was constructed using the R package “corrplot.” The network was visualized using the R package “circlize.”
Construction of the LASSO Cox Regression Model
Using the “glmnet” R package, a risk score was developed using the least absolute shrinkage and selection operator (LASSO) Cox regression analysis to predict the prognosis of HCC, which was designed for variable selection and shrinkage (Tibshirani, ; Simon et al., ). The standardized expression matrix of candidate prognostic DEGs was set as the independent variable in the regression, and the response variables were OS and patient status in the TCGA cohort. The risk score of each patient was determined based on the standardized expression level of each gene and its corresponding regression coefficients. The formula was established as follows:
The patients were then divided into the low-risk and high-risk groups according to the median value of the risk score.
Validation of the Risk Score
Principal component analysis (PCA) was performed using the “prcomp” function of the “stats” R package according to the gene expression in the signature. t-SNE was performed to verify the distribution of each group using the “Rtsne” R package. The optimal cut-off expression values for the survival analysis of different genes were determined using the “surv_cutpoint” function of the “survminer” R package. Receiver operating characteristic (ROC) curve analyses were performed using the “survivalROC” R package to evaluate the predictive power of the risk score and compare it with that of published risk scores.
Gene Set Enrichment Analysis and Functional Enrichment Analysis
Gene set enrichment analysis (GSEA) was performed using the GSEA software (v4.0.3) (https://www.gsea-msigdb.org/gsea/downloads.jsp) and the Molecular Signatures Database to generate a list of significantly different gene sets between the low- and high-risk groups (Subramanian et al., ). Gene sets with a p < 0.05 and FDR < 0.25 were considered significantly enriched. Permutations were set as 1,000.
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were carried out using the “clusterProfiler” R package according to the DEGs (|log2FC| ≥1, FDR < 0.05) between the low- and high-risk groups (Yu et al., ).
Immune-, Tumor Microenvironment-, and Tumor Stemness-Related Analyses
For immune-related analysis, 16 immune cells related infiltrating score and 13 immune-related pathways were determined using single-sample gene set enrichment analysis (ssGSEA) in the “gsva” R package (Rooney et al., ). To analyze the tumor microenvironment, we used the TCGA gene expression data. The immune and stromal ESTIMATE scores were determined using the “estimate” R package to calculate immune and stromal scores in tumor samples (Yoshihara et al., ). Meanwhile, tumor stemness scores were obtained from the Xena browser (https://xenabrowser.net/datapages/). This analysis was performed to compare tumor stemness between the low- and high-risk groups according to the TCGA gene expression data.
Development of Risk Prediction Model
A nomogram incorporating risk scores with clinical information was established to predict 1-, 2-, and 3-year OS for HCC patients using the R package “rms.” The nomogram's discrimination was evaluated using the concordance index (C-index). Age, gender, tumor–node–metastasis (TNM) staging, risk score, and the combined model were also compared using the C-index. The C-index and area under the curve (AUC) were used to compare the discriminative capability of the nomogram with that of other models. Similar comparisons were conducted to evaluate the calibration of the nomogram and measure the clinical utility of the nomogram; these comparisons were performed using the calibration plot and the decision curve analysis (DCA), respectively.
Statistical Analyses
Gene expression between tumor tissues and adjacent non-tumorous tissues was compared using Student's t-test, and differences in proportions were compared using the χ2-test. Kaplan–Meier analysis was used to compare the OS between different groups through the log-rank test, and multivariate Cox regression analyses were used to identify independent risk factors of OS. The diagnostic value of the risk score and nomogram model were evaluated using the ROC curves. All statistical analyses were conducted using R software (version 4.0.3). P with two-tailed value < 0.05 was considered statistically significant.
Results
Identification of Interaction AGs Between Differential Expression and Prognosis
The clinical characteristics of the patients in the TCGA and ICGC sets are depicted in Table 1. A total of 130 AGs were found to be associated with HCC prognosis using univariate Cox regression analysis (all p < 0.05, Figure 1A), and 116 DEGs between tumors and peritumors or normal tissues were identified (Figure 1B). Then, 72 AGs were identified as interaction genes between differential expression and prognosis (Figure 1C), among which 8 and 64 were found to be upregulated and downregulated in tumor tissues, respectively (Figure 1D).
Table 1
| TCGA cohort | LIRI-JP cohort | |
|---|---|---|
| No. of patients | 371 | 260 |
| Age | ||
| Median (range) | 61 (16–90) | 69 (31–89) |
| ≤ 65 (%) | 230 (62.0%) | 98 (37.7%) |
| >65 (%) | 141 (38.0%) | 162 (62.3%) |
| Gender (%) | ||
| Female | 120 (32.3%) | 68 (26.2%) |
| Male | 251 (67.7%) | 192 (73.8%) |
| Grade (%) | ||
| Grade 1 | 55 (14.8%) | NA |
| Grade 2 | 178 (50.0%) | NA |
| Grade 3 | 120 (32.3%) | NA |
| Grade 4 | 13 (3.5%) | NA |
| Unknown | 5 (1.4%) | NA |
| Stage (%) | ||
| I | 174 (46.9%) | 40 (15.4%) |
| II | 85 (22.9%) | 117 (45.0%) |
| III | 84 (22.6%) | 80 (30.8%) |
| IV | 4 (1.1%) | 23 (8.8%) |
| Unknown | 24 (6.5%) | 0 (0.0%) |
| T (%) | ||
| T1 | 185 (49.1%) | NA |
| T2 | 95 (25.2%) | NA |
| T3 | 81 (21.5%) | NA |
| T4 | 13 (3.4%) | NA |
| Tx | 3 (0.8%) | NA |
| N (%) | ||
| N0 | 257 (68.2%) | NA |
| N1 | 4 (1.1%) | NA |
| Nx | 116 (30.7%) | NA |
| M (%) | ||
| M0 | 272 (72.1%) | NA |
| M1 | 4 (1.1%) | NA |
| Mx | 101 (26.8%) | NA |
Clinical characteristics of the HCC patients used in this study.
Figure 1
PPI Networks and Correlational Analysis of Interaction AGs
PPI networks were determined among the 72 interaction AGs using the STRING database (Supplementary Figure 1A), and the results showed that BRCA1, CDK1, PCNA, CCNA2, POLD1, UBE2I, MAPK3, RAD51, SHC1, and HRAS were the most connected hub genes (connected with more than 10 nodes). The association between these genes is depicted in Supplementary Figure 1B. Among these, BUB1B, TOP2A, CDK1, LMNB1, and RAD51 were the most common AGs.
Construction of Prognostic Risk Score in the TCGA Training Set
Through the LASSO Cox regression analysis, the seven AGs (POLA1, CDK1, SOCS2, HDAC1, MAPT, RAE1, and EEF1E1; Table 2, Supplementary Figure 2) were used to establish a risk score to predict the OS of HCC patients in the TCGA training set (Supplementary Figure 3). The risk score of each patient was determined based on the following formula: Risk score = POLA1 × 0.015 + CDK1 × 0.041 − SOCS2 × 0.130 + HDAC1 × 0.060 + MAPT × 0.069 + RAE1 × 0.001 + EEF1E1 × 0.053.
Table 2
| ENSG ID | Symbol | Location | Expression status | Coefficient |
|---|---|---|---|---|
| ENSG00000101868 | POLA1 | Chromosome X | Upregulated | 0.015002532 |
| ENSG00000170312 | CDK1 | Chromosome 10 | Upregulated | 0.040778844 |
| ENSG00000120833 | SOCS2 | Chromosome 12 | Downregulated | −0.129910361 |
| ENSG00000116478 | HDAC1 | Chromosome 1 | Upregulated | 0.059678695 |
| ENSG00000186868 | MAPT | Chromosome 17 | Upregulated | 0.069245304 |
| ENSG00000101146 | RAE1 | Chromosome 20 | Upregulated | 0.00079437 |
| ENSG00000124802 | EEF1E1 | Chromosome 6 | Upregulated | 0.053295294 |
The 7 genes associated with the risk model in HCC.
Patients from the TCGA training set were then divided into low and high risk. The risk plot distribution in the TCGA set is shown in Figure 2A. PCA and t-SNE analysis suggested that the patients were distributed in two apparent directions (Figures 2B,C). Survival analysis showed that the OS of the low-risk group was significantly better than that of the high-risk group (p < 0.001, Figure 2D). The time-dependent curves exhibited good predictive performance in terms of the risk score for 1-, 2-, and 3-year OS, with the AUCs of 0.777, 0.763, and 0.733, respectively (Figure 2E).
Figure 2
Verification of the Risk Score in the ICGC Validation Set
Patients from the ICGC set were used to validate the current risk score, and they were also categorized into the low- and high-risk groups according to the median value using the same formula as that in the TCGA set (Figure 3A). Likewise, PCA and t-SNE analysis indicated that the patients in the two subgroups were also distributed in two discrete directions (Figures 3B,C). Similarly, survival analysis showed that low-risk patients had a prolonged OS compared with high-risk patients (p < 0.001, Figure 3D). The AUCs of the seven AGs in predicting 1-, 2-, and 3-year OS were 0.655, 0.721, and 0.705, respectively (Figure 3E). Meanwhile, the AUCs of the current risk score in predicting 1-, 2-, and 3-year OS were higher than those of the risk score of Yan et al. (), Fang and Chen (), Liang et al. (), and Tang et al. () (Supplementary Figure 4).
Figure 3
Evaluation of Clinicopathological Characteristics Using the 7-AG Model
The strip chart of clinicopathological characteristics and the heat maps of the expression profiles of the seven AGs are exhibited in Figure 4A. The corresponding scatter diagrams determined by the Wilcoxon signed-rank test showed that T stage (Figure 4B), clinical stage (Figure 4C), differentiation grade (Figure 4D), and survival status (Figure 4E) were related to the risk score.
Figure 4
Univariate and Multivariate Cox Regression Analyses of HCC Prognosis
Using the univariate Cox regression analysis, stage (HR = 2.500, 95% CI = 1.721–3.632, p < 0.001, Figure 5A) and risk score (HR = 56.069, 95% CI = 17.923–175.404, p < 0.001, Figure 5A) were identified to be related to the OS of patients in the TCGA training set. In the multivariate Cox regression model, stage (HR = 2.172, 95% CI = 1.490–3.167, p < 0.001, Figure 5B) and risk score (HR = 46.043, 95% CI = 14.413–147.082, p < 0.001, Figure 5B) were likewise confirmed to be independent prognostic factors of OS. Similarly, the risk score was also confirmed to be an independent risk factor for OS in the ICGC validation set (both p < 0.001, Figures 5C,D).
Figure 5
Functional Analysis in the TCGA Set and ICGC Set
GO enrichment and KEGG analysis were conducted using DEGs between the low- and high-risk groups to explore the biological functions and pathways related to risk score. The top 10 enriched GO terms of biological process, cellular components, and molecular function for the DEGs are depicted as a scatter diagram in Figure 6A (all adjusted p < 0.05), and as expected, almost all of the enriched terms were related to mitosis. KEGG pathway analysis revealed that DEGs were enriched in several cell cycle, cellular senescence, and retinol metabolism (all adjusted p < 0.05, Figure 6B) processes. Similar results were also identified in the ICGC set (all adjusted p < 0.05, Figures 6C,D).
Figure 6
GSEA of Risk Score–Related Signaling Pathways
GSEA was conducted to sort out significantly enriched pathways between the high- and low-risk groups in the TCGA set. A total of 100 enriched pathways with significant differences between groups of low and high risk were sorted out (FDR < 0.25, p < 0.05), which are depicted in Supplementary Table 2. Moreover, 27 enriched pathways were found to be associated with metabolism, 18 of which were enriched in the low-risk group. Pyrimidine metabolism, purine metabolism, inositol phosphate metabolism, amino sugar and nucleotide sugar metabolism, and selenoamino acid metabolism were the top five enriched metabolism-related pathways in the high-risk group (Figure 7A), whereas tryptophan metabolism, fatty acid metabolism, retinol metabolism, glycine, serine, and threonine metabolism, and propanoate metabolism were the top five enriched pathways in the low-risk group (Figure 7A). Interestingly, signaling pathways of the B cell receptor, T cell receptor, and TGF-β were enriched in group of the low risk (Figure 7B).
Figure 7
Correlations Between the Risk Score and Tumor Immunity, Tumor Microenvironment, and Tumor Stemness in the TCGA Set
The enriched scores of diverse immune cell subpopulations, related functions, or pathways were quantified using ssGSEA to further explore the associations between risk score and immune status. Interestingly, significant differences were observed between the low- and high-risk groups in terms of activated dendritic cells, macrophages, mast cells, neutrophils, natural killer (NK) cells, helper T cells, tumor-infiltrated lymphocytes, antigen process cell co-inhibition, cytolytic activity, major histocompatibility complex class I, type I interferon (IFN) response, and type II IFN response (all adjusted p < 0.05, Figures 8A,B).
Figure 8
The ESTIMATE algorithm was used for determining the stromal and immune scores to investigate the association between the low- and high-risk groups. Surprisingly, apparent differences between the two risk groups were observed in terms of the ESTIMATE score (Figure 8C), but further analysis showed that the differences existed in terms of the stromal score rather than the immune score (Figure 8C).
Furthermore, mRNA expression and DNA methylation patterns were measured in the RNA stemness score (RNAss) and DNA stemness score (DNAss) to determine the association between the risk score and tumor stemness. Interestingly, RNAss was lower in the low-risk group than in the high-risk group (p < 0.001, Figure 8D), whereas there was no between-group difference in DNAss (Figure 8D).
Development of a Nomogram in the TCGA Set
Based on the survival analysis, a nomogram incorporating the risk score and other clinicopathological information was developed to predict the OS of patients in the TCGA set, as shown in Figure 9A. The C-index was 0.528 (0.467–0.589), 0.530 (0.425–0.635), 0.548 (0.461–0.634), 0.669 (0.595–0.742), 0.717 (0.670–0.764), and 0.730 (0.681–0.778) for age, gender, grade, stage, risk score, and nomogram, respectively (Figure 9B). Good consistency was identified between the observed outcomes and the predicted 1-, 2-, and 3-year OS outcome of the nomogram (Figures 9C–E). The AUC of the nomogram to predict the OS at 1, 2, and 3 years was higher than that of age, gender, grade, stage, and risk score (Figures 9F–H). DCA showed that the nomogram exhibited better net benefits than the age, gender, grade, stage, and risk score with respect to the 1-, 2-, and 3-year OS (Figures 9I–K).
Figure 9
Discussion
The role of aging in HCC, particularly its prognostic impact, remains unclear. In this study, a risk score incorporating seven AGs was established to predict the OS of HCC patients. The results of the functional analysis revealed that mitosis- and metabolism-related pathways were enriched, and the risk score was associated with tumor immunity. To our best knowledge, this is the first study to systematically investigate AGs in HCC tumor tissues and their clinical value in HCC prognosis.
Aging is associated with impaired immunological function and an accumulated chronic inflammatory microenvironment, which might facilitate tumor formation and progression (Kuilman et al., ). Recently, a retrospective study found that elderly patients exhibited better biological characteristics than children, but the underlying mechanisms for this observation remain unknown (Atyah et al., ). AKR1B10, which plays an important role in the adoptive response to oxidative stress and damage and is considered to be involved in leading theories of aging, was identified as a crucial gene in the increase of carcinogenesis with age (Schmitz et al., ; Sato et al., ). However, the number of available data and patients, especially in the younger group, was too limited to draw definite conclusions.
In the current study, interaction AGs of DEGs and prognosis were initially identified, and seven AGs were then used to construct a risk score through the LASSO Cox regression model, which showed good performance in predicting prognosis both in the training and validation sets. Compared with low-risk patients, high-risk patients typically present with more aggressive characteristics, advanced stages, and higher mortality (all p < 0.05). Furthermore, the risk score exhibited better predictive capability than the conventional TNM stage and published risk scores (Fang and Chen, ; Liang et al., ; Tang et al., ; Yan et al., ).
In the seven AGs model, POLA1, CDK1, HDAC1, MAPT, RAE1, and EEF1E1 acted as risk factors, whereas only SOCS2 was a protective factor. CDK1 was reportedly amplified in HCC tissues and is significantly associated with poor OS (Wu et al., ). CDK1 expression is positively correlated with tumor stemness indices, and growing evidence suggests that the key biological processes of tumors, such as invasion, metastasis, and therapeutic resistance to chemotherapy and radiotherapy, are largely dependent on tumor stem cells (Matthai and Ramakrishna, ). In the current study, a significant increase in RNAss was observed in the high-risk group compared with the low-risk group, which also demonstrated the adverse effect of high RNAss on the prognosis of HCC patients. As previously reported, POLA1 was found to be a risk factor for HCC, and the underlying mechanism for this was the suppression of the antitumor immune response through IFN-γ inhibition (Starokadomskyy et al., ). Suppression of SOCS2 was found to induce tumor immune escape by inhibiting the development of Th2 cells and restricting the adaptive anti-tumor immunity of T cells in multiple tumor models. Patients with SOCS2 overexpression have better prognosis, which is consistent with the current study (Knosp et al., ; Cramer et al., ). HDAC1 plays a role in cancer initiation and progression in many cancers by removing acetyl groups (Ropero and Esteller, ). HCC patients with high HDAC1 expression reportedly have a poor survival rate (Rikimaru et al., ). RAE1 is involved in the nuclear output of poly (A)+ RNA (Bharathi et al., ) and is negatively correlated with breast cancer prognosis, which is similar to the results found for HCC in the current study (Chin et al., ; Oh et al., ). EEF1E1 is a scaffold of the macromolecule aminoacyl tRNA synthase complex, which mediates ATM/ATR-mediated p53 activation as a protective factor of prognosis (Yu et al., ; Biterge-Sut, ). MAPT plays an important role in the assembly of tubulin and stabilization of microtubules, which act as protective factors in renal and prostate cancers (Han et al., ; Sekino et al., ). However, in the current study, EEF1E1 and MAPT were found to be risk factors for prognosis in the TCGA database, as previously reported (Biterge-Sut, ; Zhang et al., ), which deserves further validation. In summary, they were first identified to represent a prognostic risk model for HCC.
Tumor cells are typically present with the activation of cyclin and upregulation of mitosis. In the current study, all seven AGs were related to the cell cycle. Classical cell cycle–specific agents such as etoposide, irinotecan, and bleomycin have been used prevalently in various cancers (Atienza et al., ), and many novel drugs such as palbociclib and CFI-402257 are also related to the cell cycle (Thu et al., ). As a member of CDKs, CDK1 plays a crucial role in the cell cycle process (Xiao et al., ). Previous studies have shown that CDKs and their regulatory factors are abnormally activated in many tumors (Abbosh et al., ). CDK1 overexpression leads to abnormal cell proliferation, and CDK activity is required in response to DNA damage during DNA replication (Müllers et al., ). In the current study, CDK1 was identified as an important gene in the prognostic model of AGs, and the upregulated genes (DEGs) were significantly enriched in the cell cycle pathways, as per the GO and KEGG analyses. The KEGG pathway analysis revealed that the cell cycle pathway was significantly activated in upregulated DEGs. These findings demonstrated that the potential mechanism of AGs might lie in the cell cycle, which might serve as a potential target in the future.
As expected, AGs in the current study were associated with cellular senescence, according to the KEGG analysis. In a study of glioma, increased senescence scores are associated with increasing age and higher malignancy, as determined by WHO, and are somewhat associated with poor prognosis (Coppola et al., ). A previous study reported that senescent cells could regulate the microenvironment by secreting matrix metalloproteinase-2, and reshape the extracellular matrix through overexpression of related proteins, which in turn helps low-grade malignant epithelial cells to form more blood vessels and promote the invasion of melanoma cells (Mo et al., ). In addition, for senescent liver cells with short telomerases, the degree of telomerase shortening was positively correlated with the degree of DNA damage (Aravinthan et al., ). These findings suggest that cellular senescence is closely related to AGs and hence deserves further study.
Generally, cancer cells require core metabolic functions to generate energy, control redox, and assimilate biomass to support cellular overgrowth (Vander Heiden and DeBerardinis, ). Cyclins and cyclin-dependent kinases (CDKs) are often involved in the majority of metabolic processes, such as glucose metabolism (Lee et al., ; Wang et al., ), lipogenesis (Zhao et al., ), amino acid metabolism (Tarrado-Castellarnau et al., ), and mitochondrial activity (Hydbring et al., ). Previous studies reported that impaired liver function is correlated with worse prognosis in HCC patients (Granito and Bolondi, ). In the current study, most of the metabolism-related signaling pathways were enriched in the low-risk group, which indicates that liver function might be impaired. Hence, impaired liver function in high-risk patients might be the reason for their unfavorable prognosis.
Interestingly, in the current study, B cell receptor and T cell receptor signaling pathways were found to be enriched in the low-risk group as per GSEA, which suggested that they may be inhibited in the high-risk group. Previous studies found that decreased TILs (Xu et al., ), T helper cells (Yan et al., ), and NK cells (Cai et al., ) were associated with poor prognosis in HCC. In this study, the high-risk group was found to be contained lower fractions of TILs, T helper cells, and NK cells through ssGSEA. Furthermore, we also identified that a high-risk score had a lower score in terms of APC co-inhibition, cytolytic activity, and type I/II IFN response. Hence, impaired antitumor immunity in high-risk patients might be the reason for their unfavorable prognosis.
Nevertheless, there are several limitations in this study. First, the current risk score was both established and validated using retrospective data from public databases, which should be verified through prospective real-world data. Second, considering that HCC is a typical polygenic disease, the intrinsic disadvantage of merely considering AGs to build a prognostic model was unavoidable. The last but not the least, the correlations between the risk score and cell cycle, metabolism, and immune activity have not yet been confirmed in vitro and in vivo experiment.
In conclusion, our study developed a novel risk score using seven AGs, which proved to be independently associated with OS in both the derivation and validation cohorts. However, the correlations between AGs and tumor immunity in HCC warrant further investigation.
Statements
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Materials.
Author contributions
Study concept and design: XTC and LW. Data analysis/interpretation: XTC and LH. Manuscript drafting: LW and LH. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by grants from Fujian Province Finance Department Project [No. (2019) 827].
Acknowledgments
We sincerely thank the researchers and study participants for their contributions toward this study.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.661988/full#supplementary-material
References
1
AbboshC.SwantonC.BirkbakN. (2018). Circulating tumour DNA analyses reveal novel resistance mechanisms to CDK inhibition in metastatic breast cancer. Ann. Oncol.29, 535–537. 10.1093/annonc/mdy017
2
AravinthanA.ScarpiniC.TachtatzisP.VermaS.Penrhyn-LoweS.HarveyR.et al. (2013). Hepatocyte senescence predicts progression in non-alcohol-related fatty liver disease. J. Hepatol.58, 549–556. 10.1016/j.jhep.2012.10.031
3
AtienzaJ.RothR.RosetteC.SmylieK.KammererS.RehbockJ.et al. (2005). Suppression of RAD21 gene expression decreases cell growth and enhances cytotoxicity of etoposide and bleomycin in human breast cancer cells. Mol. Cancer Ther.4, 361–368. 10.1158/1535-7163.MCT-04-0241
4
AtyahM.YinY.ZhouC.ZhouQ.ChenW.DongQ.et al. (2018). Integrated analysis of the impact of age on genetic and clinical aspects of hepatocellular carcinoma. Aging10, 2079–2097. 10.18632/aging.101531
5
Bernardes de JesusB.BlascoM. (2013). Telomerase at the intersection of cancer and aging. Trends Genet.29, 513–520. 10.1016/j.tig.2013.06.007
6
BharathiA.GhoshA.WhalenW.YoonJ.PuR.DassoM.et al. (1997). The human RAE1 gene is a functional homologue of Schizosaccharomyces pombe rae1 gene involved in nuclear export of Poly(A)+ RNA. Gene. 198, 251–258. 10.1016/S0378-1119(97)00322-3
7
Biterge-SutB. (2019). Alterations in eukaryotic elongation factor complex proteins (EEF1s) in cancer and their implications in epigenetic regulation. Life Sci.238:116977. 10.1016/j.lfs.2019.116977
8
CaiL.ZhangZ.ZhouL.WangH.FuJ.ZhangS.et al. (2008). Functional impairment in circulating and intrahepatic NK cells and relative mechanism in hepatocellular carcinoma patients. Clin. Immunol.129, 428–437. 10.1016/j.clim.2008.08.012
9
ChinK.DeVriesS.FridlyandJ.SpellmanP.RoydasguptaR.KuoW.et al. (2006). Genomic and transcriptional aberrations linked to breast cancer pathophysiologies. Cancer Cell10, 529–541. 10.1016/j.ccr.2006.10.009
10
CoppolaD.BalducciL.ChenD.-T.LobodaA.NebozhynM.StallerA.et al. (2014). Senescence-associated-gene signature identifies genes linked to age, prognosis, and progression of human gliomas. J. Geriatric Oncol.5, 389–399. 10.1016/j.jgo.2014.08.003
11
CramerA.de Lima OliveiraB. C.LeiteP. G.RodriguesD. H.BrantF.EsperL.et al. (2019). Role of SOCS2 in the regulation of immune response and development of the experimental autoimmune encephalomyelitis. Mediators Inflamm.2019:1872593. 10.1155/2019/1872593
12
de MagalhãesJ. (2013). How ageing processes influence cancer. Nat. Rev. Cancer13, 357–365. 10.1038/nrc3497
13
FanC.JinG.SunC.MaD.ChenC.QiaoP.et al. (2014). Clinical characteristics, treatment patterns and survival outcome of hepatocellular carcinoma patients aged 70 years or older: a single-center retrospective study from China. Aging Clin. Exp. Res.26, 123–130. 10.1007/s40520-013-0142-6
14
FangQ.ChenH. (2020). Development of a novel autophagy-related prognostic signature and nomogram for hepatocellular carcinoma. Front. Oncol.10:591356. 10.3389/fonc.2020.591356
15
FitzmauriceC.AllenC.BarberR.BarregardL.BhuttaZ.BrennerH.et al. (2017). Global, regional, and national cancer incidence, mortality, years of life lost, years lived with disability, and disability-adjusted life-years for 32 cancer groups, 1990 to 2015: a systematic analysis for the global burden of disease study. JAMA Oncol.3, 524–548. 10.1001/jamaoncol.2016.5688
16
GLOBOCAN (2020). Global Cancer Observatory. Available online at: https://gco.iarc.fr/today/ (accessed January 20, 2021).
17
GranitoA.BolondiL. (2017). Non-transplant therapies for patients with hepatocellular carcinoma and Child-Pugh-Turcotte class B cirrhosis. Lancet Oncol.18, e101–e112. 10.1016/S1470-2045(16)30569-1
18
HanX.SekinoY.BabasakiT.GotoK.InoueS.HayashiT.et al. (2020). Microtubule-associated protein tau (MAPT) is a promising independent prognostic marker and tumor suppressive protein in clear cell renal cell carcinoma. Urol. Oncol.38, 605.e609–e617. 10.1016/j.urolonc.2020.02.010
19
HydbringP.MalumbresM.SicinskiP. (2016). Non-canonical functions of cell cycle cyclins and cyclin-dependent kinases. Nat. Rev. Mol. Cell Biol.17, 280–292. 10.1038/nrm.2016.27
20
KimY.JangB.KimE.ChungW.ParkK.ChoK.et al. (2012). Hepatocellular carcinoma in the elderly: clinical characteristics, treatment, survival analysis in Korean patients older than 70 years. J. Korean Med. Sci.27, 1147–1154. 10.3346/jkms.2012.27.10.1147
21
KnospC. A.SchieringC.SpenceS.CarrollH. P.NelH. J.OsbournM.et al. (2013). Regulation of Foxp3+ inducible regulatory T cell stability by SOCS2. J. Immunol.190, 3235–3245. 10.4049/jimmunol.1201396
22
KuilmanT.MichaloglouC.VredeveldL.DoumaS.van DoornR.DesmetC.et al. (2008). Oncogene-induced senescence relayed by an interleukin-dependent inflammatory network. Cell133, 1019–1031. 10.1016/j.cell.2008.03.039
23
LeeY.DominyJ.ChoiY.JurczakM.TollidayN.CamporezJ.et al. (2014). Cyclin D1-Cdk4 controls glucose metabolism independently of cell cycle progression. Nature510, 547–551. 10.1038/nature13267
24
LiangJ.-y.WangD.-s.LinH.-c.ChenX.-x.YangH.ZhengY.et al. (2020). A novel ferroptosis-related gene signature for overall survival prediction in patients with hepatocellular carcinoma. Int. J. Biol. Sci.16, 2430–2441. 10.7150/ijbs.45050
25
MatkowskyjK.BaiH.LiaoJ.ZhangW.LiH.RaoS.et al. (2014). Aldoketoreductase family 1B10 (AKR1B10) as a biomarker to distinguish hepatocellular carcinoma from benign liver lesions. Human Pathol.45, 834–843. 10.1016/j.humpath.2013.12.002
26
MatthaiS.RamakrishnaB. (2015). Cancer stem cells in hepatocellular carcinoma–an immunohistochemical study with histopathological association. Indian J. Med. Res.142, 391–398. 10.4103/0971-5916.169195
27
MoJ.SunB.ZhaoX.GuQ.DongX.LiuZ.et al. (2013). Hypoxia-induced senescence contributes to the regulation of microenvironment in melanomas. Pathol. Res. Pract.209, 640–647. 10.1016/j.prp.2013.07.004
28
MoskalevA.ShaposhnikovM.PlyusninaE.ZhavoronkovA.BudovskyA.YanaiH.et al. (2013). The role of DNA damage and repair in aging through the prism of Koch-like criteria. Ageing Res. Rev.12, 661–684. 10.1016/j.arr.2012.02.001
29
MüllersE.Silva CascalesH.BurdovaK.MacurekL.LindqvistA. (2017). Residual Cdk1/2 activity after DNA damage promotes senescence. Aging Cell16, 575–584. 10.1111/acel.12588
30
NozawaA.KuboS.TakemuraS.SakataC.UrataY.NishiokaT.et al. (2015). Hepatic resection for hepatocellular carcinoma in super-elderly patients aged 80 years and older in the first decade of the 21st century. Surg. Today45, 851–857. 10.1007/s00595-014-0994-1
31
OhJ.HurH.LeeJ.KimY.SeoY.KimM. (2017). The mitotic checkpoint regulator RAE1 induces aggressive breast cancer cell phenotypes by mediating epithelial-mesenchymal transition. Sci. Rep.7:42256. 10.1038/srep42256
32
OishiK.ItamotoT.KohashiT.MatsuguY.NakaharaH.KitamotoM. (2014). Safety of hepatectomy for elderly patients with hepatocellular carcinoma. World J. Gastroenterol.20, 15028–15036. 10.3748/wjg.v20.i41.15028
33
RikimaruT.TaketomiA.YamashitaY.ShirabeK.HamatsuT.ShimadaM.et al. (2007). Clinical significance of histone deacetylase 1 expression in patients with hepatocellular carcinoma. Oncology72, 69–74. 10.1159/000111106
34
RooneyM.ShuklaS.WuC.GetzG.HacohenN. (2015). Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell160, 48–61. 10.1016/j.cell.2014.12.033
35
RoperoS.EstellerM. (2007). The role of histone deacetylases (HDACs) in human cancer. Mol. Oncol.1, 19–25. 10.1016/j.molonc.2007.01.001
36
SatoS.GendaT.IchidaT.MurataA.TsuzuraH.NaritaY.et al. (2016). Impact of aldo-keto reductase family 1 member B10 on the risk of hepatitis C virus-related hepatocellular carcinoma. J. Gastroenterol. Hepatol.31, 1315–1322. 10.1111/jgh.13295
37
SchmitzK.SotiropoulosG.BabaH.SchmidK.MüllerD.PaulA.et al. (2011). AKR1B10 expression is associated with less aggressive hepatocellular carcinoma: a clinicopathological study of 168 cases. Liver Int.31, 810–816. 10.1111/j.1478-3231.2011.02511.x
38
SekinoY.HanX.BabasakiT.GotoK.InoueS.HayashiT.et al. (2020). Microtubule-associated protein tau (MAPT) promotes bicalutamide resistance and is associated with survival in prostate cancer. Urol. Oncol.38, 795.e791–795 e798. 10.1016/j.urolonc.2020.04.032
39
SimonN.FriedmanJ.HastieT.TibshiraniR. (2011). Regularization paths for Cox's proportional hazards model via coordinate descent. J. Stat. Softw.39, 1–13. 10.18637/jss.v039.i05
40
StarokadomskyyP.GemelliT.RiosJ. J.XingC.WangR. C.LiH.et al. (2016). DNA polymerase-alpha regulates the activation of type I interferons through cytosolic RNA:DNA synthesis. Nat. Immunol.17, 495–504. 10.1038/ni.3409
41
SubramanianA.TamayoP.MoothaV.MukherjeeS.EbertB.GilletteM.et al. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U.S.A.102, 15545–15550. 10.1073/pnas.0506580102
42
TacutuR.ThorntonD.JohnsonE.BudovskyA.BarardoD.CraigT.et al. (2018). Human Ageing Genomic Resources: new and updated databases. Nucleic Acids Res.46, D1083–D1090. 10.1093/nar/gkx1042
43
TangB.ZhuJ.LiJ.FanK.GaoY.ChengS.et al. (2020). The ferroptosis and iron-metabolism signature robustly predicts clinical diagnosis, prognosis and immune microenvironment for hepatocellular carcinoma. Cell Commun. Signal.18:174. 10.1186/s12964-020-00663-1
44
Tarrado-CastellarnauM.de AtauriP.Tarragó-CeladaJ.PerarnauJ.YunevaM.ThomsonT.et al. (2017). De novo MYC addiction as an adaptive response of cancer cells to CDK4/6 inhibition. Mol. Syst. Biol.13:940. 10.15252/msb.20167321
45
ThuK. L.Soria-BretonesI.MakT. W.CesconD. W. (2018). Targeting the cell cycle in breast cancer: towards the next phase. Cell Cycle17, 1871–1885. 10.1080/15384101.2018.1502567
46
TibshiraniR. (1997). The lasso method for variable selection in the Cox model. Stat. Med.16, 385–39510.1002/(SICI)1097-0258(19970228)16:4<385::AID-SIM380>3.0.CO;2-3
47
Vander HeidenM.DeBerardinisR. (2017). Understanding the intersections between metabolism and cancer biology. Cell168, 657–669. 10.1016/j.cell.2016.12.039
48
VijgJ.SuhY. (2013). Genome instability and aging. Annu. Rev. Physiol.75, 645–668. 10.1146/annurev-physiol-030212-183715
49
WangH.NicolayB.ChickJ.GaoX.GengY.RenH.et al. (2017). The metabolic function of cyclin D3-CDK6 kinase in cancer cell survival. Nature546, 426–430. 10.1038/nature22797
50
WuC. X.WangX. Q.ChokS. H.ManK.TsangS. H. Y.ChanA. C. Y.et al. (2018). Blocking CDK1/PDK1/beta-Catenin signaling by CDK1 inhibitor RO3306 increased the efficacy of sorafenib treatment by targeting cancer stem cells in a preclinical model of hepatocellular carcinoma. Theranostics8, 3737–3750. 10.7150/thno.25487
51
XiaoY.LucasB.MolchoE.SchiffT.VigodnerM. (2016). Inhibition of CDK1 activity by sumoylation. Biochem. Biophys. Res. Commun.478, 919–923. 10.1016/j.bbrc.2016.08.051
52
XuX.TanY.QianY.XueW.WangY.DuJ.et al. (2019). Clinicopathologic and prognostic significance of tumor-infiltrating CD8+ T cells in patients with hepatocellular carcinoma: a meta-analysis. Medicine98:e13923. 10.1097/MD.0000000000013923
53
YanH.ChenY.WangK.YuL.HuangX.LiQ.et al. (2020). Identification of immune landscape signatures associated with clinical and prognostic features of hepatocellular carcinoma. Aging12, 19641–19659. 10.18632/aging.103977
54
YanJ.LiuX.XiaoG.LiN.DengY.HanL.et al. (2014). Prevalence and clinical relevance of T-helper cells, Th17 and Th1, in hepatitis B virus-related hepatocellular carcinoma. PLoS ONE9:e96080. 10.1371/journal.pone.0096080
55
YoshiharaK.ShahmoradgoliM.MartínezE.VegesnaR.KimH.Torres-GarciaW.et al. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun.4:2612. 10.1038/ncomms3612
56
YuG.WangL.HanY.HeQ. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. Omics16, 284–287. 10.1089/omi.2011.0118
57
YuX.ZhengH.ChanM. T.WuW. K. (2017). HULC: an oncogenic long non-coding RNA in human cancer. J. Cell. Mol. Med.21, 410–417. 10.1111/jcmm.12956
58
ZhangZ.LiJ.HeT.OuyangY.HuangY.LiuQ.et al. (2019). Two predictive precision medicine tools for hepatocellular carcinoma. Cancer Cell Int.19:290. 10.1186/s12935-019-1002-z
59
ZhaoX.FengD.WangQ.AbdullaA.XieX.ZhouJ.et al. (2012). Regulation of lipogenesis by cyclin-dependent kinase 8-mediated control of SREBP-1. J. Clin. Investig.122, 2417–2427. 10.1172/JCI61462
Summary
Keywords
aging, gene, hepatocellular carcinoma, tumor immunity, metabolism, prognosis, risk score
Citation
Chen X, Wang L, Hong L, Su Z, Zhong X, Zhou H, Zhang X, Wu J and Shao L (2021) Identification of Aging-Related Genes Associated With Clinical and Prognostic Features of Hepatocellular Carcinoma. Front. Genet. 12:661988. doi: 10.3389/fgene.2021.661988
Received
01 February 2021
Accepted
04 May 2021
Published
23 June 2021
Volume
12 - 2021
Edited by
Soichiro Yamamura, University of California, San Francisco, United States
Reviewed by
Jun Liu, Yuebei People's Hospital, China; Yun Hak Kim, Pusan National University, South Korea
Updates
Copyright
© 2021 Chen, Wang, Hong, Su, Zhong, Zhou, Zhang, Wu and Shao.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Junxin Wu junxinwufj@aliyun.comLingdong Shao 1065900966@qq.com
This article was submitted to Genetics of Aging, a section of the journal Frontiers in Genetics
†These authors have contributed equally to this work
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.