A novel telomere-associated genes signature for the prediction of prognosis and treatment responsiveness of hepatocellular carcinoma

Abstract

Background

Hepatocellular carcinoma (HCC) is a prevalent malignancy worldwide, characterized by its high malignancy and poor prognosis. Telomeres, crucial components of eukaryotic chromosomes, have been increasingly recognized for their involvement in tumorigenesis, development, and impact on the prognosis of cancer patients. However, the precise role of telomere-associated genes in HCC remains incompletely elucidated.

Methods

The Cancer Genome Atlas (TCGA) database was utilized to download data from 374 HCC and 50 normal liver tissue samples. Differential genes were screened and intersected with 2093 telomere-related genes (TRGs) in GeneCards, resulting in the identification of 704 TRGs exhibiting survival differences. Through univariate Cox regression analysis, multivariate Cox regression analysis, and LASSO regression, a prognostic model consisting of 18 TRGs for HCC risk assessment was developed. The single-cell and spatial transcriptomics were utilized to analyze the expression and distribution of 18 TRGs in HCC. Subsequently, Mendelian randomization (MR) analysis confirmed a causal relationship between ASF1A and alcoholic HCC among the identified 18 TRGs. The expression and functional significance of ASF1A in HCC cell lines were investigated through colony formation assays, Transwell migration assays, and wound healing experiments.

Results

We developed a prognostic risk model for HCC incorporating 18 TRGs. Kaplan–Meier analysis demonstrated that the overall survival (OS) rate of the high-risk group was significantly inferior to that of the low-risk group. Cox regression analysis identified age (HR = 1.017, 95% CI: 1.002–1.032, P = 0.03), stage (HR = 1.389, 95% CI: 1.111–1.737, P = 0.004), and risk score (HR = 5.097, 95% CI: 3.273–7.936, P < 0.001) as three independent risk factors for HCC patients. The five-year receiver operating characteristic curve (ROC) and multivariate Cox regression analysis further validated the accuracy of our model. Time-dependent ROC results revealed that the 1-year, 3-year, and 5-year AUC values were AUC = 0.801, AUC = 0.734, and AUC = 0.690, respectively. The expression and distribution of 18 TRGs in HCC were further validated through single-cell and spatial transcriptomics data. Additionally, immune subtype analysis indicated a significantly lower proportion of C3 and C4 subtypes in the high-risk TRG group compared to the low-risk group. Meanwhile, tumor immune dysfunction and exclusion (TIDE) were significantly higher in the high-risk group than in the low-risk group. Furthermore, we observed differences in IC50 values among nine chemotherapeutic drugs across different TRG risk subtypes which partially confirmed our model’s predictive efficacy for immunotherapy. Amongst these eighteen TRGs analyzed by MR analysis, ASF1A was found to be associated with alcoholic HCC pathogenesis. We further confirmed ASF1A was significant overexpression in HCC by Western blotting. We also explored it’s the carcinogenic role of ASF1A in HCC via the transwell, wound healing, and clone formation experiments.

Conclusion

In this study, we developed a novel prognostic model comprising 18 TRGs for HCC, which exhibited remarkable accuracy in predicting HCC patients’ prognosis. Additionally, through MR analysis, we have successfully established a causal relationship between ASF1A and alcoholic HCC for the first time, which also provided a new theoretical foundation for the management of alcoholic HCC.

Introduction

Hepatocellular carcinoma (HCC), which accounts for more than 75% of incidence of liver cancer worldwide, has received the majority of attention in the global public health concern [1, 2]. The pathogenesis of HCC is intricate and involves significant risk factors such as hepatitis viruses and alcohol consumption [3, 4]. Moreover, aberrant expression or mutations in various genes have been closely associated with the initiation and progression of HCC. However, studies investigating telomere-related genes (TRGs) remain relatively scarce.

Telomeres, repetitive sequences of nucleotides (TTAGGG), are located at the ends of eukaryotic chromosomes and play a crucial role in determining cell lifespan and division capacity. In normal human cells, telomere length gradually decreases with each cell division [5]. Telomere shortening is directly associated with various diseases, including tumorigenesis [6, 7]. However, the significance of telomere shortening in tumors remains controversial as it may inhibit tumor formation and progression by impeding cell proliferation but can also lead to widespread genomic instability, increasing the risk of tumorigenesis [8,9,10]. Studies have shown that variations in telomere structure were correlated with the development of ovarian, prostate, breast colorectal, and lung cancers [11]. Furthermore, shortened telomeres elevated the risk of breast and pancreatic cancer development [12, 13]. However, previous studies have primarily focused on the role of telomere length and stability in tumor prognosis, with a limited investigation into the impact of TRGs in HCC [14, 15].

In this study, we aimed to establish a groundbreaking risk model for telomere-related genes to predict prognosis in HCC, while also assessing its potential implications in tumor immune escape and therapeutic drug selection. Furthermore, through MR analysis, we provided compelling evidence that supporting a causal relationship between the ASF1A and alcoholic HCC. Additionally, we validated the expression of ASF1A in HCC and elucidated its pivotal role. Therefore, we identified novel prognostic markers for HCC, progressed the study of TRGs in HCC, and offered new insights into early diagnosis, treatment strategies optimization, as well as improved prognosis management of HCC.

Materials and methods

Data acquisition and screening of TRGs

The flowchart of this study was illustrated in Supplementary Fig. 1. We obtained data on 374 HCC patients and 50 normal liver tissue samples from the TCGA database [16], which included clinical information and raw transcriptome expression data (https://portal.gdc.cancer.gov/repository). From GeneCards (https://www.genecards.org/) [17], we retrieved a total of 2093 TRGs (Supplementary Table 1). By intersecting these TRGs with differentially expressed genes identified between tumor and normal groups in the TCGA dataset, we identified 949 TRGs that were differentially expressed in tumors. Finally, through univariate COX regression analysis, we further narrowed down to 704 TRGs that exhibited survival differences among HCC patients (Supplementary Table 2).

Construction and validation of the prognostic model for TRGs

Using the LASSO regression method, we applied a penalty coefficient proportional to the shrinkage of the regression coefficient to construct a prognostic risk model for HCC using 704 differentially expressed TRGs. This approach enabled us to identify TRGs with prognostic value. Based on this methodology, we developed an HCC risk prognostic model consisting of 18 TRGs (Supplementary Table 3). The TCGA training group and GSE14520 patients were categorized into low-risk and high-risk subgroups based on the median risk score. Kaplan–Meier survival analysis was conducted to compare OS between different risk subtypes. Additionally, univariate Cox analysis, multivariate Cox analysis, and ROC curve analysis were performed to evaluate the specificity and sensitivity of our model in predicting survival outcomes for HCC patients. To assess whether our prognostic model can serve as an independent predictor of OS in HCC patients, we utilized a nomogram. Finally, we represented the predicted value of the histogram as a ROC using the R package “timeROC”.

Acquisition and analysis of single-cell data

This study incorporated single-cell transcriptome data from GSE146115 (encompassing 16 samples from 4 liver cancer patients)and GSE179795 (consisting of 24 liver cancer samples), and analyzed these single-cell transcriptome data using the R package ‘Seurat (v4.3.0)’. To ensure the quality of the data, we performed data filtering using CreateSeuratObject: min.cells = 3, min. features = 200. Utilizing the “harmony” R package to address batch effects in cancer and adjacent tissue samples, followed by normalization of single-cell RNA sequencing data using the “Seurat” R package to identify highly variable genes and reduce data dimensions through principal component analysis. Utilize the “FindAllMarkers” function from the “scran” R package for identification of genes exhibiting specific differential expression, and manually annotate cell subclusters using the Cellmarker database. The Harmony algorithm was employed to correct batch effects, uniform manifold approximation and projection (UMAP) was utilized for reducing dimensionality, and KNN clustering was applied to choose cell subclusters. Subsequently, the AUCell algorithm assessed the scores of telomere-related genes across various annotated cell types in HCC.

Trajectory analysis

Following the utilization of the “Monocle3” R package for trajectory analysis, the study examined the trajectory patterns of particular cell clusters. This approach can demonstrate alterations in gene expression and precisely pinpoint genes that are under regulation throughout the complete duration of the trajectory.

Concurrent investigation of single-cell and spatial transcriptomics

We acquired spatial transcriptomic data for HCC from the Mendeley Data repository and conducted data retrieval, standardization, dimensionality reduction, cluster analysis, and identification of highly variable genes. SPOTlight(v1.6.7) integrates data from single-cell and spatial transcriptomics (ST) to discern the cellular composition and distribution within each designated area.

By utilizing the logNormCounts function, we standardized the gene expression data for individual cells and subsequently developed a genetic variation model to identify genes exhibiting high variability. Following this, we identified 3,000 genes with significant variation for further analysis. Subsequently, non-ribosomal and non-mitochondrial genes were filtered out and their marker scores were computed. After rigorous screening and sorting, only marker genes with an average AUC value exceeding 0.7 were retained and organized into a structured dataset. Finally, by integrating spatial transcriptomics data, cell types were inferred.

Analysis of immune cell infiltration levels and immune subtypes

The CIBERSORT algorithm was employed to assess the infiltration of immune cell populations, including CD4⁺ T cells, CD8⁺ T cells, B cells, and macrophages within the tumor microenvironment of HCC [18]. Furthermore, we investigated the correlation between immune cell infiltration levels and their respective frequencies. Immune subtypes of HCC patients in the TCGA database were obtained using the TCGAbiolinks package in R language. Additionally, online tools were utilized to calculate the disparity in TIDE scores between high-risk and low-risk groups based on TRGs (p < 0.001). Moreover, we explored the association between the TIDE score and the TRGs risk score.

Potential relationship between TRGs and multiple targeted therapeutic agents

We further investigated the role of TRGs in the response to drug therapy for HCC. By analyzing the correlation between the TRGs expression levels and the median inhibitory concentration (IC50) of targeted therapeutic drugs used in HCC patients, we aimed to gain deeper insights into the predictive role of TRGs in HCC drug therapy response. Additionally, we utilized the oncoPredict package in R language for further analysis.

MR analysis of causative TRGs in alcoholic HCC

MR approach is an emerging analytical method that utilizes genetic variants as instrumental variables (IVs) to evaluate the causal effects of modifiable exposures on outcomes [19]. Initially, data were obtained from a genome-wide association study (GWAS), with a significance threshold of p < 5*e-08, to screen for single nucleotide polymorphisms (SNPs) strongly associated with alcoholic HCC. Subsequently, multivariate testing was conducted using both the MR-Egger and inverse variance-weighted IVW methods. A p-value greater than 0.05 indicated horizontal multivariate validity. Heterogeneity was assessed using Cochran’s Q statistic, and a Q-value above 0.05 suggested no heterogeneity among studies included in the analysis. The statistical tests were performed using the Mendelian Randomization package.

Cell culture and low expression vector transfection

The normal hepatocyte cell line LO2 and HCC cell lines 97H, LM3, Hep G2, and 7721 were cultured in DMEM medium supplemented with 10% fetal bovine serum (Eisentec), 100 u/ml penicillin, and 100 ug/ml streptomycin (Solarbio). The cells were incubated at 37 °C in a CO₂ incubator. The medium was changed every other day. The cells were divided into the vector group and the ASF1A knockdown group. ASF1A vector was purchased from Wuhan Sanying Biotechnology Co. For knockdown experiments, a total of 2.5ug knockdown da vector and 5 μl Lipofectamine 3000 (Invitrogen, USA) were diluted with 125 μl Opti-MEM.

MTT experiment

Transfected Hep G2 and 7721 cells were seeded into 96-well plates containing 100 μL of cell suspension (5000 cells per well) and incubated for 3 days. MTT solution (20μL per well) was added, followed by a further incubation period of 4 h. Subsequently, an additional volume of 150 μL DMSO was introduced to each well, and the wells were vigorously shaken for 10 min. Optical density values at a wavelength of 450 nm were measured using an enzyme marker (Thermo Fisher).

Wound scratch test

The Hep G2 and 7721 cells were cultured in 12-well plates and manually induced to form wounds by scratching. The rate of wound healing or cell migration was assessed by measuring the area devoid of cells at 0 and 24 h, respectively.

Colony formation experiment

The HCC cells (1000 per well) were inoculated into 6-well plates. After 7 days of culture, the cells were stained with a solution of 0.1% crystal violet. Cell clones containing more than 50 cells were observed and quantified manually using a light microscope.

Transwell migration experiment

We assessed HCC cell migration using the Transwell method, where 3 × 10*4 cells were placed in the upper chamber with a serum-free medium. Serum medium was added to the lower chamber and incubated at 37 °C for 24 h. Afterward, non-invasive cells on top of the membrane were washed away, while migrating cells were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet. The total cell count was determined under a light microscope.

Western blot experiment

Cells were lysed in cell lysis buffer and their protein concentration was determined using the BCA kit. Subsequently, samples were loaded onto a sodium dodecyl sulfate–polyacrylamide gel and electrophoretically transferred to a PVDF membrane. The PVDF membranes were blocked with 5% skimmed milk for 1 h, followed by overnight incubation at 4 °C with primary antibody. After washing three times with Tris-buffered saline containing Tween (TBST) for 10 min each time, the membranes were incubated with HRP-labeled secondary antibody at room temperature for 1 h. Finally, the membranes were visualized using Clarity Western ECL reagent and radiographic autoradiographic film exposure. Antibodies used included ASF1A (11224-1-AP; Protein tech, Wuhan, China; dilution, 1:400) and Tubulin (11224-1-AP; Protein tech, Wuhan, China; dilution, 1:200).

Statistical analysis

The statistical analyses were performed using the R software package. TRG risk scores were identified as independent prognostic factors for HCC through univariate Cox regression analysis. Differences in OS between different risk subtype groups were compared using Kaplan–Meier analysis, and ROC curves were utilized to evaluate the prognostic value of risk scores with the assistance of the R package. Statistical significance was defined as p-values < 0.05.

Results

Identification of TRGs in HCC

Based on the TCGA database, we performed a comprehensive analysis using 374 HCC samples and 50 paracancerous tissue samples. Through the utilization of PCA, we investigated the distribution of differential genes between tumor and non-tumor groups. Our findings demonstrated distinct distributions of these genes in both groups (Fig. 1A), indicating that PCA is an effective method for distinguishing between tumor and non-tumor groups. Subsequently, we intersected these differential genes with 2093 telomere-associated genes and identified those exhibiting significant differences in expression between HCC and normal liver tissues (FDR < 0.05, |log2 FC|> 1.5). This analysis yielded a total of 935 up-regulated genes and 14 down-regulated genes. A heat map displaying the top 100 differentially expressed genes was generated (Fig. 1B). Furthermore, a one-way Cox analysis was conducted to explore the prognostic significance of these telomere-associated genes, resulting in the identification of 704 survival-related TRGs (Fig. 1C). To enhance their prognostic value, LASSO regression analysis was employed to identify the most informative predictive markers among these differential genes. Based on an optimal penalty parameter value λ, we selected a subset of 18 telomere-associated genes for model construction (Fig. 1D). Finally, utilizing scores derived from these selected TRGs with non-zero coefficients, a risk score model was developed (Fig. 1E). Risk score = Risk score = (-0.0288* ACADS) + (0.1574* MMP1) + (-0.0078* ACAT1) + (0.0344 *ACOT7) + (-0.0623* ACSL6) + ( 0.0744* AGPAT5) + (0.0283* FLAD1.) + (0.0214* PDSS1.) + (0.2112* HSPD1.) + (0.0317 * FKBP1A) + (0.0162 *AKR1B10) + (-0.0879* PDE2A) + (0.0605 *HDAC1) + (0.0482 *HDAC2) + (0.2339 *MAPT) + (-0.0611 *PON1).

Fig. 1

The construction of a risk-prognosis model associated with telomeres in a HCC cohort. A PCA comparing TCGA HCC tumor and non-tumor groups. B Heatmap displaying the top 100 differentially expressed TRGs in HCC. C Volcano plot highlighting prognostically significant differentially expressed TRGs in HCC. D Determination of the optimal parameter (λ) through LASSO regression analysis. E Spectrum of LASSO regression coefficients for 18 telomeric-associated genes

Prognostic risk modeling to assess clinical characteristics of HCC patients

To investigate the independent prognostic potential of the risk score, we initially examined the expression levels of 18 model genes in HCC and adjacent tissues. The results depicted in Fig. 2A illustrated a significant upregulation of these 18 genes in HCC tumor samples compared to adjacent tissues (p < 0.01). Subsequently, based on the median risk score, patients from the TCGA HCC training set were categorized into high-risk and low-risk groups, each consisting of 185 individuals. Kaplan–Meier analysis revealed that within the TCGA cohort, patients classified as low-risk exhibited significantly better OS than those classified as high-risk (p < 0.001) (Fig. 2B). Subsequently, the test group GSE14520 was divided into high-risk (n = 89) and low-risk (n = 139) groups according to their respective RS values obtained from the training set. Kaplan–Meier survival curves generated for this test group demonstrated a significantly longer predicted survival time for individuals assigned to the low-risk group compared to those in the high-risk group (p < 0001) (Fig. 2C). Univariate Cox regression analysis of TCGA training group identified Age (HR = 1.017, 95% CI:1.002–1.032, p = 0.03), stage (HR = 1.389, 95% CI:1.111–1.737, p = 0.004), and risk score (HR = 5 0.097,95%CI:3 0.273 -7 0.936, p < 0 0.001) as three independent prognostic factors for HCC patients (Fig. 2D). The feasibility of the risk score as a prognostic factor was further explored through multivariate Cox regression analysis. The results demonstrated that the TRGs risk score, could serve as an independent prognostic factor for patients with HCC (HR = 4.888, 95% CI: 3.324–7.188, p < 0.001, Fig. 2E). In addition, we re-selected a set of GEO datasets and chose GSE76427 as the test group for univariate and multivariate Cox regression analysis. Univariate Cox regression analysis showed that the risk score was significantly associated with the prognosis of HCC patients (HR: 16.153, 95% CI: 1.442–180.908, p = 0.024, Supplementary Fig. 2A). Multivariate Cox regression analysis revealed that the risk score was an independent prognostic factor for patients with HCC (HR: 20.003, 95% CI: 1.441–277.728, p = 0.026, Supplementary Fig. 2B). Furthermore, the time-dependent ROC at five years was constructed to evaluate the relationship between the telomere-associated risk model and clinical characteristics such as age, gender, tumor grade, and stage. The results indicated that our telomere-associated risk model (AUC = 0.801) outperformed traditional clinical factors including age (AUC = 0.531), gender (AUC = 0.509), tumor grade (AUC = 0.499), and stage (AUC = 067) (Fig. 2F). The time-dependent receiver operating characteristic curve in the test group revealed a one-year AUC of 0.801, a three-year AUC of 0.734, and a five-year AUC of 0.690, indicating that the risk profile consistently performs well in predicting overall survival in HCC patients (Fig. 2G). These findings underscored the superior accuracy of our prognostic risk model. 

Fig. 2

The prognostic validation of 18 TRGs. A Expression levels of these 18 TRGs were examined in both HCC tissues and adjacent tissues. B Kaplan–Meier curve analysis was performed to compare OS among different risk subtypes in the TCGA cohort. C Kaplan–Meier curve analysis demonstrated differences in OS between various risk subtypes within the GSE14520 cohort. D Forest plot displays the results of univariate Cox regression analysis conducted in the TCGA cohort. E Forest plot presents findings from multivariate Cox regression analysis carried out in the TCGA cohort. F ROC curve analysis was employed to evaluate the predictive value of clinicopathological parameters and risk scores in the TCGA cohort. G The area under the time-dependent ROC curve for 1-, 2-, and 3-year OS was calculated for evaluation purposes within the GSE14520 cohort

Constructing column-line graphs to predict TRGs

To enhance the accuracy of predicting survival in HCC patients, we applied the characterization of telomere-related genes into clinical utility. This was achieved by integrating risk and clinical characterization scores to generate column line plots that incorporated independent risk factors such as age, sex, RS, grading, and staging. These plots were used to evaluate and quantify patient survival after 1, 3, and 5 years (Supplementary Fig. 3A). Furthermore, Nomogram modeling demonstrated almost identical survival rates at 1-, 3-, and 5-year intervals between actual observations and predicted outcomes (Supplementary Fig. 3B).

Enrichment score of telomere-related genes in scRNA-seq

The preprocessed single-cell sequencing data was further used to assess the distinctions between hepatocellular carcinoma and normal samples. PCA with a resolution of 0.1 was initially employed for dimensionality reduction (Supplementary Fig. 4), generating 9 cellular subclusters (Fig. 3A). Subsequently, the cell subpopulations were annotated based on the known marker genes (Fig. 3B). The heatmap showed that these subpopulations were in NKT cells, liver cells, hematopoietic stem/progenitor cell-like cells, macrophages, and B lymphocytes (Fig. 3C). Utilizing the UMAP algorithm, we conducted dimensionality reduction analysis on the TRGs of individual cells post-scoring to ascertain their distribution across various cell types. Our findings revealed the elevated TRG scores in three specific cell types: B cells, macrophages, and NKT cells (Fig. 3D). Ultimately, an examination was conducted on the expression of 18 pivotal TRGs across various cell types. As indicated by the Fig. 3E, ASF1A exhibited predominantly elevated expression in NKT cells, while demonstrating diminished expression in Hepatocytes.

Fig. 3

Single-cell RNA sequencing analysis of HCC. A Examination of single-cell RNA sequencing in HCC tissue shows that all cells are grouped into 9 distinct subclusters, as illustrated in the UMAP plot. B Heatmap of the correlation between the marker gene and 5 cell types. C The cells have been labeled as NKT cells, hepatocytes, HPC-like cells, macrophages, and B cells, among several others. D Following the scoring of TRGs, a UMAP visualization map depicting cell distribution was generated. E A visual heat map depicting the expression patterns of 18 crucial TRGs across five distinct cell types

Fig. 4

The investigation into the differentiation pathway of NKT cells. A NKT cells have been clustered into 5 subpopulations in the UMAP plot. B NKT cells undergo differentiation into five unique temporal cell population states. C A scatter plot of different state distribution in two dimensions. D Graph of Pseudo-Time Trajectories. E The temporal changes in the expression levels of RUVBL1, SLC7A11, TCOF1, and ASF1A, along with the two-dimensional scatter plot depicting pseudotime trajectories

Analysis of differentiation trajectories of TRGS and NKT cells

To investigate the relationship between 18 TRGs and NKT cells, we conducted a separate dimensionality reduction analysis on NKT cells and processed them at a resolution of 0.1 (Supplementary Fig. 5). Through pseudotime trajectory analysis, we have identified that NKT cells can be categorized into 5 distinct temporal state cell clusters (Fig. 4A-B). Based on the findings presented in Fig. 4C-D, within these 5 differentiating temporal state cell clusters, 4 were situated in an early differentiation stage (specifically states 1, 2, 3, and 5), while one was positioned in a late differentiation stage (state 4). Following this, an analysis was conducted on the expression of 18 TRGs at various stages of NKT cell differentiation. The findings indicated that RUVBL1, SLC7A11, and TCOF1 genes exhibited high levels of expression during the terminal differentiation stage (Fig. 4E-G), while other genes like ASF1A were predominantly expressed during the early stages of NKT cell differentiation (Fig. 4H).

Analysis of TRGs expression at the ST level

We conducted single-cell sequencing data analysis of 18 TRGs in HCC, and subsequently investigated the expression of these 18 key genes in HCC using ST data. In the ST data of HCC, we observed the comprehensive gene expression profile (Fig. 5A). Following the exclusion of ribosomal and mitochondrial genes, we applied the SCTransform method for data standardization. Subsequently, we further elucidated the dimensionality reduction fractions through PCA (Supplementary Fig. 6) and elbow plot (Supplementary Fig. 7) analysis. This approach facilitated dimensionality reduction clustering, leading to the identification of 9 distinct cell subpopulations (Fig. 5B-C). Following this, FindMarkers was employed to identify 9 genes with spatial variability, and their expression levels on the tissue section are illustrated in Fig. 5D. Additionally, the expression patterns of 18 TRGs in the ST dataset are shown in Fig. 6E-F. The violin plot provides a clear depiction of the spatial transcriptome expression patterns for these 18 genes across different regions. It is evident from the plot that PPM1G, SMG5, and TALDO1 exhibited the highest levels of expression (Fig. 5G-H).

Fig. 5

Assessment of the expression of TRGs in the ST of HCC. A The overall gene expression profile in spatial transcriptomics data of HCC. B After undergoing clustering and dimensionality reduction using UMAP, the spatial transcriptomics data of liver cancer were classified into 9 distinct categories. C Displaying the 9 primary regions visually. D The genetic expression of nine highly mutable genes within ST. E, F The expression and spatial distribution of 18 TRGs in the ST data of HCC. G, H The violin plot illustrates the expression profiles of 18 TRGs

Fig. 6

Analysis of immune infiltration associated with TRGs. A Evaluation of immune cell composition in the TME of HCC across different risk subtypes. B Investigation into the correlation between risk characteristics and infiltrated immune cells. *: p < 0.05 **: p < 0.01 ***: p < 0.001. C Assessment of immune subtypes within distinct risk groups of TRGs. D Comparison of TIDE scores among different risk subtypes. E Examination of the relationship between TIDE scores and TRGs risk scores

Analysis of the correlation between TRGs and the TME in HCC

To further investigate the relationship between 18 TRGS and the TME in HCC, we used SPOTlight to perform inverse convolution on single-cell data to infer the main distribution of five TME cell types at various spatial locations. The expression profile of five types of TME cells was depicted in Supplementary Fig. 8A. The heatmap of correlations showed a strong positive relationship between macrophages and hepatocytes and between NKT cells and macrophages. Conversely, B cells correlated negatively with HPC-like cells (Supplementary Fig. 8B). Supplementary Fig. 8C vividly illustrated the proportion relationship between different cell types, from which we can observe that the proportion of Hepatocytes to B cells, HPC-like cells, macrophages, and NKT cells was relatively limited. Furthermore, the network graph analyzed the interconnections between the five TME cells (Supplementary Fig. 8D). Ultimately, we conducted a visual analysis of the expression of TME cells on the ST (Supplementary Fig. 8E). By comparing the spatial distribution of 18 TRGs and 5 TME markers on ST, it was observed that ASF1A exhibits predominant expression in HPC-like cells.

Relationship between TRGs and immune infiltration of HCC

To gain further insights into the distinct composition of the TME between high-risk and low-risk groups, we employed CIBERSORT to analyze the distribution of 22 immune cell types within HCC TME. The results obtained from the CIBERSORT algorithm were visually presented as a stacked plot (Fig. 6A). Moreover, we investigated the relationship between key genes and risk scores derived from 18 TRGs in HCC, focusing on their association with immune cell infiltration. Notably, Fig. 6B demonstrated significant correlations between a majority of these key genes and CD4⁺ T cells as well as macrophages exhibiting an M0 phenotype. Next, we analyzed the disparities in immune phenotypes between the high-risk group (N = 182) and the low-risk group (N = 178), comprising 360 HCC patients from the TCGA database. The immune subtypes were classified into six main groups as follows: C1, associated with wound healing; C2, characterized by IFN-g dominance; C3, linked to inflammation; C4, indicative of lymphocyte depletion; C5, representing immune quiescence; and C6, marked by TGF-β dominance. In the high-risk group for TRGs, subtype C4 accounted for 43% prevalence while subtype C3 represented 27%. Subtypes C1 and C2 constituted 12% and 18%, respectively. Conversely, in the low-risk TRGs group, there was a higher occurrence of subtypes C3 at 48% and subtype C4 at 45%, with only minimal representation of subtypes C1 at 1% and subtype C2 at 6%. Notably, p = 0.001 indicated a significant disparity in immunophenotypic composition between the low-TRGs risk group and the high-TRGs risk group (Fig. 6C). The correlation between TIDE scores and TRG risk scores was subsequently examined, revealing a significant elevation in TIDE scores within the high-TRG risk group compared to the low-risk group (p < 0.001, Fig. 6D). Furthermore, a robust positive correlation between these two variables was observed (R = 0.47, p < 2.2e -16, Fig. 6E).

Importance of TRGs in chemotherapy and immunotherapy

To gain further insights into the variations in drug resistance levels among different TRG risk subgroups of HCC patients, we conducted a comparison of IC50 levels for nine conventional chemotherapeutic agents between low and high-risk subgroups. Notably, significant differences in IC50 values were observed across risk groups for all nine representative drugs. Specifically, higher IC50 values were identified for 5-Fluorouracil, Dasatinib, Fulvestrant, Gefitinib, and Ipatasertib in the low-risk group, suggesting their potential suitability as treatment options for patients with lower RS in HCC TRG (Supplementary Fig. 9A-E). Conversely, Sorafenib, Axitinib, Cisplatin and Entospletinib exhibited particular sensitivity in the high-risk group which may render them more applicable to HCC patients with higher RS (Supplementary Fig. 9F-I).

MR identifies the involvement of the ASF1A gene in the pathogenesis of alcoholic HCC

Among the 18 TRGs selected by this risk model, we further investigated the TRG genes associated with alcoholic HCC development and identified a potential involvement of ASF1A in the pathogenesis of alcoholic HCC. To validate our hypothesis, MR analysis was conducted to examine the causal relationship between the ASF1A and alcoholic HCC. Five SNPs were obtained from the GWAS website for predicting the susceptibility of the ASF1A to alcoholic HCC occurrence (p < 5e-8, r² > 0.001, kb = 10,000, F > 10) (Fig. 7A, Table 1). The leave-one-out method demonstrated that excluding each SNP did not significantly alter the overall error line, indicating strong reliability and stability of the ASF1A gene in promoting alcoholic HCC occurrence (Fig. 7B). The MR-Egger intercept test showed p = 0.50 > 0.05, suggesting no significant pleiotropy of the ASF1A gene in relation to alcoholic HCC occurrence (Table 2). Subsequently, sensitivity analysis using Cochran’s Q test MR Egger regression (Q = 0.94) and Inverse variance weighted method (Q = 0.92) indicated less biased results caused by SNPs (Fig. 7C). Finally, the IVW method revealed a significant association between the increased expression level of the ASF1A and higher incidence of alcoholic HCC with an odds ratio of OR = 3.325 and a confidence interval for OR: 0.014–0.126; p = 0.02 (Table 3). The scatter plot visually depicted the causal relationship between ASF1A and the occurrence of alcoholic HCC (Fig. 7D). These findings collectively support a causal link between the ASF1A and the development of alcoholic HCC.

Fig. 7

Mendelian randomization analysis model demonstrating a direct association between ASF1A and alcoholic HCC. A Forest plot of the Mendelian randomization analysis assessing the causal relationship between ASF1A and alcoholic HCC. B Residual sensitivity analysis investigating the impact of ASF1A SNP on alcoholic HCC. C Funnel plot of the Mendelian randomization analysis evaluating the causal relationship between ASF1A and alcoholic HCC. D Scatter plot illustrating the Mendelian randomization analysis examining the causal relationship between ASF1A and alcoholic HCC

Table 1 Five SNPs associated with both HCC and ASF1A

Table 2 The Egger regression test evaluates potential multivariate interaction effects

Table 3 Five methodologies were employed to evaluate the causal association between ASF1A and the susceptibility to HCC

Validation of ASF1A in vitro and in clinical samples

To investigate the biological function of the TRG gene ASF1A in HCC progression, a series of in vitro experiments were conducted for validation purposes. Firstly, the expression of ASF1A was examined in the normal hepatocyte cell line LO2 and four common HCC cell lines. The results depicted in Fig. 8A demonstrated a significant elevation of ASF1A expression specifically in HCC cell lines Hep G2 and 7721. Subsequently, siRNA was employed to silence ASF1A expression in these two aforementioned cell lines. The specific sequences for ASF1A-specific siRNAs are provided in Supplementary Table 4. Western blot analysis performed 48 h after transfection revealed that both si-ASF1A-#1 and -#2 effectively inhibited the expression of ASF1A (Fig. 8B). Based on these experimental findings, these two siRNAs were ultimately selected for subsequent studies. Next, we investigated the impact of ASF1A knockdown on HCC cell proliferation. The MTT assay results demonstrated significant inhibition of proliferation in Hep G2 and 7721 cells transfected with si-ASF1A-#1 or -#2 (Fig. 8C). Furthermore, the Transwell assay revealed that si-ASF1A reduced the migratory and invasive capabilities of Hep G2 and 7721 cells (Fig. 8D). A 24-h wound healing assay confirmed that ASF1A inhibition effectively attenuated the migration ability of HCC cells (Fig. 8E). Additionally, the colony formation assay indicated that si-ASF1A suppressed HCC cell growth (Fig. 8F). In conclusion, these findings collectively demonstrate that the downregulation of ASF1A significantly inhibits both proliferation and motility in HCC cells.

Fig. 8

Biological function of ASF1A in HCC. A Expression levels of ASF1A in normal hepatocytes and HCC cell lines. B Western blot analysis demonstrated effective inhibition of ASF1A expression by si-ASF1A-#1 or -#2. C MTT assay revealed that si-ASF1A suppressed the proliferation of Hep G2 and 7721 cells. D Transwell assay showed that si-ASF1A inhibited the migration and invasion abilities of Hep G2 and 7721 cells. E Wound healing assay demonstrated reduced migration and invasion capacities of Hep G2 and 7721 cells upon treatment with si-ASF1A. F In vitro colony formation assay indicated growth inhibition of Hep G2 and 7721 cells by si-ASF1A treatment. Magnification at 100 × , * p < 0.05, ** p < 0.01

To better investigate the potential clinical application of ASF1A in HCC, we further explored the expression of ASF1A in clinical samples. We found that the expression levels of ASF1A mRNA were significantly increased in HCC tissues from the TCGA database. We further analyzed the ASF1A protein levels from the CPTAC database, and the results demonstrated that the ASF1A protein was elevated in HCC. Moreover, by using KM-plotter, we explored that elevated expression of ASF1A was associated with an inferior overall survival (OS), progression free survival (PFS), and relapse free survival (RFS) for the HCC patients (Supplementary Fig. 10). The results indicated that ASF1A maybe potential biomarker for HCC patients.

Discussion

In this study, we identified differentially expressed TRG genes in HCC by analyzing transcriptomic gene expression data from HCC samples and normal liver tissue samples in the TCGA database. Subsequently, univariate analysis and LASSO regression analysis were employed to identify the most significant TRGs for predicting the prognosis of HCC. Ultimately, a prognostic model consisting of 18 selected TRGs was constructed. This study represents the first investigation into the prognostic value of TRGs in HCC and demonstrates that this model can serve as a foundation for selecting therapeutic agents targeting HCC. Immune subtypes C1-6 are a panel of surface markers associated with the immune system and commonly utilized for discerning distinct types of immune cells. Accumulated evidence has demonstrated that different immune subtypes exhibit varying prognoses in diverse tumors [20]. Specifically, the C3 subtype demonstrates the most favorable prognosis, while the C4 and C6 subtypes display the poorest prognosis [21]. Our findings solely investigated immunosubtypes C1-4, revealing that immunosubtype analysis indicated significantly lower levels of C3 in the high-risk group of TRGs compared to the low-risk group. This outcome suggests that patients in the high-TRGs group have a poorer prognosis than those in the low-TRGs group. Furthermore, TIDE was significantly higher in the high-risk group of TRGs than their low-risk counterparts, indicating a potential close association between tumor immune escape and this high-risk subgroup. Subsequently, we investigated the variations in resistance levels among different subtypes and our findings demonstrated that distinct risk subtypes exhibited varying sensitivities to nine conventional chemotherapeutic agents for HCC, thereby offering potential therapeutic targets for patients with diverse subtypes of HCC.

Telomeres, which are specialized structures located at the ends of chromosomes, play a crucial role in safeguarding the stability and integrity of chromosomes [22]. Telomere-associated genes encompass telomerase, telomeric repeat sequences, telomerase-activating receptors, and telomerase inhibitors [9, 23]. Aberrant expression of these genes is closely linked to the occurrence and progression of various neoplastic diseases such as renal cancer [14], lung adenocarcinoma [15], pancreatic cancer [13], and glioma [24].

The research on the impact of TRGs on HCC remains grossly inadequate. Within our risk model, we have successfully identified 18 TRGs that exhibit a profound association with HCC, showcasing their pivotal roles across various tumor types. As an integral ribosomal protein, ASF1A assumes a crucial function in cellular protein synthesis. A plethora of accumulating evidence has underscored the significant involvement of ASF1A in tumorigenesis and disease progression. Notably, the expression level of ASF1A correlates directly with cancer differentiation, invasion potential, and overall prognosis; for instance, in cases of lung cancer, there is a remarkable upregulation of ASF1A expression which closely aligns with an unfavorable prognosis for patients [25]. In triple-negative breast cancer, ASF1A is believed to be associated with resistance to the chemotherapeutic drug doxorubicin [26]. A study conducted by Wu et al. revealed a significant up-regulation of ASF1A expression in HCC and indicated that high levels of ASF1A predicted a poorer prognosis for HCC patients [27], which aligns with our findings and partially validates the accuracy of the TRGs risk model in predicting HCC patient prognosis. Furthermore, ASF1A may play a role in immune evasion from tumors. Li et al. demonstrated that ASF1A could inhibit M1-like macrophage polarization as well as T cell activation and proliferation, thereby diminishing the immune system’s capacity to target tumor cells [28]. Additionally, ASF1A has been shown to promote clone formation and metastasis of tumor cells, consequently augmenting tumor spreading and invasiveness [29, 30].

In our groundbreaking study, we have successfully confirmed through MR imaging for the first time that ASF1A was aberrantly expressed in alcoholic HCC and closely associated with its pathogenesis. Furthermore, our in vitro results also demonstrated a significant upregulation of ASF1A expression in HCC cell lines. Notably, functional experiments on cells reveal that downregulating the expression of ASF1A remarkably inhibited both proliferation and migration of HCC cells. These compelling findings strongly suggested that ASF1A played a pivotal role in the development of HCC, however, further research was warranted to elucidate its precise molecular mechanism.

The CDCA8 gene plays a pivotal role in orchestrating cell development and proliferation throughout the intricate dance of the cell division cycle [31]. It has been unequivocally demonstrated that aberrant expression of CDCA8 is intricately intertwined with the genesis of various tumors, encompassing lung adenocarcinoma [32], thyroid cancer [33], prostate cancer [34], and HCC [35], thus endowing it with immense potential as a prognostic biomarker. Moreover, CDCA8 has been unambiguously shown to augment migratory prowess and invasiveness across diverse tumor cells, spanning pancreatic cancer cells [36], osteosarcoma cells [37], and bladder cancer cells [38]. In lung adenocarcinoma, HMMR has been identified as a pivotal oncogenic driver closely linked to patient prognosis in this particular malignancy [39, 40]. Notably, Guo et al. have convincingly demonstrated that HMMR exerts additional potent effects on prostate cancer progression and metastasis through modulation of the AURKA/mTORC2/E2F1 signaling cascade [41]. Moreover, IPO13 overexpression significantly contributes to both enhanced proliferative capacity and increased metastatic potential in non-small cell lung carcinoma [42] as well as endometrial carcinoma cases [43]. Furthermore, aberrant MT3 expression acts as an influential oncogenic factor promoting uncontrolled growth and invasive behavior in bladder neoplasms [44]. Prognostic risk models are not only applied in liver cancer, but have also been extensively developed and utilized in various other tumors. For example, cancer-testis antigens can more accurately predict the survival rate of STAD patients, and ELOVL4 has been identified as a potential therapeutic target for gastric cancer [45]. Additionally, cancer-associated fibroblasts can effectively predict the clinical prognosis of ovarian cancer patients, as well as the tumor immune microenvironment and the response to immune checkpoint inhibitors [46]. Meanwhile, another study showed that neurociliin 1 (NRP1), a key gene in tumor-associated fibroblasts, is upregulated in cervical cancer tissue, influencing tumor progression and differentiation processes [47]. Bao et al.’s study suggests that the risk prognostic model of PANoptosis associated long non-coding ribonucleic acidscan be used to predict clinical outcomes in patients with lung adenocarcinoma and provide a theoretical basis for personalized treatment of patients, such as immunotherapy [48]. These signature risk models provide better precision and personalized medicines, thereby enhancing patients’ overall treatment outcomes and quality of life.

However, our study still possesses certain limitations. Firstly, the HCC telomere-associated risk model we constructed using the TCGA database lacks translation and follow-up time for validation purposes. Secondly, although the risk model can predict the prognosis of HCC patients to a certain extent, it is based on transcriptomic features which adds intricacy to its operation. Thirdly, further fundamental in vitro and in vivo experiments are required to refine the specific roles of the 18 TRGs in HCC within the constructed risk model. Lastly, additional clinical experiments are needed to refine the utilization of 9 chemotherapeutic agents that exhibit sensitivity toward different risk subtypes.

Conclusion

The present study has successfully developed a prognostic model consisting of 18 TRGs, which accurately predicted the prognosis and immunotherapy effectiveness for HCC patients. Moreover, it established a causal relationship between ASF1A and alcoholic HCC. This groundbreaking discovery not only enhanced our understanding of HCC development but also provided new insights and targets for preventing and treating HCC patients.