The 365 HCC patients whose RNA sequencing data and relevant clinical information were obtained from The Cancer Genome Atlas (TCGA) database are regarded as the training cohort. The expression and information of 231 HCC patients were downloaded from the ICGC database as a validation cohort.

As presented in Supplementary Table 1 , the 1,136 mitochondria-related genes were extracted from MitoCarta3.0 (14). The differentially expressed genes (DEGs) were conducted by the “limma” R package in HCC and the relevant normal tissues, with the threshold of |log₂FC| >2 and FDR <0.05. The “VennDiagram” R package was used to identify prognostic mitochondria-related genes.

Univariate Cox regression was used to explore the potential prognostic mitochondria-related genes in HCC patients and the hazard ratio (HR) was shown. The mitochondria-based prognostic model was established by Lasso Cox regression and the R package “glmnet,” and 10-fold cross validation was employed to test the feasibility of the signature (15). The independent and dependent variables were candidate genes for the prognostic model and HCC patients’ survival times, respectively. According to the optimum parameter, the risk score was derived from the expression of candidate genes with standardized regression coefficients. The concrete formula is risk score = (MTHFD1L expression ∗ 0.109) + (NT5DC2 expression ∗ 0.156) + (POLQ expression ∗ 0.385) + (RECQL4 expression ∗ 0.017) + (TOMM40L expression ∗ 0.243) + (TXNRD1 expression ∗ 0.221). Then, HCC patients were classified into low- and high-risk groups because of the median risk score. Subsequently, the R packages “Rtsne” and “ggplot2” were used to examine the characteristics of principal component analysis (PCA) and t-distributed random neighbor embedding (t-SNE). The prognostic value of the model was estimated by the R packages “survminer” and “timeROC.” In addition, the predictive probability of the mitochondria-based prognostic model was observed from nomogram, which established a risk score and clinicopathological features using the R packages “rms” and “regplot.” The calibration curves and C-index value were introduced to verify the accuracy of the nomogram and risk score.

Multivariate Cox regression was investigated to judge the model’s independent prognostic value and the relationship between the risk score and clinical features (age, gender, grade, and tumor stage) in the TCGA and ICGC cohorts, and the visualized results were presented by the R package “forestplot.”

The R package “ClusterProfiler” was determined to the functional enrichment including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) between the low- and high-risk groups, with the threshold of Mann–Whitney test, |log2FC| ≥1, and FDR <0.05.

To identify the expression of immune cells and immune functions in low- and high-risk groups, we further conducted using the R package “GSVA” (16). The composition of 22 tumors associated with immune infiltrate cells was investigated by the R package “CIBERSORT” in risk groups. Spearman’s regression analyzed the infiltration of immune cells for the stromal score and immune score. The relationship between tumor stem cells and risk scores was evaluated by Spearman’s correlation. Moreover, we explored immune checkpoints and immune subtypes in low- and high-risk groups using the R package “ggpubr” and were tested by one-way ANOVA analysis.

According to the consensus cluster analysis, the R package “ConsensusClusterPlus” was used to estimate prognosis and immune characteristics in HCC patients (17). The R packages “Rtsne” and “ggplot2” were used to compare the characteristics of PCA and t-SNE among clusters. The Kaplan–Meier curve in the “survival” R package visualizes survival times in different clusters. The ESTIMATE score, immune score, and stromal score were obtained from each cluster’s R package “ESTIMATE.” In the TIMER database, we analyzed immune cell infiltration in HCC patient clusters based on various algorithms.

We collected 10 pairs of HCC and adjacent non-tumor tissues from Fuyang Hospital of Anhui Medical University, which were authorized by HCC patients and the Ethics Committee of the Hospital.

The human HCC cell lines (7721 and HepG2) and hepatic normal cells (LO2) were donated by the First Hospital Affiliated to Anhui Medical University. The cell culture conditions were as follows: DMEM with high glucose (HyClone), 10% fetal bovine serum (VivaCell, Shanghai, China), 5% CO₂, and 37°C.

The total RNA of tissues and cells was extracted using the TRIzol reagent (Takara). According to the instructions, reverse transcription was applied to the PrimeScript™ kit (Takara). Subsequently, the SYBR Green qPCR Mix was used to detect the relative expression of signature genes using the 2^−ΔΔCt method. The primer sequences used for the experiment are listed in Supplementary Table 2 .

Immunohistochemistry was used to detect the protein expression of six mitochondria-related genes. We analyzed the images of five-gene protein expression in “tissue” and “pathology” of modules in the HPA database. Patients and image serial numbers are included in this study. All images were re-judged by two pathologists. The regents and citations are as follows: MTHFD1L: Atlas Antibodies Cat#HPA029041, RRID: AB_2672880, dilution: 1:350; NT5DC2: Atlas Antibodies Cat#HPA050683, RRID: AB_2681215, dilution: 1:40; RECQL4: Atlas Antibodies Cat#HPA025821, RRID: AB_1856168, dilution: 1:40; TOMM40L: Atlas Antibodies Cat#HPA051304, RRID: AB_2681432, dilution: 1:75; TXNRD1: Santa Cruz Biotechnology Cat#sc-28321, RRID: AB_628405, dilution: 1:2500; POLQ: Atlas Antibodies Cat#HPA026762, RRID: AB_10601248, dilution: 1:450.

The one-way analysis of variance in HCC and normal tissues examined the DEGs. Cox regression with hazard ratios (HRs) and 95% confidence intervals (CIs) was employed to estimate prognosis. Those experimental results were processed using one-way ANOVA analysis using GraphPad Prism 8.0 software. The R software (version 4.0) was applied to statistical analyses, with p <0.05 considered significant.

As shown in Figure 1A , there were 106 prognostic genes differentially expressed in HCC and the corresponding normal tissues. Combined with 32 DEGs related to mitochondria in HCC, we screened eight mitochondria-related genes by Venn diagram ( Figure 1B ). The heatmap analysis suggested that the eight mitochondria-related genes were up-regulated in HCC tissues compared to normal tissues ( Figure 1C ). We identified the eight mitochondria-related genes as prognostic factors, which are regarded as hazards with a risk ratio greater than 1 ( Figure 1D ). Figure 1E presents the correlation among these genes.

We constructed a prognostic signature of six mitochondria-related genes using LASSO COX regression analysis, which was judged by optimum λ values. The risk score = 0.109 ∗ (MTHFD1L expression) + 0.156 ∗ (NT5DC2 expression) + 0.385 ∗ (POLQ expression) + 0.017 ∗ (RECQL4 expression) + 0.243 ∗ (TOMM40L expression) + 0.221 ∗ (TXNRD1 expression) ( Figures 1F , G ). In Table 1 , we suggested that the risk score was related to gender and tumor grade in the high-risk group of HCC patients. According to median risk scores, the patients with HCC were divided into low- and high-risk groups in the TCGA cohort ( Figure 2A ). As shown in Figure 2B , the PCA and t-SNE analyses, there are obvious clusters between low- and high-risk groups. The results of patients in the high-risk group show a higher risk score and shorter survival time ( Figure 2C ). The ROC curves of the 1-year, 2-year, and 3-year OS rates demonstrated the sensitivity and specificity of the model with values of 0.761, 0.685, and 0.677, respectively. Furthermore, the Kaplan–Meier analysis suggested the two groups was a statistical difference in survival in the TCGA cohort (P = 5.22e−04). The expression of MTHFD1L, NT5DC2, POLQ, RECQL4, TOMM40L, and TXNRD1 in the low-expression group had a longer survival time ( Supplementary Figure 1 ). In addition, 151 HCC patients were referred to as the validation set in the ICGC cohort. Based on the same median risk score as the TCGA cohort, the ICGC cohort was classified into low- and high-risk groups with 11 and 220 HCC patients, respectively ( Figure 2D ). The PCA and SNE analyses reconfirmed the difference in distribution between the two groups ( Figure 2E ). Moreover, compared with high-risk group, patients in low-risk groups had a better OS (P = 1.641e−04). The AUC value of the ROC analysis was 0.726 at 1-year, 0.700 at 2-year, and 0.727 at 3-year survival ( Figure 2F ).

To better evaluate the prognosis of HCC patients, we developed a Nomo diagram, including mitochondria-related gene signatures, grade, age, gender, risk score, and stage, to forecast the survival probability at 1, 3, and 5 years ( Figure 3A ). In Figure 3B , the OS predicted by nomogram was in good agreement with those observed OS in calibration curves. Moreover, the C-index demonstrated the prognostic model possessed greater superiority than other characteristics ( Figure 3C ).

To identify whether the mitochondria-related gene model acted as an independent prognostic indicator, we further investigated it using the multivariable Cox regression analysis. Figures 3D , E indicated that the risk score was observably related to OS in the TCGA cohort (P <0.001, HR = 2.718 (1.502–4.921)) and ICGC cohort (P <0.001, HR = 2.680 (1.874–3.832)). Combined with clinicopathological characteristics, there were significantly different tumor stages I–II and tumor stage III–IV in patients in the TCGA cohort (P = 0.042) and the ICGC cohort (P = 0.008) ( Figures 3F –L ). According to the TCGA cohort, the risk score was also correlated with tumor grade (P = 3.2e−06). The above results demonstrated that the mitochondria-based six-gene model could be an adverse prognostic model for patients with HCC.

To evaluate the immune activity of the mitochondria-related gene signature, we first investigated the immune cells and immune functions using ssGSEA in low- and high-risk groups. As shown in Figures 4A , B , the nine immune cells and seven immune functions were statistical significance between two groups in the TCGA cohort. However, in the ICGC cohort, only the two (NK_cells and Th2_cells) immune cells were associated with risk groups ( Figures 4C , D ). Compared with the high-risk group, the activity of the Type_II IFN response was upregulated in the low-risk groups.

Then, according to 22 types of immune cells, we draw a detailed immune cell compositions in HCC. The CD8+ T cells, M0 macrophages, and CD4+ T cells naive were the prime cell compositions in the high-risk group, whereas plasma cell in the low-risk group ( Figure 4E ). We further analyzed tumor stemness, immune, and stromal cell scores to illuminate the role of the mitochondria-related gene signature in the tumor microenvironment. The results indicated that risk score was positively correlated with RNA stemness score (R = 0.28, P = 1.2e−07), while stromal scores were negatively correlated (R = −0.160.28, P = 0.003) ( Figure 4F ). Moreover, the immune checkpoint analysis suggested that almost all the relevant genes in high-risk groups were significantly higher than those in low-risk groups ( Figure 4G ).

Subsequently, we explored the relationship between risk groups and partial immune checkpoints (CTLA4, CXCR2, GPC3, HAVCR2, PDCD1, and PDL1). As shown in Supplementary Figure 2 , the partial immune checkpoints except PDL1 were elevated significantly in the high-risk group. The expression of these immune checkpoints was positively correlated with the risk score. In addition, based on TCGA, the HCC patients were classified into four immune subtypes. The C3 and C4 subtypes are the prime immune subtypes, with a proportion of 82%. Most patients in the low-risk group are C3 (n = 80, 45%), C4 (n = 82, 46%) subtypes, while C4 (n = 77, 44%) is less common ( Figure 4H ).

With the expression profiles of six mitochondria-related genes, we regrouped HCC patients into three clusters (cluster 1, cluster 2, and cluster 3) by consensus clustering analysis ( Figure 5A ). We verified the distributions of the above three clusters using PCA and tSNE analysis, implying that these clusters were distinguished ( Figure 5B ). Next, Kaplan–Meier analysis suggested that the survival of cluster 1 was longer than that of other clusters ( Figure 5C ). Moreover, we visualized the relationship between clusters and risk groups in an alluvial diagram ( Figure 5D ). Cluster 1 was primarily associated with the low-risk group, while clusters 2 and 3 were related to the high-risk group.

Concerning immunity, clusters 2 and 3 were closely associated with immune cells infiltrating ( Figure 5E ). The ESTIMATE score indicated there were no significant differences among the three clusters ( Figure 5F ). Cluster 1 had a higher stromal score compared to cluster 2. Cluster 3 had a higher immune score, signifying a more diverse tumor microenvironment than cluster 1. Besides, the expression of immune checkpoints was more abundant in clusters 2 and 3, indicating that they were susceptible to immunotherapy ( Figure 5G ).

To explore potential biological functions, we carried out GO and KEGG analyses in low- and high-risk groups. In Figures 6A , B , we found that the richest biological processes were related to organelle fission, nuclear division, mitotic nuclear division, and chromosome segregation. The first three cellular components were the spindle, the chromosome region, and the condensed chromosome. Moreover, in terms of molecular functions, microtubule binding, oxidoreductase activity (acting on paired donors with incorporation or reduction of molecular oxygen and acting on the CH-OH group of donors), and heme binding. In addition, the KEGG analysis showed that the risk groups jointly participated in the cell cycle, metabolism of xenobiotics, drug metabolism, oocyte meiosis, lipid atherosclerosis, and retinol metabolism. The above results concluded that its functions were closely associated with oxidation, metabolism, and the cell cycle ( Figures 6C , D ).

To further improve the therapeutic effect on HCC patients, we investigated the relationship between risk groups and clusters and drug sensitivity. We identified significant differences in the IC50 of drugs between low- and high-risk groups, including Imatinib, Tipifarnib, Epothilone, Doxorubicin, and Bosutinib ( Supplementary Figures 3A–E ). Moreover, the results of IC50 associated with various drugs among clusters were consistent with those in the risk groups ( Supplementary Figures 3F–K ). Therefore, considering the mitochondria-related model, the prognosis of HCC patients is expected to be greatly improved in the future.

The protein expression of 6-gene was analyzed by immunohistochemistry in the HPA database. In Figures 7A–F , we found that the 6-gene protein were positively expression in HCC tissues compared to normal tissues. Moreover, MTHFD1L, NT5DC2, POLQ, TOMM40L, and TXNRD1 were in cytoplasmic or membranous, while RECQL4 was in nuclear. Based on the above results, we further investigated the expression of six mitochondria-related genes by RT-qPCR analysis. Compared with normal tissues, MTHFD1L, NT5DC2, POLQ, RECQL4, TOMM40L, and TXNRD1 expression were upregulated in the 10 pairs of HCC tissues ( Figures 7G –L ). In terms of HCC cell lines (7721 and HepG2), these six genes were identified with the results of HCC tissues. So, the mitochondria-related genes could be a potential biomarker for HCC ( Figures 7M –R ).