Characterization of immune features and immunotherapy response in subtypes of hepatocellular carcinoma based on mitophagy

Mitophagy is suggested to be involved in tumor initiation and development; however, mitophagy heterogeneity in hepatocellular carcinoma (HCC) and its association with immune status and prognosis remain unclear. Differentially expressed genes (DEGs) were identified using expression profiles acquired from The Cancer Genome Atlas (TCGA). Mitophagy-related subtypes were identified using the ConsensusClusterPlus software. The differences in prognosis, clinical characteristics, and immune status, including immune cell infiltration, immune function, immune-checkpoint gene expression, and response to immunotherapy, were compared between subtypes. A mitophagy-related gene signature was constructed by applying least absolute shrinkage and selection operator regression to the TCGA cohort. The International Cancer Genome Consortium cohort and the cohort from Peking Union Medical College Hospital were utilized for validation. Carbonyl cyanide m-chlorophenylhydrazone was used to induce mitophagy in HCC cell lines to obtain our own mitophagy signature. Real-time polymerase chain reaction was used for the experimental validation of the expression of model genes. Two mitophagy-related subtypes with distinct prognoses, clinical characteristics, immune states, and biological function patterns were identified based on the mitophagy-related DEGs. The subtype that showed higher mitophagy-related DEG expression had worse survival outcomes, suppressed immune function, higher immune-checkpoint gene expression, and a better response to immunotherapy, indicating that this subpopulation in HCC may benefit from immune-checkpoint blockade therapy and other immunotherapies. A risk model consisting of nine mitophagy-related genes was constructed and its performance was confirmed in two validation cohorts. The risk score was an independent risk factor even when age, sex, and tumor stage were considered. Our study identified two distinct mitophagy subtypes and built a mitophagy signature, uncovering mitophagy heterogeneity in HCC and its association with immune status and prognosis. These findings shed light on the treatment of HCC, especially with immunotherapy.

Introduction

Hepatocellular carcinoma (HCC) accounts for the majority of primary liver cancer cases, which is ranked the fifth in cancer-related death (1, 2). Despite much progress in the diagnosis and treatment of HCC, the prognosis of patients with HCC remains poor, with a median survival time of 9 months (3). For patients with HCC at early stage, curative treatments such as radiofrequency ablation and liver section can achieve a 40%-70% 5-year survival rate; and palliative treatments such as transarterial chemoembolization has been shown to improve median OS of intermediate stage HCC to approximately 20 months (4). For HCC at advanced stage or terminal stage, survival outcomes are still unsatisfactory even with the help of molecular therapy. Immunotherapies, such as immune-checkpoint blockade (ICB), have shown strong antitumor activity and lead to a substantial prolonged survival for advanced HCC, whereas only a subset of patients can benefit from these therapies (5). Therefore, there is an urgent need to explore the underlying molecular mechanisms of HCC and provide new targets and strategies for treatment.

Major breakthroughs in a mechanism called mitophagy have recently gained considerable attention (6). Mitophagy, also known as mitochondrial autophagy, eliminates denatured or damaged mitochondria, preventing the accumulation of mitochondrial DNA mutations and maintaining mitochondrial quality (7). Hence, mitophagy plays a vital role in regulating energy metabolism and removing excessive cytotoxic reactive oxygen species (8). Mitophagy plays a dual role in the development of cancer by suppressing tumors at an early stage and promoting tumors at an advanced stage (9). The ubiquitin-dependent PINK1/Parkin pathway is the most common mitophagy cascade, and some core genes within this pathway, such as PINK1 and PARK2, can predict prognosis in patients with papillary renal cell cancer (10, 11). However, the role of mitophagy-related genes (MRGs) in HCC is not fully understood. Several studies have reported heterogeneity in autophagy in other types of cancer (12, 13). As a specific form of autophagy, the heterogeneity of mitophagy likely influences the development and prognosis of HCC. Research focused on mitophagy may help to concretize problems. Therefore, we aimed to investigate the role of MRGs and mitophagy-related subtypes in HCC, focusing on the association with immune status and response to immunotherapy, as there has been massive interest in immunotherapy for HCC (14).

A flowchart of the study design is shown in Figure 1. In this study, we first screened differentially expressed MRGs (DEMs) between tumor and normal tissues of patients with HCC. Based on the expression profile of the DEMs, we classified the patients into two subtypes and explored their prognoses, clinical characteristics, immune states, and drug sensitivities. Subsequently, based on the MRG signature, a prognostic model was constructed and validated in two cohorts. Moreover, we explored the differences in biological functions between these subtypes and risk groups. Cell experiment and qPCR were performed to validate our results. Our findings are helpful in accurately characterizing HCC and providing personalized treatment for patients.

FIGURE 1

Figure 1 Study flowchart. DEMs, Differentially expressed mitophagy-related genes; ICB, immune-checkepoint blockade.

Materials and methods

Data acquisition

We obtained two gene sets from public databases by searching the keyword “mitophagy”: one is Mitophagy-animal pathway (Entry: hsa04137), which contains 72 genes from the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (https://www.kegg.jp/kegg/pathway.html), the other one is REACTOME-MITOPHAGY gene sets (source: R-HSA-5205647), which contains 29 genes from the C2:CP : REACTOME in Molecular Signatures Database with Gene Set Enrichment Analysis (GSEA) (http://www.gsea-msigdb.org). Since 20 genes overlapped in two gene sets, we eventually acquired a total of 81 MRGs for subsequent analyses. The RNA-seq and clinical information of HCC samples were downloaded from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/), which included 50 normal samples and 374 cancer samples. We also acquired two HCC cohort datasets for validation: one was downloaded from the International Cancer Genome Consortium (ICGC) database (https://dcc.icgc.org/), which included 243 cancer samples; the other was collected from Peking Union Medical College Hospital (PUMCH) and included 20 patients with HCC (Supplementary Table 1). The cohort from our center was approved by the Ethics Committee of PUMCH and CAMS (Chinese Academy of Medical Sciences) & PUMC (Peking Union Medical College), and written informed consent was obtained from all patients.

Identification and analysis of DEMs

DEMs between tumor tissues and normal tissues in TCGA were screened using the “limma” R package (15) based on the following criteria: |log₂fold change| > 0.5 and false discovery rate < 0.05. The protein-protein interaction network of DEMs was obtained from the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (https://string-db.org), and the interaction between core genes was visualized using Cytoscape software (version 3.8.2) (16).

Consensus cluster analysis

We used the ConsensusClusterPlus software (17) of R to perform unsupervised consensus clustering of TCGA dataset based on the expression of DEMs. The optimal cluster number k was determined by evaluating the delta area, consensus cumulative distribution function, and consensus matrix. Principal component analysis was used to verify the results of the cluster analysis. The correlation between clusters and clinical variables was tested using Chi-square test.

Immune status analysis

To explore the impact of mitophagy on patient immune status, two mitophagy-related subtypes were compared in terms of differences in infiltrating immune cells, immune function, immune-checkpoint gene expression levels, and response to immunotherapy. We quantified the relative abundance of immune cell types and the activity of immune function in each sample using single sample GSEA algorithms through the R package “GSVA” (18). The expression levels of immune-checkpoint genes can reflect the response to ICB treatment; thus, the following well-known immune-checkpoint genes were chosen for expression level evaluation in each subtype: CTLA4, CXCL9, CD8A, TBX2, PDCD1, LAG3, HAVCR2, IFNG, TNF, and CD274. Moreover, Tumor Immune Dysfunction and Exclusion (TIDE) scores (19) of each HCC sample were calculated online (http://tide.dfci.harvard.edu/) and compared between subtypes to verify the differences in response to immunotherapy.

Drug sensitivity analysis

To discover potential drugs for patients with different mitophagy-related subtypes of HCC, we evaluated their responses to various antitumor drugs using the “pRRophetic” R package (20), which is based on the Genomics of Drug Sensitivity in Cancer database.

Construction and validation of risk model

First, we performed univariate Cox regression analysis on MRGs to screen for survival-related prognostic genes in the TCGA cohort. We then obtained genes for model construction by intersecting prognostic genes with DEMs, followed by least absolute shrinkage and selection operator (LASSO) regression using the “glmnet” R package (21) to form the final gene signature for the risk model. The risk score was formulated as follows:

$$\mathrm{risk score = \beta 1 \times gene1 expression + \beta 2 \times gene2 expression + \dots + \beta n \times gene expression}$$

where β represents the coefficient value of each gene. Patients in the training and validation cohorts were divided into high- and low-risk groups based on the median risk scores. The receiver operating characteristic (ROC) curves in each cohort were plotted using the R package “timeROC” (22), and the time-dependent area under the curve values were measured to evaluate the performance of the model. Univariate and multivariate Cox regression analyses were used to evaluate whether the risk score was an independent predictor.

Functional analysis

KEGG and Gene Ontology (GO) enrichment analyses were utilized for functional annotation of DEMs and differentially expressed genes (DEGs) between subtypes using the “ClusterProfiler” R package (23). Significant GO terms and pathways were selected with a p-value cutoff of < 0.01. The biological functions enriched in the high- and low-risk groups of the TCGA cohort were derived by GSEA of KEGG pathways.

Acquisition of mitophagy signature by inducing mitophagy in HCC cell lines

Human HCC cell line Huh 7 was purchased from Procell Life Science & Technology Co. Ltd. (Wuhan, China). MHCC97H was from the Liver Cancer Institute (Zhongshan Hospital, Fudan University, China). SNU398 was from ATCC (Manassas, VA, USA). Huh7 and MHCC97H were cultured in Dulbecco’s modified Eagle medium, supplemented with 10% fetal bovine serum and 1% antibiotics in 5% CO₂ at 37°C. SNU398 was routinely cultured in RPMI-1640 medium supplemented with 10% fetal bovine serum and 1% antibiotics in 5% CO₂ at 37°C. Then three HCC cell lines were treated with 10μM carbonyl cyanide m-chlorophenylhydrazone (CCCP, selleck,China), which was used to induce mitophagy, for 24 hours (24, 25). RNA sequencing for Huh 7, MHCC97H, SNU398 cells treated with or without CCCP by Beijing Auwigene Tech, Ltd (Beijing, China) using the Illumina second-generation high-throughput PE150 sequencing platform (Illumina, Inc., CA, United States). Between cell lines treated with and without CCCP, top 100 differentially expressed genes ranked by |log₂fold change| were considered as mitophagy signature for validation.

Real-time polymerase chain reaction

Total RNA was extracted from liver tumors and peritumoral tissues using TRIzol (Invitrogen, Thermo Fisher, Waltham, MA, USA), following the manufacturer’s instructions. RNA (1 µg) was reverse transcribed using the Hifair^® II 1st Strand cDNA Synthesis SuperMix for qPCR (gDNA digester plus) (Yeasen Biotechnology, Shanghai, China) in a 20 μl reaction. After 20-fold dilution, 4 μl of cDNA was used as a template in a 20 μl real-time polymerase chain reaction (PCR). For real-time PCR, amplification was performed for 40 cycles using BlasTaq™ 2X qPCR MasterMix (Applied Biological Materials Inc., Richmond, BC, Canada). Primers were designed on exon junctions to prevent co-amplification of genomic cDNA; the sequences are presented in the Supplementary Table 2.

Statistical analyses

All statistical analyses were performed using R (version 4.0), and relevant packages were applied for processing and visualization. The Wilcoxon test was used to compare differences in continuous variables. The overall survival was evaluated using the Kaplan–Meier analysis, and the survival curve was plotted using the R package “survminer” (http://cran.r-project.org/). The log-rank test was used to examine the differences between subtypes or groups. If not specifically stated, bilateral p < 0.05 was considered statistically significant.

Results

Profile and functional annotation of DEMs between tumor and normal samples

We collected 81 MRGs from the database gene sets, and 51 DEMs were identified between tumor and normal samples from the TCGA cohort. In order to verify DEMs acquired from one public database, we performed RNA sequencing on 8 pairs of HCC samples and peritumoral tissues, and the expression matrix can be found in Supplementary Table 1. Of 42 DEMs identified in our own samples, 35 DEMs overlapped with 51 DEMs from TCGA cohort, and 39 of 40 DEMs identified in ICGC cohort overlapped with those in TCGA cohort, indicating the reliability of DEMs identified in TCGA cohort. A heatmap of the 51 DEMs is shown in Figure 2A. Fifty of the DEMs were upregulated in tumors, primarily including genes involved in the PINK1/Parkin pathway (PINK1, PARK2, ATG family, and TOMM family) and receptor-mediated mitophagy (FUNDC1, PGAM5, and ULK1). In addition, oncogenes such as TP53, KRAS, and HRAS were also upregulated, as these genes may be related to hypoxic stress. Only one gene, JUN, was downregulated. We then performed protein-protein interaction analysis to determine the interactions between DEMs and the core network, as shown in Figure 2B; the TOMM family, MFN1, VDAC1, PARK2, and RPS27A were hub genes, and oncogenes such as HRAS, KRAS, and TP53 participated in the core network of mitophagy. The KEGG analysis in Figure 2C showed that, besides mitophagy and autophagy, these DEMs were also enriched in the PD-1/PD-L1 checkpoint pathway, which is related to the response to ICB treatment of HCC. Pathways involving hepatitis B and apoptosis were also enriched, and they were shown to be involved in tumorigenesis and the development of HCC (26, 27).

FIGURE 2

Figure 2 Expression level, interactions and functional enrichment analysis of DEMs between tumor samples and normal samples. (A) Profile of DEMs based on sample type. Color represents expression level (blue to red). (B) Hub protein-protein interaction network among DEMs. Color represents confidence (blue to red). (C) KEGG analysis of DEMs. DEMs: Differentially expressed mitophagy-related genes, KEGG: Kyoto Encyclopedia of Gene and Genome.

Subpopulation of HCC on account of expression pattern of DEMs

Consensus clustering was applied to identify HCC subtypes based on the expression levels of the DEMs acquired from the previous step. We determined the k value as 2, at which point the relative change in area under the cumulative distribution function reached an approximate maximum and the consensus matrix showed a clear boundary simultaneously (Figures 3A–C). Therefore, two clearly distributed subtypes were classified; these were denoted as cluster 1 (containing 211 samples) and cluster 2 (containing 163 samples). To further verify the clustering result, principal component analysis was performed, and the principal component distribution was in accordance with the consensus matrix, ensuring the stability of consensus clustering (Figure 3D). The DEM expression and clinical features of each sample grouped by cluster are shown in a heat map (Figure 3E). Cluster 2 generally had higher DEM expression than cluster 1. Tumor stage and grade were found to be correlated with subtype. Cluster 2 had a higher proportion of tumors with advanced stage and high grade. Moreover, the Kaplan–Meier survival analysis showed that the two subtypes had significant differences in OS (Figure 3F). Cluster 2 tended to have worse outcomes than cluster 1 (p = 2.618e-04), with 5-year survival rates of 37.6% and 56.0%, respectively. The difference in prognosis between the two clusters was in accordance with their differences in tumor stage and grade. These findings confirm the existence of mitophagy heterogeneity in HCC and its impact on the development and prognosis of HCC.

FIGURE 3

Figure 3 Identification of two mitophagy-related HCC subtypes with different prognosis. (A–C) Consensus matrix of HCC samples co-occurrence proportion for k = 2 (A), relative change in area under the CDF curve for k from 2 to 7 (B), consensus clustering CDF for k from 2 to 7 (C). (D) Principal component analysis of HCC samples grouped by subtype. (E) A heatmap showing the association of mitophagy-related subtypes with clinical characteristics and expression of DEMs. (F) The Kaplan–Meier plot showing the overall survival differences between the two subtypes. The asterisks represent the p value (*p < 0.05; **p < 0.01; ***p < 0.001). DEMs, Differentially expressed mitophagy-related genes. CDF, Cumulative distribution function.

Characterization of immune status and drug sensitivity affected by mitophagy heterogeneity in HCC

To compare the immune characteristics of the two mitophagy-related subtypes, we first estimated immunocyte infiltration and immune function using single sample GSEA algorithms. As shown in Figure 4A, compared with cluster 1, cluster 2 showed higher infiltration of activated dendritic cells (aDCs), immature dendritic cells (iDCs), macrophages, follicular helper T cell (Tfh), T helper 2 cell (Th2) and regulatory T cells (Treg). Regarding immune function in Figure 4B, cytolytic activity and type II interferon response were increased in cluster 1, while cluster 2 had strong antigen-presenting cell (APC) co-stimulation and major histocompatibility complex (MHC) class I reactions.

FIGURE 4

Figure 4 The comparison of immune status between mitophagy-related subtypes. Box plots showing the differences of infiltrating immunocyte abundance (A), immune reaction activity (B), expression of immune-checkpoint genes (C) and violin plots showing Tumor Immune dysfunction and Exclusion (TIDE) score (D). The asterisks represent the p value (*p < 0.05; **p < 0.01; ***p < 0.001). ns, not significant.

To explore the effect of mitophagy on the response to ICB treatment, we compared the expression levels of immune-checkpoint genes in each subtype. As shown in Figure 4C, all immune-checkpoint genes were consistently overexpressed in cluster 2, indicating that cluster 2 tended to be more sensitive to ICB treatment. Furthermore, we calculated the TIDE score of every sample and the scores were significantly lower in cluster 2 than in cluster 1, further verifying that patients in cluster 2 may be more likely to benefit from immunotherapy (Figure 4D). In contrast, with higher TIDE scores, cluster 1 was more likely to achieve tumor immune escape and exhibit a lower response rate to ICB treatment.

We also evaluated the drug sensitivity of each subtype to identify potential chemotherapeutic drugs. Lower IC50 values indicate higher sensitivity. As shown in Figure 5, compared with cluster 2, cluster 1 was more sensitive to AKT inhibitor III, epidermal growth factor receptor inhibitors such as erlotinib, gefitinib, and lapatinib, and vascular endothelial growth factor receptor inhibitors such as axitinib and sunitinib. Conversely, cluster 2 had a higher response rate to AZD8055 (mTOR inhibitor), bleomycin, cisplatin, etoposide, sorafenib, and methotrexate.

FIGURE 5

Figure 5 The difference of chemo drugs sensitivity between subtypes, including metformin (A), sorafenib (B), sunitinib (C), methotrexate (D), AKT.inhibitor.VIII (E), bleomycin (F), axitinib (G), AZD8055 (H), etoposide (I), lapatinib (J), cisplatin (K), erlotinib (L).

Functional annotation of DEGs between subtypes

To reveal the differences in biological functions between the two subtypes, we conducted GO and KEGG enrichment analysis on the DEGs between the two subtypes with a cutoff of |log₂fold change|> 1 and false discovery rate < 0.05. A total of 260 DEGs met the criteria, and the results showed that complement and coagulation cascades, the peroxisome proliferators-activated receptors signaling pathway, bile secretion, chemical carcinogenesis, and drug and compound metabolism were significantly enriched in the KEGG pathway analysis (Figure 6A). The results of GO enrichment analysis are shown in Figure 6B.

FIGURE 6

Figure 6 Functional enrichment analysis of DEGs between two mitophagy-related subtypes. Bar plots showing the biological function of DEGs using KEGG (A) and GO (B) enrichment. DEGs, Differentially expressed genes; KEGG, Kyoto encyclopedia of genes and genomes; GO, Gene ontology.

Constructing prognosis model of MRGs

Defining the TCGA dataset as a training cohort, we performed univariate Cox regression analysis on 81 MRGs, 35 of which were significantly associated with the OS of patients with HCC (p < 0.05). After intersection with 51 DEMs, we obtained 24 genes for model construction (Figure 7A). All 24 genes were risk genes with a hazard ratio of > 1 (Figure 7B). The LASSO regression model was then utilized, and nine genes were screened to build the prognostic risk model (Figures 7C, D). The risk score was calculated using the corresponding coefficients and gene expression. Finally, the risk score model was formulated as follows:

$$\begin{array}{l}ATG9A\times 0.2414855+ATG12\times 0.0294790\\ +HRAS\times 0.0977971+MFN1\times 0.1418254\\ +NRAS\times 0.2199555+PGAM5\times 0.2023505\\ +SQSTM1\times 0.1721932+TOMM22\\ \times 0.1709305+TOMM5\times 0.1709305\end{array}$$

FIGURE 7

Figure 7 Construction of a LASSO regression model and correlation between subtypes and risk groups. (A) Venn diagram showing intersection between DEMs and prognostic genes. (B) Forest plots showing the results of univariate Cox regression analysis of overlapped genes. (C, D) LASSO regression analysis of the overlapped genes. (E) The Sankey diagram showing the distribution of patients in mitophagy-related subtypes, risk groups and survival outcomes. LASSO, Least absolute shrinkage and selection operator.

Additionally, we created a Sankey diagram to show the connection among mitophagy subtypes, risk scores, and survival (Figure 7E).

Validation of model efficiency

To verify the performance of the risk model, we performed Kaplan–Meier survival analysis in the training and validation cohorts. The survival curves showed that improved survival rates of low-risk patients continued for nearly 7 years in the TCGA training cohort (p = 9.707e-04), and this advantage existed in the ICGC validation cohort (p = 1.749e-04) (Figures 8A, B). In addition, we used an HCC cohort (n=20) registered in our center to validate the risk model, and the difference in OS was still significant (p = 3.924e-02) (Figure 8C). Regarding model accuracy, the 1-year, 3-year, and 5-year AUC of the model for OS was 0.781, 0.690, and 0.650, respectively, in the TCGA training cohort (Figure 8D), and 0.709, 0.749, and 0.716, respectively, in the ICGC validation cohort (Figure 8E). The AUC of the model in PUMCH cohort was still satisfactory (Figure 8F). We ranked the risk scores of patients with HCC in all cohorts from low to high, and the survival status and time of each patient were plotted according to the risk score (Figures 8G–L). The plot revealed that high-risk patients generally had poorer survival rates than low-risk patients.

FIGURE 8

Figure 8 The risk model performance in training cohort and two validation cohorts. (A–C) Kaplan-Meier curves for the OS of patients in the high- and low-risk group in the TCGA cohort (A), ICGC cohort (B) and PUMCH cohort (C). (D–F) AUC of time-dependent ROC curves in the TCGA cohort (D), ICGC cohort (E) and PUMCH cohort (F). (G–I) The distribution and median value of the risk scores in the TCGA cohort (G), ICGC cohort (H) and PUMCH cohort (I). (J–L) The distributions of risk scores, survival states and survival outcomes in the TCGA cohort (J), ICGC cohort (K) and PUMCH cohort (L). (M, N) Forest plots showing the univariate (M) and multivariate (N) Cox regression analyses regarding OS in the TCGA cohort. OS, Overall survival, AUC, area under the curve, ROC, receiver operating characteristic, TCGA, The Cancer Genome Atlas, ICGC, International Cancer Genome Consortium, PUMCH, Peking Union Medical College Hospital.

To determine whether the risk score is an independent risk factor for the prognosis of patients with HCC, we performed univariate Cox regression on the risk score and clinical variables (Figure 8M). The results showed that only the stage and risk scores were significantly associated with OS (p < 0.001). Next, these variables were included in the multivariate Cox regression analysis. After correction for other confounding factors, including age, sex, stage, and grade, the risk score was still significantly associated with OS, implying that the risk score was an independent risk factor (p < 0.001) (Figure 8N).

Functional enrichment analysis based on the risk score

We performed GSEA on the TCGA cohort, and the most significantly enriched KEGG pathways are shown in Figure 9. The cell cycle, mTOR signaling pathway, NOTCH signaling pathway, endocytosis, and pathways in cancer were enriched in the high-risk group. Primary bile acid biosynthesis, drug metabolism, cytochrome P450, fatty acid metabolism, glycine serine and threonine metabolism, and linoleic acid metabolism pathways were enriched in the low-risk group.

FIGURE 9

Figure 9 Gene-set enrichment analysis identifying KEGG pathways enriched in the high- and low-risk group. KEGG, Kyoto Encyclopedia of Gene and Genome.

Validation of mitophagy heterogeneity through cell experiment

In order to validate mitophagy-related subtypes obtained using public mitophagy gene sets, we acquired mitophagy signature through inducing mitophagy in HCC cell lines for same analyses in TCGA HCC dataset. The expression matrix of HCC cell lines before and after mitophagy induction was demonstrated in Supplementary Table 1 and the expression levels of mitophagy signature genes were applied for clustering of TCGA HCC dataset. As shown in Figure 10A, tumor grade and stage were still correlated with clusters. And cluster 2 had significantly worse survival outcome than cluster 1 (p = 0.006) (Figure 10B). Similar to results in Figure 4A, cluster 2 had higher infiltration of aDCs, iDCs, macrophages, Tfh, Th2, and Treg than cluster 1 (Figure 10C). In terms of immune function, type II interferon response was still suppressed in cluster 2, while check-point, APC co-stimulation, APC co-inhibition, HLA, MHC class I, and proinflammation exhibited higher levels in cluster 2 (Figure 10D), which confirmed the association between mitophagy and immune status in HCC. In addition, regarding response to ICB treatment, all immune-checkpoint genes except CXCL9 were significantly overexpressed in cluster 2 (Figure 10E), and TIDE scores remained lower in cluster 2 than in cluster 1 (Figure 10F), indicating that patients in cluster 2 tended to benefit from ICB treatment. These findings following cell experiment further validated our results from using public MRGs.

FIGURE 10

Figure 10 The validation of mitophagy-related HCC subtypes using mitophagy signature obtained from cell experiment. (A) A heatmap showing association between subtypes and clinical characteristics. (B) The Kaplan–Meier plot showing distinct prognosis between two subtypes. (C–E) Box plots showing the differences of infiltrating immunocyte abundance (C), immune reaction activity (D), and expression of immune-checkpoint genes (E) between two subtypes. (F) Violin plots comparing the Tumor Immune dysfunction and Exclusion (TIDE) score of each subtype. The asterisks represent the p value (*p < 0.05; **p < 0.01; ***p < 0.001). ns, not significant.

Validation of expression of model genes in tissue

To verify the reliability of the results acquired from the public database, we further validated the expression levels of the nine genes consisting of the mitophagy signature in five pairs of HCC tissues and peritumoral tissues. As shown in Figure 11, the tumor expressed significantly higher mRNA levels of all genes (ATG9A, ATG12, HRAS, MFN1, NRAS, PGAM5, SQSTM1, TOMM22, and TOMM5) than peritumoral tissue, which was consistent with the public database analysis.

FIGURE 11

Figure 11 The experimental validation of nine genes consisting the risk model using Real-Time PCR, including TOMM5 (A), SQSTM1 (B), TOMM2 (C), PGAM5 (D), NRAS (E), MFN1 (F), ATG9A (G), ATG12 (H), HRAS (I). The asterisks represent the p value (*: p < 0.05; **: p < 0.01; ***: p < 0.001). PCR, polymerase chain reaction.

Discussion

Emerging immunotherapy, especially ICB treatment, has become an effective and promising option to treat HCC (28). However, only a portion of patients respond to immunotherapy; thus, it is important to determine which groups of patients can benefit from immunotherapy, facilitating the progress of personalized treatment. Recently, mitophagy has attracted the attention of researchers as a potential therapeutic target for cancer. Hence, this study aimed to investigate mitophagy heterogeneity in HCC and its association with immune status, identify two mitophagy subtypes with distinct clinical and immune characteristics, and offer more detailed insights into immunotherapy or combination therapy for HCC.

HCC has been confirmed to exhibit high molecular heterogeneity (29). We identified two mitophagy subtypes in TCGA HCC samples based on the expression levels of DEMs, showing that these two subtypes with different mitophagy patterns were characterized by significantly different tumor stages and prognoses. This verified that mitophagy heterogeneity is associated with HCC development and has prognostic value in HCC, although the underlying mechanisms are still not well understood. Mitophagy appears to be tumor-promoting or tumor-suppressive, depending on the tumor type and intrinsic stage (30). For instance, PARK2, which encodes a core mitophagy protein, Parkin, was found to be inactivated in colon and lung cancer (31). Parkin-null mice are susceptible to spontaneous HCC (32). In contrast, upregulation of the mitochondrial inner membrane protein STOML2 can amplify PINK1/Parkin-mediated mitophagy and facilitate the migration and invasion of HCC cells, thus promoting HCC growth and metastasis (33). Previous studies have also found that hyperactivated mitophagy can induce sorafenib resistance in HCC under hypoxic stress (34). Therefore, the dual role of mitophagy may be involved in HCC heterogeneity. In our study, nearly all DEMs were upregulated in HCC samples compared with normal tissues, and cluster 2 had generally higher expression levels of DEMs but worse prognosis than cluster 1. Based on the above evidence, cluster 2 is likely to be characterized by higher mitophagy activity, which results in a more advanced tumor and worse survival outcome. Furthermore, regulation of various mitophagy pathways, such as PINK1/Parkin-mediated mitophagy and BNIP3/BNIP3L/FUNDC1-mediated mitophagy, may also be involved in HCC heterogeneity, which warrants further study.

To determine whether mitophagy heterogeneity has an impact on the tumor immune microenvironment, we evaluated the differences in immune characteristics between the two subtypes. Immune cell infiltration is closely related to clinical outcomes, and immune cells can serve as an immunotherapy target (35). Our single sample GSEA results indicated that cluster 2 had a higher abundance of regulatory T cells and macrophages, which are considered to be HCC promoting (36). Moreover, cluster 2 was characterized by higher expression levels of immune-checkpoint genes. Overexpression of immune-checkpoint genes can suppress the antitumor immune response so that tumor cells can easily evade immune surveillance. These findings explain the poor survival outcomes of cluster 2.

ICB therapy can restore dysfunctional immune system and has achieved remarkable results in cancer treatment. ICB agents against programmed cell death protein 1 (PD-1) and cytotoxic T lymphocyte antigen 4 (CTLA-4) have been approved for HCC by the FDA (37). However, the limited response rate makes it important to screen patients who are sensitive to ICB therapy. Cluster 2 showed higher expression levels of immune check-point genes, indicating a better response to ICB treatment. The TIDE algorithm is believed to perform better than the expression level of immune check-points in predicting the survival outcome of cancer patients treated with ICB agents (19). Corresponding to the prediction based on immune-checkpoint expression, the TIDE results revealed that cluster 2 was more likely to respond to ICB treatment. Therefore, ICB treatment may help reverse the poor prognosis of cluster 2. Taken together, mitophagy heterogeneity in HCC may influence immune status and can predict the response rate to ICB agents, revealing the association between mitophagy and immunity. This result enhanced our understanding of the heterogeneity of HCC, promoting personalized therapy in clinical practice and inspiring immunotherapy development in scientific research and trials. Furthermore, our study explored potential drugs for subpopulations with different mitophagy patterns, providing ideas for synergistic combination of ICB agents and targeted therapies. Systemic therapy in HCC should be explored to improve clinical efficacy (38).

To more precisely predict the prognosis of patients with HCC, we constructed a risk model based on a mitophagy signature. Notably, all nine genes were risk factors for HCC. Of these genes, ATG9A and ATG12 are core regulators of autophagy (39). MFN1, also known as mitofusin-1, was analyzed both in vivo and in vitro and its effects on HCC metastasis were revealed (40). PGAM5 is an atypical mitochondrial serine/threonine phosphatase that dephosphorylates FUNDC1 to activate mitophagy. Previous studies have reported that depleting PGAM5 inhibits tumor development and enhances the 5-fluorouracil sensitivity of HCC cells (41, 42). The TOMM complex (translocase of the outer mitochondrial membrane) imports nearly all mitochondrial proteins from the cytoplasm into the mitochondria, and TOMM22 functions as a central receptor (43). Although no significant increase in the expression of TOMM genes was observed in prostate cancer compared to normal tissues, our results demonstrated that this protein was elevated in HCC and may be a good candidate biomarker; this requires validation (44). Through ROC analysis in different cohorts, we found that our risk model showed better efficacy in predicting prognosis compared to models constructed based on other gene signatures, such as pyroptosis in HCC (45). Remarkably, the performance of the model was consistently stable and even better in the validation cohort of ICGC and our own cohort than in the training cohort, which verifies the robustness of our risk model.

Functional analyses revealed that various metabolic pathways were enriched in the mitophagy subgroups and risk groups. Metabolic reprogramming is a hallmark of tumor growth and progression (46). Mitophagy plays a critical role in the metabolic adaptation of cancer cells so that these cells can survive under stress factors produced in the tumor microenvironment, and these adaptions are closely related to the acquisition of metastatic potential and chemoresistance (47). Therefore, some metabolic regulators or pathways related to mitophagy may serve as new therapeutic targets for cancer. Additionally, metabolic pathway regulation can affect immune cell function and fate, leading to a connection to the immune microenvironment (48). This crosstalk between metabolic reprogramming and the immune microenvironment adds further layers to the search for novel therapeutic strategies, regardless of forthcoming challenges. Combining existing evidence and our results, we hypothesize that a metabolism-mitophagy-immunity network exists in HCC, which needs to be explored and validated in future studies.

Nevertheless, our study has several limitations. First, our study focused on the expression of genes, lacking multi-omics data, such as copy number variations and DNA methylation. Second, the study was conducted retrospectively based on data from a public database rather than using a prospective cohort. Furthermore, HCC cell lines and our own HCC cohort used for validation had limited sample sizes, though the results are still reliable. Finally, the mechanisms underlying mitophagy, metabolism, and immunity in tumors warrant further study.

Conclusion

In summary, we identified two prognostically and clinically relevant mitophagy subtypes in HCC. These two subtypes differed in multiple aspects, including immune characteristics, responses to immunotherapy, and biological functions. We also constructed a mitophagy-related risk model that exhibited stable efficiency and performed better than models based on other signatures. The expression of these model genes was subsequently validated using laboratory results. These findings suggest mitophagy as a potential treatment target and shed new light on the strategy of immunotherapy in HCC.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement

The studies involving human participants were reviewed and approved by ethics committee of PUMCH and CAMS (Chinese Academy of Medical Sciences) & PUMC (Peking Union Medical College). The patients/participants provided their written informed consent to participate in this study.

Author contributions

LS designed the study. YW, BP, and CN performed the statistical analysis and wrote the first manuscript. YW and SH performed laboratory experiment. YM and HY collected the data and provided the specimen. LS revised the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by National Natural Science Foundation of China (81602456, 81930053).

Acknowledgments

We are sincerely acknowledged the contributions from the TCGA project and the ICGC project.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.966167/full#supplementary-material

References

1. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer statistics, 2022. CA Cancer J Clin (2022) 72(1):7–33. doi: 10.3322/caac.21708

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Geh D, Manas DM, Reeves HL. Hepatocellular carcinoma in non-alcoholic fatty liver disease-a review of an emerging challenge facing clinicians. Hepatobiliary Surg Nutr (2021) 10(1):59–75. doi: 10.21037/hbsn.2019.08.08

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Kim DW, Talati C, Kim R. Hepatocellular carcinoma (HCC): Beyond sorafenib-chemotherapy. J Gastrointest Oncol (2017) 8(2):256–65. doi: 10.21037/jgo.2016.09.07

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Couri T, Pillai A. Goals and targets for personalized therapy for HCC. Hepatol Int (2019) 13(2):125–37. doi: 10.1007/s12072-018-9919-1

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Sangro B, Sarobe P, Hervás-Stubbs S, Melero I. Advances in immunotherapy for hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol (2021) 18(8):525–43. doi: 10.1038/s41575-021-00438-0

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Onishi M, Yamano K, Sato M, Matsuda N, Okamoto K. Molecular mechanisms and physiological functions of mitophagy. EMBO J (2021) 40(3):e104705. doi: 10.15252/embj.2020104705

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Youle RJ, Narendra DP. Mechanisms of mitophagy. Nat Rev Mol Cell Biol (2011) 12(1):9–14. doi: 10.1038/nrm3028

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Srinivas US, Tan BWQ, Vellayappan BA, Jeyasekharan AD. ROS and the DNA damage response in cancer. Redox Biol (2019) 25:101084. doi: 10.1016/j.redox.2018.101084

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Xu HM, Hu F. The role of autophagy and mitophagy in cancers. Arch Physiol Biochem (2019) 128(2):281–9. doi: 10.1080/13813455.2019.1675714

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Eiyama A, Okamoto K. PINK1/Parkin-mediated mitophagy in mammalian cells. Curr Opin Cell Biol (2015) 33:95–101. doi: 10.1016/j.ceb.2015.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Simon AG, Tolkach Y, Esser LK, Ellinger J, Stöhr C, Ritter M, et al. Mitophagy-associated genes PINK1 and PARK2 are independent prognostic markers of survival in papillary renal cell carcinoma and associated with aggressive tumor behavior. Sci Rep (2020) 10(1):18857. doi: 10.1038/s41598-020-75258-4

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Zhang MY, Huo C, Liu JY, Shi ZE, Zhang WD, Qu JJ, et al. Identification of a five autophagy subtype-related gene expression pattern for improving the prognosis of lung adenocarcinoma. Front Cell Dev Biol (2021) 9:756911. doi: 10.3389/fcell.2021.756911

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Yoon JY, Brezden-Masley C, Streutker CJ. Autophagic heterogeneity in gastric adenocarcinoma. Front Oncol (2021) 11:555614. doi: 10.3389/fonc.2021.555614

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Johnston MP, Khakoo SI. Immunotherapy for hepatocellular carcinoma: Current and future. World J Gastroenterol (2019) 25(24):2977–89. doi: 10.3748/wjg.v25.i24.2977

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res (2015) 43(7):e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res (2003) 13(11):2498–504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics (2010) 26(12):1572–3. doi: 10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf (2013) 14:7. doi: 10.1186/1471-2105-14-7

CrossRef Full Text | Google Scholar

19. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med (2018) 24(10):1550–8. doi: 10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Geeleher P, Cox N, Huang RS. pRRophetic: an r package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One (2014) 9(9):e107468. doi: 10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Software (2010) 33(1):1–22.

Google Scholar

22. Blanche P, Dartigues JF, Jacqmin-Gadda H. Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Stat Med (2013) 32(30):5381–97. doi: 10.1002/sim.5958

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Yu G, Wang LG, Han Y, He QY. clusterProfiler: An r package for comparing biological themes among gene clusters. Omics (2012) 16(5):284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Yapryntseva MA, Zhivotovsky B, Gogvadze V. Induction and detection of mitophagy. Methods Mol Biol (2022) 2445:227–39. doi: 10.1007/978-1-0716-2071-7_14

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Mauro-Lizcano M, Esteban-Martínez L, Seco E, Serrano-Puebla A, Garcia-Ledo L, Figueiredo-Pereira C, et al. New method to assess mitophagy flux by flow cytometry. Autophagy (2015) 11(5):833–43. doi: 10.1080/15548627.2015.1034403

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Shen Z, Zhou H, Li A, Wu T, Ji X, Guo L, et al. Metformin inhibits hepatocellular carcinoma development by inducing apoptosis and pyroptosis through regulating FOXO3. Aging (Albany NY) (2021) 13(18):22120–33. doi: 10.18632/aging.203464

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Wang S, Shi H, Liu T, Li M, Zhou S, Qiu X, et al. Mutation profile and its correlation with clinicopathology in Chinese hepatocellular carcinoma patients. Hepatobiliary Surg Nutr (2021) 10(2):172–9. doi: 10.21037/hbsn.2019.09.17

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Lan XB, Papatheodoridis G, Teng YX, Zhong JH. The upward trend in the immunotherapy utilization for hepatobiliary cancers. Hepatobiliary Surg Nutr (2021) 10(5):692–5. doi: 10.21037/hbsn-21-342

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Craig AJ, von Felden J, Garcia-Lezana T, Sarcognato S, Villanueva A. Tumour evolution in hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol (2020) 17(3):139–52. doi: 10.1038/s41575-019-0229-4

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Poole LP, Macleod KF. Mitophagy in tumorigenesis and metastasis. Cell Mol Life Sci (2021) 78(8):3817–51. doi: 10.1007/s00018-021-03774-1

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Veeriah S, Taylor BS, Meng S, Fang F, Yilmaz E, Vivanco I, et al. Somatic mutations of the parkinson's disease-associated gene PARK2 in glioblastoma and other human malignancies. Nat Genet (2010) 42(1):77–82. doi: 10.1038/ng.491

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Fujiwara M, Marusawa H, Wang HQ, Iwai A, Ikeuchi K, Imai Y, et al. Parkin as a tumor suppressor gene for hepatocellular carcinoma. Oncogene (2008) 27(46):6002–11. doi: 10.1038/onc.2008.199

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Zheng Y, Huang C, Lu L, Yu K, Zhao J, Chen M, et al. STOML2 potentiates metastasis of hepatocellular carcinoma by promoting PINK1-mediated mitophagy and regulates sensitivity to lenvatinib. J Hematol Oncol (2021) 14(1):16. doi: 10.1186/s13045-020-01029-3

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Wu H, Wang T, Liu Y, Li X, Xu S, Wu C, et al. Mitophagy promotes sorafenib resistance through hypoxia-inducible ATAD3A dependent axis. J Exp Clin Cancer Res (2020) 39(1):274. doi: 10.1186/s13046-020-01768-8

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Hao X, Sun G, Zhang Y, Kong X, Rong D, Song J, et al. Targeting immune cells in the tumor microenvironment of HCC: New opportunities and challenges. Front Cell Dev Biol (2021) 9:775462. doi: 10.3389/fcell.2021.775462

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Deng L, He K, Pan Y, Wang H, Luo Y, Xia Q. The role of tumor-associated macrophages in primary hepatocellular carcinoma and its related targeting therapy. Int J Med Sci (2021) 18(10):2109–16. doi: 10.7150/ijms.56003

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Pinter M, Jain RK, Duda DG. The current landscape of immune checkpoint blockade in hepatocellular carcinoma: A review. JAMA Oncol (2021) 7(1):113–23. doi: 10.1001/jamaoncol.2020.3381

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Lee YH, Tai D, Yip C, Choo SP, Chew V. Combinational immunotherapy for hepatocellular carcinoma: Radiotherapy, immune checkpoint blockade and beyond. Front Immunol (2020) 11:568759. doi: 10.3389/fimmu.2020.568759

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Li X, He S, Ma B. Autophagy and autophagy-related proteins in cancer. Mol Cancer (2020) 19(1):12. doi: 10.1186/s12943-020-1138-4

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Zhang Z, Li TE, Chen M, Xu D, Zhu Y, Hu BY, et al. MFN1-dependent alteration of mitochondrial dynamics drives hepatocellular carcinoma metastasis by glucose metabolic reprogramming. Br J Cancer (2020) 122(2):209–20. doi: 10.1038/s41416-019-0658-4

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Cheng J, Qian D, Ding X, Song T, Cai M, Dan X, et al. High PGAM5 expression induces chemoresistance by enhancing bcl-xL-mediated anti-apoptotic signaling and predicts poor prognosis in hepatocellular carcinoma patients. Cell Death Dis (2018) 9(10):991. doi: 10.1038/s41419-018-1017-8

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Liu L, Sakakibara K, Chen Q, Okamoto K. Receptor-mediated mitophagy in yeast and mammalian systems. Cell Res (2014) 24(7):787–95. doi: 10.1038/cr.2014.75

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Kravic B, Harbauer AB, Romanello V, Simeone L, Vögtle FN, Kaiser T, et al. In mammalian skeletal muscle, phosphorylation of TOMM22 by protein kinase CSNK2/CK2 controls mitophagy. Autophagy (2018) 14(2):311–35. doi: 10.1080/15548627.2017.1403716

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Asmarinah A, Paradowska-Dogan A, Kodariah R, Tanuhardja B, Waliszewski P, Mochtar CA, et al. Expression of the bcl-2 family genes and complexes involved in the mitochondrial transport in prostate cancer cells. Int J Oncol (2014) 45(4):1489–96. doi: 10.3892/ijo.2014.2576

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Zheng S, Xie X, Guo X, Wu Y, Chen G, Chen X, et al. Identification of a pyroptosis-related gene signature for predicting overall survival and response to immunotherapy in hepatocellular carcinoma. Front Genet (2021) 12:789296. doi: 10.3389/fgene.2021.789296

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Sun L, Suo C, Li ST, Zhang H, Gao P. Metabolic reprogramming for cancer cells and their microenvironment: Beyond the warburg effect. Biochim Biophys Acta Rev Cancer (2018) 1870(1):51–66. doi: 10.1016/j.bbcan.2018.06.005

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Ferro F, Servais S, Besson P, Roger S, Dumas JF, Brisson L. Autophagy and mitophagy in cancer metabolic remodelling. Semin Cell Dev Biol (2020) 98:129–38. doi: 10.1016/j.semcdb.2019.05.029

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Pearce EL, Pearce EJ. Metabolic pathways in immune cell activation and quiescence. Immunity (2013) 38(4):633–43. doi: 10.1016/j.immuni.2013.04.005

PubMed Abstract | CrossRef Full Text | Google Scholar