Gene patterns linked to copper-induced cell death and the tumor immune environment in liver cancer

CRGs and clinical information of patients with HCC were obtained from the TCGA database (https://portal.gdc.cancer↗. Gov/) (Tomczak et al., 2015). A total of 371 HCC tissues and 50 paracancerous tissues were collected and used in this study. The data used in this study were standardized data per million transcripts, and the data distribution was close to normal distribution. R software (v4.0.3) was used to extract gene expression data, and a data matrix was constructed for further analysis.

Nineteen CRGs were first extracted from the TCGA expression matrix. According to the consistent clustering of these 19 genes, the R software package Consensus-ClusterPlus (v1.54.0) was used for the analysis of consistency, and the maximum number of clusters was 6 (Wilkerson and Hayes, 2010). The TCGA-LIHCREAD cases were then divided into 2 clusters based on the expression profiles of CRGs. This process was repeated 500 times to ensure the stability and reproducibility of the classification.

DEGs between C1 and C2 subgroups were obtained using Limma package (version 3.40.2) in R software (Ritchie et al., 2015). Adjusted p values were analyzed in the TCGA database to correct for false positive results. “Adjusted p < 0.05 and log2 (fold change) > 1 or log2 (fold change) < −1″ was defined as the criteria for screening the differential expression of mRNA. GeneMANIA (http://www.genemania.org↗), which is a software that elucidates the relationship between genes and datasets by building gene interaction networks (Warde-Farley et al., 2010), was used to visualize the gene network of CRGs through physical interaction, co-expression, prediction, co-localization, and genetic interaction, and to evaluate their functions.

The GO function and the enrichment of KEGG pathways were analyzed by the “clusterProfiler” package in R software (Yu et al., 2012). Potential biological pathways were also identified by GSEA (http://software.broadinstitute.org/gsea/index.jsp↗) (Powers et al., 2018). According to the data of TCGA, the DEGs were divided into upregulated and downregulated groups. Each analysis performed 10,000 gene combination permutations to identify signaling pathways with significant changes; the genes were considered as enriched for meaningful signaling pathways when p.adjust <0.05 and FDR (false discovery rate) < 0.25. Statistical analysis and graphing were performed using the R package clusterProfiler (3.18.0).

The R software package immunedeconv was used for immune score assessment to compare the degree of immune cell infiltration in C1 (N = 252) and C2 (N = 119) subgroups by Wilcoxon test, by integrating six state-of-the-art algorithms, including TIMER (Li et al., 2016), xCell (Aran et al., 2017), MCP-counter (Becht et al., 2016), CIBERSORT (Newman et al., 2015), EPIC (Racle et al., 2017) and quantTIseq (Finotello et al., 2019; Sturm et al., 2020). The expression of some immune checkpoint-related genes was also analyzed. The results were visualized by the R (v4.0.3) packages “ggplot2″ and “pheatmap".

The expression of CRGs in 371 HCC tissues and 50 paired adjacent tissues was analyzed using the TCGA database. In addition, univariate Cox regression analysis was used to investigate the effect of CRGs on the prognosis of HCC. Kaplan-Meier curves, p values, and hazard ratios (HR) with 95% confidence intervals (CIs) were obtained by logrank test and univariate Cox regression. Five CRGs (CDKN2A, DLAT, DLST, GLS, and PDHA1) with higher hazard ratios were screened from the Cox regression analysis plot. The relationship between the prognostic CRGs and the OS rate in HCC patients was also analyzed, and the AUC (Area Under Curve) under the ROC curve was calculated.

The following HCC cell lines and hepatocyte were used in this study: Hep3B, Huh7, HepG2, and L02 (Chinese Academy of Medical Sciences, Beijing, China). All cell lines were maintained in Dulbecco’s modified Eagle’s medium (DMEM; Gibco, Grand Island, NY, United States) supplemented with 10% fetal bovine serum (Gibco, Grand Island, NY, United States of), 100 U/ml penicillin, and 100 U/ml streptomycin (Invitrogen, Carlsbad, CA, United States) and incubated at 37°C in a humidified atmosphere with 5% CO₂.

First, TRIzol reagent (Invitrogen, CA, United States) was used to extract total RNA from the samples. Then, the total RNA concentrations were measured using a NanoDrop 2000c (Thermo Fisher Scientific, Inc.). Next, cDNA was synthesized by reverse transcription using a miScript II RT kit (Qiagen, Germany). Finally, a miScript SYBR Green PCR kit (Qiagen, Germany) was used to detect the expression of target genes on a Lightcycler 480 Real-Time PCR System (Roche Diagnostics GmbH, Mannheim, Germany). The 20-μl reaction mixture included 10 μl 2X QuantiTect SYBR Green PCR Master Mix, 0.8 μl 10X miScript universal Primer, 2 μl cDNA template, 0.8 μl specific miRNA primer, and 7.4 μl RNase-free water. Additionally, the thermocycling steps were as follows: 95°C for 15 min, followed by 40 cycles of 94°C for 15 s, 55°C for 30 s, and 72°C for 30 s. The relative standard curve method (2^−△△CT) was employed to determine the relative mRNA expression, with the glyceraldehyde 3-phosphate dehydrogenase (GAPDH) gene as the reference. Supplementary Table S1 lists the polymerase chain reaction primers used in this study.

A prognostic model was constructed using LASSO-Cox regression analysis based on the above CRGs associated with the prognosis of HCC patients. According to the results of the multivariate Cox regression analysis, the prognostic CRG risk score was calculated as follows: Risk score = ∑i Cofficient(mRNAi) × Expression (mRNAi). Then, TCGA-LIHC patients were divided into low-risk and high-risk subgroups according to the average risk score, Kaplan-Meier analysis was used to compare the OS rates of the two subgroups, and timeROC analysis was performed to predict the accuracy of the model. Each variable (including the p-value, and HR with 95% CI) was displayed by univariate and multivariate cox regression analysis and using forest plots by the “forestplot” package. The “rms” package was used to build a Nomogram model for predicting 1, 3, and 5-year OS and PFS based on the results of multivariate cox proportional hazards analysis.

GEPIA is a newly developed interactive website for analyzing RNA sequencing expression data from 9736 tumors and 8587 normal samples from TCGA and Genotype-Tissue Expression databases, such as differential expression analysis of tumor/normal tissues, tumor type or pathological stage, patients’ survival analysis, similar gene detection analysis, and dimensionality reduction analysis (Tang et al., 2017). The expression data of prognostic CRGs in different stages of HCC were obtained using the “Stage Plot” module in the GEPIA2 database (http://gepia2.cancer-pku.cn/#index↗). The UALCAN database (http://ualcan.path.uab.edu/↗) is an interactive portal containing TCGA RNA transcriptome data and clinical data from 31 cancer types. TISIDB (http://cis.hku.hk/TISIDB↗) is a portal for assessing tumors and the immune system that integrates multiple heterogeneous data types including genomics, transcriptomics, and clinical data from TCGA of 30 cancer types (Ru et al., 2019). The UALCAN and TISIDB databases were used to confirm the association between prognostic CRGs and the clinical stage of HCC patients. p < 0.05 was considered statistically significant.

Immunohistochemical images of CRGs were collected from the Human Protein Atlas (https://www.proteinatlas.org↗) database to assess the differences in CRG expression between HCC and adjacent tissues (Sun et al., 2021).

cBioPortal (http://www.cbioportal.org/↗) provides visualization tools for analyzing cancer genetic data. The genomic maps of CRGs were analyzed using cBioPortal based on the TCGA database to understand the mutation frequency in HCC (Gao et al., 2013). In this study, 372 LIHC samples was selected to explore the impact of prognostic CRGs on survival of HCC patients. In addition, the effects of CRG mutations on cancer signaling pathway expression, TMB, MSIsensor scores, and hypoxia-related scores (Ragnum Hypoxia Score, Buffa Hypoxia Score, and Winter Hypoxia Score) were examined. p < 0.05 was considered statistically significant.

The effects of prognostic CRGs on the abundance of infiltrating immune cells in tumors were analyzed by the TIMER (https://cistrome.shinyapps.io/timer/↗) database (Li et al., 2017). The detected immune cells included tumor purity, B cells, CD4⁺ T cells, CD8⁺ T cells, neutrophils, macrophages, and dendritic cells. In addition, the infiltration of immune cell types was quantified by the R package “GSVA” single sample gene set enrichment analysis (ssGSEA) (Hanzelmann et al., 2013). Next, the Spearman method was used to explore the correlation of prognostic CRGs with the level of immune cell infiltration. Finally, the relationship between prognostic CRGs and checkpoints (CD274, CTLA4, and PDCD1) was examined.

The Spearman’s method was used to analyze the correlation between prognostic CRGs and MSI in HCC. The chemotherapeutic response was predicted for each sample using the GDSC database (https://www.cancerrxgene.org/↗) (Yang et al., 2013). The R package pRRophetic was used to predict the half maximum inhibitory concentration (50% inhibition of the concentration, IC50) of the chemotherapeutic drugs, in which the IC50 of the sample was estimated by ridge regression. All parameters were set by default, using the batch effect of combat and the organization types of all, and summarizing the replicate gene expression as an average. Drug sensitivity and gene expression profiling data from cancer cell lines in the GDSC database were integrated in this study.

The correlation between 19 CRGs and m6A-related genes expression and prognostic in 371 HCC samples was analyzed by R (4.0.3) packages “ggplot2.” The difference in m6A-related gene expression between C1 group (N = 252) and C2 group (N = 119) was investigated. The m6A-related genes analyzed included METTL3, YTHDC1, YTHDC2, METTL14, RBM15, RBM15B, IGF2BP1, IGF2BP2, IGF2BP3, VIRMA, WTAP, YTHDF1, YTHDF2, YTHDF3, ZC3H13, HNRNPA2B1, HNRNPC, RBMX, FTO, and ALKBH5.

The effect of prognostic CRGs on the expression of single cell subsets in the TME was investigated using TISCH (http://tisch.comp-genomics.org/↗) (Sun et al., 2021). TISCH is a scRNA-seq database focused on the TME, providing detailed cell-type annotation at the single-cell level, which is beneficial for exploring the TME in different cancer types. In this dataset, three main cell types are present, such as immune cells, stromal cells, and malignant cells. In this study, the t-distributed stochastic neighborhood embedding (t-SNE) map of LIHC_GSE125449_aPDL1aCTLA4 and the heatmap of LIHC_GSE125449_aPDL1aCTLA4 were displayed through the TISCH database to demonstrate the effect of CRGs on the TME of HCC.

Potential miRNA targets of prognostic CRGs were predicted using the ENCORI (http://starbase.sysu.edu.cn/↗) database (Li et al., 2014), RNAInter (http://www.rnainter.org/↗) (Lin et al., 2020) and RNA22 (https://cm.jefferson.edu/rna22/interactive↗) database (Loher and Rigoutsos, 2012), and the prognostic value of these potential miRNA targets in HCC was confirmed using ENCORI, Kaplan-Meier Plotter and TCGA-LIHC cohort. Then, the lncRNAs potentially binding to prognostic miRNAs were predicted by the miRNet (http://www.mirnet.ca/↗) database (Chang et al., 2020). Subsequently, a miRNA-lncRNA regulatory network was established by Cytoscape (version 3.7.1; http://www.cytoscape.org/↗) software (Shannon et al., 2003). In addition, the prognostic value of these potential lncRNA targets in HCC were further verified. Finally, a lncRNA-miRNA-mRNA regulatory network was established.

The flowchart of the study is illustrated in Figure 1. A total of 371 HCC carcinoma samples were clustered in the TCGA database using consensus clustering T = to identify potential CRG clusters. All tumor samples were divided into k (k = 2–6) distinct clusters based on the expression of 19 CRGs in HCC (Figures 2A,B). Then, the number of clusters was selected as 2 according to the cluster analysis results, indicating that the patients with HCC were accurately divided into two clusters (Cluster 1 had 252 patients and Cluster 2 had 119 patients) (Figure 2C). Figure 2D shows the heat map of cuproptosis-related gene expression in different subgroups.

The study included 421 liver HCC samples consisting of 371 tumor samples and 50 adjacent normal samples. The DEGs identified between the C1 and C2 subtypes contained 5044 upregulated genes and 52 downregulated genes. Then, a volcano map (Figure 3A) and a heat map (Figure 3B) were constructed based on these DEGs. The results of GeneMANIA revealed that the functions of the significant genes co-expressed in this network (DLD, GLS, NFE2L2, DLST, LIAS, ATP7B, PDHA1, SLC31A1, FDX1, DBT, GCSH, LIPT1, CDKN2A, DLAT, NLRP3, ATP7A, PDHB, LIPT2, and MTF1) were related to several processes, including oxidoreductase complex, tricarboxylic acid cycle enzyme complex, cellular amino acid catabolic process, dihydrolipoyl dehydrogenase complex, acyl-CoA metabolic process, acetyl-CoA biosynthetic process and thioester metabolic process (Figure 3C). The identified up and downregulated DEGs were then further subjected to GO and KEGG enrichment analysis. The results of biological process (BP) analysis showed that the DEGs were mainly enriched in extracellular matrix tissues, regulation of DNA metabolism, response to transforming growth factor β, activation of protein kinase activity, and cell cycle G1/S phase transition. The results of cellular component (CC) analysis showed that the DEGs were mainly enriched in focal adhesions, cell-matrix adhesion junctions, extracellular matrix components, collagen-containing extracellular matrix, and fibrous collagen trimers. The results of molecular function (MF) analysis showed that the DEGs were mainly enriched in several functions including cell adhesion molecule binding, extracellular matrix structural components, growth factor binding, protein serine/threonine kinase activity, and collagen binding. The KEGG analysis showed that the DEGs were mainly enriched in focal adhesion, extracellular matrix-receptor interaction, PI3K-Akt signaling pathway, TGF-β signaling pathway, HCC, cell cycle, NOD-like receptor signaling pathway, PD-L1 expression and the PD-1 checkpoint pathway in cancer (Figure 3D). The results of Gene Set Enrichment Analysis (GESA) showed that CRGs were closely related to cancer pathways in HCC, including focal adhesions, cell cycle, T cell receptor signaling pathway, JAK STAT signaling pathway and MAPK signaling pathway. The activation of these pathways increased the risk of tumorigenesis and tumor progression (Figure 3E).

The TCGA-LIHC cohort was downloaded to explore the difference in immune cell infiltration between the C1 and C2 subtypes to explore the role of CRGs in the immune response in HCC. Six state-of-the-art algorithms were then integrated for a reliable immune score assessment. quantTIseq, EPIC, MCP-counter, CIBERSORT, xCell and TIMER showed the expression of CRGs was correlated with uncharacterized cell, Neutrophil, NK cell, Macrophage M1, T cell regulatory (Tregs), Macrophage M2, B cell (Figures 4A,B); Macrophage, uncharacterized cell, T cell CD4⁺, NK cell (Supplementary Figure S1A); NK cell, B cell, Myeloid dendritic cell, Monocyte, Macrophage/Monocyte, Neutrophil, Endothelial cell (Supplementary Figure S1B); Macrophage M2, T cell CD8⁺, T cell follicular helper, Mast cell activated, Neutrophil (Supplementary Figure S1C); Macrophage M2, Endothelial cell, Mast cell, T cell CD8⁺ naive, T cell CD4⁺ memory, Plasmacytoid dendritic cell, T cell CD4⁺ naive, Common lymphoid progenitor, T cell NK, T cell CD4⁺ Th1 (Supplementary Figure S1D); T cell CD4, Neutrophil, Macrophage, T cell CD8⁺, B cell, Myeloid dendritic cell (Supplementary Figure S1E) respectively. Finally, the difference between the two subgroups in the expression of eight immune checkpoint-related genes was explored, and the results showed that CD274 (p < 0.01), HAVCR (p < 0.01), PDCD1LG2 (p < 0.01), TIGIT (p < 0.05) and SIGLEC15 (p < 0.05) were highly expressed in the C1 subgroup than in the C2 subgroup (Figure 4C).

The expression of 19 CRGs in HCC and normal tissues was investigated using the TCGA dataset; among them, 15 CRGs were up- or downregulated. The expression of ATP7A, SLC31A1, LIAS, LIPT1, LIPT2, DLD, DLAT, PDHA1, PDHB, MTF1, GLS, CDKN2A, and DLST was upregulated, while NLRP3 and DBT were downregulated in cancer tissues than in the normal tissues (Figure 5A). Most of the 19 CRGs in the HCC samples were positively correlated, while 6 pairs were negatively correlated, showing a strong correlation between them (Figure 5B). Univariate Cox regression analysis was then performed to screen CRGs with a prognostic value (Figure 5C). The mRNA expression of five prognostic CRGs was significantly upregulated (p < 0.001) in HCC cell lines (Hep3B, Huh7 and HepG2) compared to their expression in the hepatocyte cell line (L02) (Figure 5D). Finally, five genes with a prognostic value were identified from the results. The Kaplan-Meier survival curves are shown in Figure 5E. HCC patients had a lower survival rate when CDKN2A (HR = 1.75, p = 0.002), DLAT (HR = 1.7208, p = 0.0024), DLST (HR = 1.5484, p = 0.0141), GLS (HR = 1.5001, p = 0.0222) and PDHA1 (HR = 1.4981, p = 0.0226) were highly expressed.

A prognostic gene model based on prognostic CRGs was constructed by LASSO Cox regression analysis (Figures 6A,B). The risk score was calculated as follows: (0.1242)*CDKN2A+ (0.2058)*DLAT + (0.0795)*GLS. HCC patients were divided into two groups according to the risk score distribution, survival status, and expression of GLS, DLAT, and CDKN2A (Figures 6C,D). Patients had an increased risk of death and decreased survival time as the risk scores increased (Figure 6C). The Kaplan-Meier curves showed that HCC patients with high-risk scores had a lower probability of OS than patients with low-risk scores (median time = 3.4 years vs. 5.8 years, p = 0.00106) (Figure 6D). The area under the ROC curve (AUC) in the 1-year, 3-year and 5-year ROC curves was 0.737, 0.646 and 0.634, respectively (Figure 6E). The results of our constructed cuproptosis-related risk profile showed a significant association with survival in patients with HCC.

Considering the clinicopathologic features and prognostic CRGs, a predictive nomogram was also built to predict the survival probability. The results of the univariate and multivariate analysis showed that CDKN2A expression and pT stage were independent factors affecting the prognosis of HCC patients (Figures 7A,B; Supplementary Figure S2). The predictive nomogram suggested that the 3-year and 5-year OS rates and PFS rates were predicted enough accurately compared with an ideal model in the entire cohort (Figures 7C,D; Supplementary Figure S2).

Based on the above analysis, the correlation between the expression of CRGs and the pathological stage of HCC was analyzed by Gene Expression Profiling Interactive Analysis (GEPIA) 2 database. The results showed that the expression of CDKN2A, DLAT, GLS and PDHA1 were significantly correlated with the pathological stage of HCC (p < 0.05), while the correlation between DLST and the pathological stage was not significant (Supplementary Figure S3). The analysis of UALCAN database revealed that the expression of prognostic CRGs in HCC patients in stages 1–3 was significantly upregulated compared with their expression in the normal group. The analysis of TISIDB database revealed that the expression of CDKN2A, DLST, and GLS in different stages and scores of HCC patients was significantly different, while DLAT and PDHA1 was slightly changed although not significantly (Supplementary Figure S3).

The protein expression of CRGs in 371 HCC and adjacent tissues was measured using the Human Protein Atlas to determine the difference in protein expression of the prognostic CRGs in HCC tissues and normal liver. The results of the immunohistochemical staining showed that the prognostic CRGs were moderately or highly expressed in HCC tissues but were low in paracancerous tissues (Figures 8A,B).

The genetic alteration of the prognostic CRGs were analyzed using the cBio Cancer Genomics Portal (cBioPortal) online tool for HCC. CDKN2A, DLAT, DLST, GLS, and PDHA1 were altered in 80%, 25%, 48%, 37%, and 32% of the 347 samples of the 372 sequenced patients, as shown in Figure 9A. The mRNA expression z-scores of the prognostic CRGs relative to normal samples are shown in Figure 9B. The survival analysis results showed that the genetic alteration in prognostic CRGs was significantly associated with shorter disease free survival (DFS) (p = 0.0377, Figure 9D) and PFS (Figure 9E, p = 0.018) of HCC patients. However, these two groups have little correlation with OS (p = 0.0723, Figure 9C). In addition, the clinical analysis showed that HCC tissues with prognostic CRG mutations upregulated the expression of cancer-related signaling pathways (p = 1.051 × 10⁻³, Figure 9F), tumor burden mutation (TMB) (p = 1.490 × 10⁻³, Figure 9G), MSIsensor Score (p = 1.426 × 10⁻³, Figure 9H), as well as hypoxia scores, including Buffa Hypoxia Score (p = 2.329 × 10⁻³, Figure 9I), Ragnum Hypoxia Score (p = 3.249 × 10⁻⁵, Figure 9J) and Winter Hypoxia Score (p = 9.727 × 10⁻⁴, Figure 9K). Therefore, the mutation of prognostic CRGs in HCC tissues significantly affected the prognosis of HCC patients.

The above results demonstrated that CRGs play a crucial role in the development of tumor immune microenvironment in HCC. The present study further investigated the correlation between the expression of CRGs (CDKN2A, DLAT, DLST, GLS, and PDHA1) and immune infiltration in HCC using the TIMER database and TCGA database. TIMER data showed that CDKN2A, DLAT, GLS and PDHA1 were strongly and positively correlated with tumor purity, B cells, CD4⁺ T cells, CD8⁺ T cells, neutrophils, macrophages, and dendritic cells, while DLST was only positively correlated with tumor purity, CD4⁺ T cells, neutrophils, and dendritic cells (Figures 10A–E). A total of 24 immune cell types infiltrating HCC were next identified using single-sample GSEA, followed by Spearman analysis to investigate the association between prognostic CRGs and immune cell infiltration. The high expression of prognostic CRGs was positively correlated with several cell types, including T helper cells, Th1 cells, TFH, and aDC, while negatively correlated with cytotoxic cells, DC, and pDC, as shown in Supplementary Figure S4. In addition, the correlation between prognostic CRG expression and important immune checkpoints (CD274, PD-L1 and CTLA4) was further investigated. The results showed that CDKN2A was significantly and positively correlated with CD274 (p = 2.41 × 10⁻⁰⁵, cor = 0.217), PD-L1 (p = 1.24 × 10⁻⁰³, cor = 0.167) and CTLA-4 (p = 3.62 × 10⁻⁰⁵, cor = 0.213), and GLS was also positively correlated with CD274 (p = 3.02 × 10⁻⁰⁹, cor = 0.302), PD-L1 (p = 4.25e−07, cor = 0.259) and CTLA-4 (p = 8 × 10⁻⁰⁵, cor = 0.203), while DLAT (p = 3.27e−13, cor = 0.366), DLST (p = 3.71e−10, cor = 0.318) and PDHA1 (p = 1.77e−07, cor = 0.267) were significantly and positively correlated only with CD274 (Supplementary Figure S5). These results indicated a significant correlation between prognostic CRGs and tumor immune infiltration and showed that prognostic CRGs might be predictive markers of anti-CD274 therapy.

Microsatellite instability (MSI) and TMB can be used as indicators to predict the response of various tumors to immunotherapy (Rizzo et al., 2021). The above results showed that the expression of CRGs was closely related to tumor immune cell infiltration in HCC. Next, the association of prognostic CRGs with TMB and MSI in HCC was analyzed to determine whether CRGs can be considered as biomarkers for immunotherapy. The results showed that PDHA1 was positively correlated with MSI (p = 0.044) (Figure 11A). Furthermore, CDKN2A (p = 1.64 × 10⁻⁰⁴) was positively correlated with TMB, and GLS (p = 7.74 × 10⁻⁰⁵) and negatively correlated with TMB (Figure 11B). Finally, the data on the gene expression profiles of cancer cell lines were integrated in the Genomics of Drug Sensitivity in Cancer (GDSC) database for drug sensitivity to fully explore the potential value of CDKN2A, DLAT, DLST, GLS, and PDHA1 genes as novel therapeutic targets. A Pearson correlation analysis of the data was performed, which showed that the expression of CDKN2A, DLAT, DLST, GLS, and PDHA1 was positively correlated with Saracatinib, Selumetinib, Serdemetan, Rucaparib, Roscovitine, Refametinib, Proteasome, Palbociclib, Motesanib, Lapatinib, Imatinib and Erlotinib, but negatively correlated with Tretinoin, Tipifarnib, Nilotinib, Navitoclax, Doramapimod, CCT018159 and AICA_ribonucleotide (Figure 11C).

It is well known that m6A methylation modification affects gene expression by regulating RNA metabolism and plays an important role in the occurrence and development of malignant tumors. The potential correlation of CRGs with m6A modification in HCC was evaluated, and the differential expression of 20 m6A-related genes (RBM15B, VIRMA, IGF2BP29, HNRNPA2B1, IGF2BP1, YTHDF3, IGF2BP3, HNRNPC4, RBM15, RBMX, METTL14, YTHDC2, METTL3, ZC3H13, WTAP, YTHDF1, YTHDC1, FTO, YTHDF2) were found between the C1 and C2 subgroups. The results showed that except IGF2BP1, other m6A-related genes were significantly different between the two groups (p < 0.01, Figure 12A). The C1 group with higher expression of CRGs has also higher expression of m6A methylation-related genes. In addition, the correlation between prognostic CRGs and m6A-related genes was analyzed by the TCGA dataset, and the results showed a significant positive correlation between prognostic CRGs and m6A-related genes (Figure 12B). The above results suggested that CRGs were closely related to m6A modification in HCC.

The TME is composed of extracellular matrix, cancer associated fibroblasts (CAFs), myofibroblasts, immune cells and other factors. The prognostic CRGs were assessed at a single-cell level expression by Tumor Immunity Single Cell Center (TISCH, http://tisch.comp-genomics.org/↗) to examine the relationship between prognostic CRGs and the infiltration of immune cells, stromal cells, and malignant cells. The HCC single-cell GSE dataset was explored (LIHC_GSE125449_aPDL1aCTLA4). T cells, B cells, plasma cells, macrophages, endothelial cells, fibroblasts, malignant cells, and hepatic progenitors were annotated by single-cell RNA sequencing analysis (Figure 13A). The results showed that CDKN2A, DLAT, DLST, GLS, and PDHA1 were correlated in single cell subsets such as immune cells, stromal cells, and malignant cells, and were significantly expressed in macrophages, fibroblasts, endothelial cells, and hepatic progenitor cells (Figures 13B,C). CSFs and tumor-associated macrophages (TAMs) play an important role in the occurrence and development of cancer (Rizzo et al., 2021). Therefore, the association between prognostic CRGs and biomarkers associated with CAFs and TAMs was further explored. The results showed that prognostic CRGs were significantly and positively correlated with CAF markers such as FAP, PDPN, S100A8, S100A9, TGFB1, TGFB2, ACTA2, PALLD, TNC, and COL11A1 (Figure 13D). Prognostic CRGs were also highly correlated with TAM-related markers such as CCL2, CD68, and IL 10 (Figure 13E).

Next, the upstream regulatory miRNAs of prognostic CRGs were predicted using a comprehensive miRNA-related database. Firstly, 27 pairs of CDKN2A-miRNAs, 63 pairs of DLAT-miRNAs and 172 pairs of GLS-miRNAs were obtained by the intersection of ENCORI and RNAInter databases. A total of 110 pairs of DLST-miRNAs and 39 pairs of PDHA1-miRNAs were obtained by the intersection of ENCORI databases and RNA22 databases (Figure 14A). Then, according to the classical mechanism of miRNAs negatively regulating mRNA expression, a negative correlation between the mRNA and the predicted miRNA was expected. The ENCORI database was used to screen the correlation, prognosis and expression of these candidate miRNAs in HCC. Among these miRNA-mRNA interactions, the results showed that 1 pair of miRNA-CDKN2A, 6 pairs of miRNAs-DLAT, 3 pairs of miRNAs-DLST, 18 pairs of miRNAs-GLS, and 2 pairs of miRNAs-PDHA1 were significantly and negatively correlated (Figure 14B; Supplementary Figure S6). Theoretically, miRNAs that bind to prognostic CRGs should be downregulated and indicate poor prognosis in HCC. Therefore, the prognostic role and expression of these potential miRNAs in HCC were further confirmed using the ENCORI database. The results showed that only the low expressed hsa-miR-125b-5p, hsa-miR-101-3p and hsa-miR-23c corresponded to a poor prognosis, and hsa-miR-125b-5p, hsa-miR-101-3p and hsa-miR-23c were significantly lower in HCC tissues than in normal tissues (Figure 14C). All these results indicated that CDKN2A-hsa-miR-125b-5p, GLS-hsa-miR-125b-5p, GLS-hsa-miR-101–3, GLS-hsa-miR-23c and PDHA1-hsa-miR-125b-5p might represent key pathways that mediated the occurrence and development of HCC and were related to the prognosis of patients.

Previous studies showed that lncRNAs bind to miRNAs and mediate the expression of target genes, thereby exerting biological roles (Militello et al., 2017; Zhang and Lou, 2020). Based on the above results, lncRNAs upstream of miRNAs were also predicted to construct the miRNA-lncRNA axis. The lncRNAs potentially binding to hsa-miR-101-3p, hsa-miR-125b-5p and hsa-miR-23c were predicted by the intersection of ENCORI and miRNet databases, and the results revealed 27 lncRNAs targeting hsa-miR-101-3p, 47 lncRNAs targeting hsa-miR-125b-5p and 56 lncRNAs targeting hsa-miR-23c (Figure 15A). A miRNA-lncRNA regulatory network was established by Cytoscape software for a better visualization (Figure 15B). According to the ceRNA hypothesis, lncRNAs can increase mRNA expression by competitively binding to miRNAs. Therefore, lncRNAs were negatively correlated with miRNAs or positively correlated with mRNAs. The correlation between lncRNAs and hsa-miR-101-3p, hsa-miR-125b-5p and hsa-miR-23c was detected by the ENCORI database, and the results showed that 14 lncRNAs were significantly associated with hsa-miR-101-3p and GLS, 9 lncRNAs were significantly associated with hsa-miR-125b-5p and CDKN2A, 10 lncRNAs were significantly associated with hsa-miR-125b-5p and GLS, 4 lncRNAs were significantly associated with hsa-miR-125b-5p and PDHA1, and 12 lncRNAs were significantly associated with hsa-miR-23c and GLS (Supplementary Tables S2–6; Supplementary Figure S7). Subsequently, the prognostic value of lncRNAs in HCC was assessed by the Kaplan-Meier plotter, and the expression of these lncRNAs in HCC were detected using the TCGA-LIHC cohort. The results in combination with the survival analysis and expression analysis, showed that GSEC, PTPRG-AS1, CYTOR, DANCR, TRAF3IP2-AS1 and DLEU2 were significantly upregulated in HCC, and their upregulation was associated with poor prognosis in HCC patients (Figure 15C). Finally, a key mRNA-miRNA-lncRNA triple regulatory network associated with HCC prognosis was established, which included 3 mRNAs (CDKN2A, GLS, and PDHA1), 3 miRNAs (hsa-miR-125b-5p, hsa-miR-101-3p and hsa-miR-23c) and 6 lncRNAs (GSEC, PTPRG-AS1, CYTOR, DANCR, TRAF3IP2-AS1 and DLEU2). Eight lncRNA-miRNA-mRNA networks were also constructed, including lncRNA CYTOR/miR-125-5p/CDKN2A, lncRNA CYTOR/miR-125-5p/GLS, lncRNA DANC R/miR-125-5p/GLS, lncRNA DANCR/miR-125-5p/PDHA1, lncRNA DLEU2/miR-23c/GLS, lncRNA GSEC/miR-101-3p/GLS, lncRNA PTPRG-AS1/miR-101-3p/GLS and lncRNA TRAF3IP2-AS1/miR-23c/GLS regulatory axes. (Figure 15D).

The datasets are available in TCGA database (https://portal.gdc.cancer.gov/↗), GDSC database (https://www.cancerrxgene.org/↗), GeneMANIA (http://www.genemania.org↗), GSEA (http://software.broadinstitute.org/gsea/index.jsp↗), Human Protein Atlas database (https://www.proteinatlas.org↗), cBioPortal (http://www.cbioportal.org/↗), GDSC database (https://www.cancerrxgene.org/↗), GEPIA2 database (http://gepia2.cancer-pku.cn/#index↗), TISIDB (http://cis.hku.hk/TISIDB↗), TISCH database (http://tisch.comp-genomics.org/↗), ENCORI database (http://starbase.sysu.edu.cn/↗), RNAInter database (http://www.rnainter.org/↗), miRNet database (http://www.mirnet.ca/↗), RNA22 database (https://cm.jefferson.edu/rna22/interactive↗) as well as TIMER database (https://cistrome.shinyapps.io/timer/↗), UALCAN database (http://ualcan.path.uab.edu/↗).