Mitophagy-related gene signature predicts prognosis, immune infiltration and chemotherapy sensitivity in colorectal cancer

Abstract

BACKGROUND

Mitophagy plays essential role in the development and progression of colorectal cancer (CRC). However, the effect of mitophagy-related genes in CRC remains largely unknown.

AIM

To develop a mitophagy-related gene signature to predict the survival, immune infiltration and chemotherapy response of CRC patients.

METHODS

Non-negative matrix factorization was used to cluster CRC patients from Gene Expression Omnibus database (GSE39582, GSE17536, and GSE37892) based on mitophagy-related gene expression. The CIBERSORT method was applied for the evaluation of the relative infiltration levels of immune cell types. The performance signature in predicting chemotherapeutic sensitivity was generated using data from the Genomics of Drug Sensitivity in Cancer database.

RESULTS

Three clusters with different clinicopathological features and prognosis were identified. Higher enrichment of activated B cells and CD4⁺ T cells were observed in cluster III patients with the most favorable prognosis. Next, a risk model based on mitophagy-related genes was developed. Patients in training and validation sets were categorized into low-risk and high-risk subgroups. Low risk patients showed significantly better prognosis, higher enrichment of immune activating cells and greater response to chemotherapy (oxaliplatin, irinotecan, and 5-fluorouracil) compared to high-risk patients. Further experiments identified CXCL3 as novel regulator of cell proliferation and mitophagy.

CONCLUSION

We revealed the biological roles of mitophagy-related genes in the immune infiltration, and its ability to predict patients’ prognosis and response to chemotherapy in CRC. These interesting findings would provide new insight into the therapeutic management of CRC patients.

Key Words: Colorectal cancer; Mitophagy; Tumor microenvironment; Immunotherapy; Prognosis

Core Tip: We for the first time revealed the biological roles of mitophagy-related genes in the immune infiltration, and its prediction of patients’ prognosis and chemotherapy response in colorectal cancer. These novel findings would provide new insight into the therapeutic management of colorectal cancer patients.

INTRODUCTION

Colorectal cancer (CRC) is one of the most common malignant diseases in the world[1]. Though the past decades have witnessed a great improvement in treatment programs, the 5-year survival for patients with CRC after tumor resection remains poor and approximately 50% of patients will develop distant metastasis[2]. Thus, effective prognostic markers and novel potential therapeutic targets are necessary to help individualize therapy for CRC patients in clinical practice.

AS a specific type of autophagy, mitophagy is essential in eliminating damaged or excessive mitochondria through autophagolysosomes[3,4]. The process of mitophagy is frequently activated under oxidative and energetic stress, which occur frequently during cancer initiation and progression[5-7]. Previous studies have reported the correlation between mitophagy dysregulation and colorectal cancer formation and invasion. Yin et al[8] revealed that the mitophagy protein PTEN induced putative kinase 1 suppresses CRC cell proliferation by rewiring metabolism via activating p53 and reducing acetyl-CoA production. In addition, induction of mitophagy was shown to contribute to drug resistance in CRC stem cells[9]. A recent study showed that in intestinal epithelial cells, adaptive immunity could be activated by mitophagy during CRC initiation, confirming the correlation between mitophagy and anti-tumor immunity in colorectal cancer[10]. Therefore, targeting mitophagy in CRC cells represents a promising anticancer strategy. Therefore, we aimed to reveal the underlying association between mitophagy and prognosis, immune infiltration and the response to chemotherapy in CRC patients by analyzing the transcriptomics of mitophagy-associated genes in CRC comprehensively.

In this study, we analyzed three datasets (GSE39582, GSE17536, and GSE37892) from the GEO database to elucidate the pattern of mitophagy in CRC. Three mitophagy-related phenotypes in CRC were identified. A prognostic signature was established to quantify the death risk, immune infiltrating characteristics and therapeutic response. Our findings may provide new insights for treatment strategies in CRC based on the mitophagy-related gene signatures.

MATERIALS AND METHODS

Data preparation

GSE39582, GSE17538, and GSE37892 microarray datasets were retrieved from the GEO database, whose website is https://www.ncbi.nlm.nih.gov/geo/. The three abovementioned datasets were generated using the Affymetrix HG-U133 plus 2.0 platform and the Robust Multichip Average method was used to perform data normalization. When multiple datasets were combined as a whole cohort, the ComBat method was used to remove the batch effects. This method was executed in the SVA R package (version 3.2.1).

Dataset selection

The inclusion criteria for dataset selection were: (1) All datasets were generated using the Affymetrix platform (Affymetrix HG-U133 plus 2.0); (2) sample size > 100; (3) patients were pathologically diagnosed with colorectal cancer; (4) basic clinical information for patients was available for analysis; and (5) patients had complete follow-up information. Based on the inclusion criteria, the CRC datasets GSE39582, GSE17538 and GSE37892 were included in this study. The three datasets were combined together and randomly separated into training sets and validation sets at a ratio of 1:1.

Mitophagy-related genes

the Gene Set Enrichment Analysis (http://www.gsea-msigdb.org/gsea/index.jsp) and h GeneCards (https://www.genecards.org/) were utilized to generate a total of 51 genes related to mitophagy. The gene list is shown in Supplementary Table 1.

Clustering by non-negative matrix factorization

We performed consensus clustering “Nonnegative Matrix Factorization (NMF)” with the nrun parameter being set to 50 for CRC patients classification based on the expression levels of mitophagy-related genes[11,12]. This procedure was performed using the “NMF” package in R (version 0.24.0).

Gene set variation analysis

Based on gene expression data, Gene set variation analysis (GSVA) has been widely used to evaluate the activity of signaling pathways via the unsupervised method. The “GSVA” R package (version 1.42.0) was used to quantify the activity of the signaling pathways in each patient. For GSVA processing, the MSigDB database (https://www.gsea-msigdb.org/gsea/index.jsp) was used to download gene set information.

Immune infiltration analysis

We performed a CIBERSORT deconvolution approach to evaluate the relative abundance of 22 tumor-infiltrating immune cells, which included naïve B cells, naïve T cells, memory B cells, activated CD4⁺ memory T cells, plasma cells, resting CD4⁺ memory T cells, CD8⁺ T cells, regulatory T cells (Tregs), CD4⁺ T cells, activated NK cells, follicular helper T cells, resting NK cells, gamma delta T cells, activated mast cells, activated dendritic cells, monocytes, M2 macrophages, M0 macrophages, resting mast cells, M1 macrophages, resting dendritic cells, neutrophils, and eosinophils in each CRC sample[13].

Construction of a mitophagy-related prognostic gene signature

To develop the mitophagy-related prognostic gene signature, we overlapped differentially expressed genes (DEGs) of each cluster identified by “NMF” function to explore the mitophagy related gene expression pattern. The DEGs were identified by lusing LIMMA package (version 3.26.9). Then univariate Cox regression analysis was adopted to identify those significantly associated with overall survival. Next, the least absolute shrinkage and selection operator method (LASSO) Cox regression model was carried out at 10-fold cross-validation to identify the most meaningful mitophagy related genes[14]. The LASSO Cox model was perfomed by lusing glmnet package (version 2.0.5). Finally, the gene signature was established based on selected 36 target genes to calculate risk_score, which was calculated as follows: Risk_score = Σ (Expi × coefi).

Based on the median value of individual risk score, the patients in the training and validation sets were categorized into low and high-risk groups and a Kaplan-Meier survival analysis was conducted to analyze survival difference. The subsequent survival receiver operator characteristic (ROC) as utilized for the evaluation of the prognostic accuracy of the gene signature.

Drug susceptibility analysis

We trained our prognostic signature in pharmacogenomics database (Genomics of Drug Sensitivity in Cancer, https://www.cancerrxgene.org/) to compare the difference of chemotherapeutic sensitivities between low and high-risk patients. The R package “pRRophetic” (version 0.5) was used to conduct the prediction procedure and a ridge regression model was used to estimate the IC₅₀ value of each chemotherapeutic drug.

Statistical analysis

All statistical analysis was performed based on R (version 4.1.2). A Wilcoxon rank-sum test was carried out for data displaying a skewed distribution. A Student’s t-test was performed for data with a normal distribution to compare the differences between groups. A Kruskal-Wallis test was adopted to compare the differences between the three groups. A P value < 0.05 was considered to be significant statistically.

RESULTS

Identification of mitophagy-related subtypes in CRC

A total of 836 CRC patients from the three GEO datasets were identified in this study. The expression correlation among mitophagy-related genes was shown in Figure 1A and survival analysis found that 18 of them were prognostic for CRC (Supplementary Table 2). To further characterize the expression profiling patterns mediated by mitophagy in CRC, we used the “NMF” method and stratified patients into three clusters based on the optimal clustering number (Figure 1B and C). 435 samples were grouped in cluster I, 172 samples in cluster II, and 229 samples in cluster III. The three clusters could be notably discriminated according to PCA analysis (Figure 1D). Survival analysis showed that patients in cluster III exhibited the best prognosis (Figure 1E; P = 0.0057). Further correlation analysis did not observe significant connection between Cluster III and early stage of CRC, suggesting that this identified molecular cluster may be an independent prognostic factor (Figure 1F).

Figure 1 Grouping of mitophagy-related clusters in colorectal cancer. A: Expression correlation among mitophagy-related genes; B: The cophenetic plot of NMF clustering analyses; C: Heatmap showing three clusters were identified, along with the optimal value for consensus clustering; D: Principal component analysis of three mitophagy-related clusters; E: Kaplan-Meier curves for survival rate of three clusters showing patients in cluster III have the best prognosis (P = 0.0057). The P value was calculated using the log-rank test; F: Alluvial diagram of clusters shows patients in cluster III tend to be early-stage. NMF: Nonnegative Matrix Factorization.

Pathway enrichment and immune infiltration analysis among mitophagy-related clusters

Pathway enrichment analysis based on GSVA showed that the metastasis related pathways, such as epithelial mesenchymal transition, KRAS signaling, myogenesis, coagulation, and angiogenesis were mainly enriched in cluster I and cluster II compared to cluster III (Figure 2A and B). To further characterize the immune infiltration among the three clusters, we employed CIBERSORT[15] analyses and found that activated dendritic cells, memory B cells, and activated memory CD4⁺ T cells (Figure 2C) were predominantly enriched in cluster III.

Figure 2 Characterization of the three clusters. A: Differences in metastasis-related pathways scored by gene set variation analysis (GSVA) between cluster I and cluster III; B: Differences in pathway activities scored by GSVA between cluster II and cluster III; C: Relative abundance of 22 tumor-infiltrating immune cells in each of the three clusters. The statistical differences between the three clusters as determined through the Kruskal-Wallis H test. ^aP < 0.05; ^bP < 0.01; ^cP < 0.001.

Construction and validation of the mitophagy-related gene signature

To investigate the prognostic prediction of mitophagy-related genes, we overlapped the DEGs of three clusters and obtained a total of 3112 genes (Figure 3A). Further GO enrichment indicated that these genes were associated with mitochondrial biological function, including mitochondrial gene expression (Figure 3B). Next, univariate Cox analysis was performed and we found 347 prognostic genes with a P value < 0.001 (Supplementary Table 3). The least absolute shrinkage and selection operator (LASSO) regression model was applied to the training set and the minimize λ method resulted in 36 mitophagy-related DEGs (Figure 3C and D). The risk score of each patient was then calculated for the training set based on the established formula. The coefficients of the selected 36 mitophagy-related DEGs were shown in Table S4. Basing on the median risk score, we separated the patients in the training set into low-risk (n = 209) or high-risk (n = 209) groups. The distribution patterns of risk scores and patients’ survival status were shown in Figure 3E and it was noticed that the risk of death increased gradually with the incremental risk scores. The same analyses were then conducted in the validation set and a similar distribution was found (Figure 3F).

Figure 3 Construction and validation of the prognostic risk_score. A: Venn diagram showing overlapped genes between the three clusters; B: Gene ontology enrichment (GO) analysis of mitophagy-related genes; C and D: Thirty-six candidate genes screened out by least absolute shrinkage and selection operator method analysis with minimal lambda; E and F: Ranked dot and scatter plots showing risk_score distribution, patient survival status, survival time and expression of 36 mitophagy-related genes in training set and validation set.

Predictive capacity of the mitophagy-related gene signature

Based on the median mitophagy risk score, we explored patients’ survival rates between the low- and high-risk groups. We observed that high-risk patients showed a significantly worse survival rate than low risk patients in both training (Figure 4A) and validation sets (Figure 4B). Next, we calculated the area under the curve (AUC) values of survival rates in two cohorts. The results showed that the 5-year AUC values all reached 0.7 in the training (Figure 4C) and validation sets (Figure 4D).

Figure 4 Mitophagy-related gene signatures in training and validation sets. A and B: Kaplan–Meier survival analysis between high-risk and low-risk patients; C and D: Receiver operator characteristic curves and area under the curve to predict the sensitivity and specificity of 3-, 5-, and 7-year survival according to the risk_score. AUC: Area under the curve.

After other clinicopathological parameters were adjusted, the mitophagy-related gene signature remained an independent prognostic factor in the whole cohort (Figure 5A). To further facilitate the clinical use of prognostic signatures, a nomogram was developed based on the multivariate result (Figure 5B). The mitophagy-related gene signature, TNM stage and tumor location were integrated into nomogram. The developed nomogram exhibited outstanding prognostic accuracy and displayed superior performance compared to any other factors alone in predicting survival (Figure 5C and D).

Figure 5 Nomogram based on gene signature. A: Hazard ratios of the mitophagy-related gene signature after adjusting for other clinicopathological variables by multivariate Cox regression analysis; B: Nomogram integrating risk_score and clinical characteristics; C: Receiver operator characteristic curves of the nomogram, gene signature, stage and tumor location at 3 years; D: Receiver operator characteristic curves of the nomogram, gene signature, stage and tumor location at 5 years. AUC: Area under the curve.

External validation of mitophagy-related gene signature in TCGA cohort

To further externally validate the prognostic and predictive value of the establised mitophagy-related gene signature, CRC patients were identified from TCGA database. Consistent with the result in training and internal validation sets, the risk of death increased gradually with the incremental risk scores in TCGA cohort (Supplementary Figure 1A). By applying the prognostic signature in TCGA cohort, it was found that patients were successfully separated into two groups with significantly different overall survival (Supplementary Figure 1B). In addition, survival ROC showed that the developed mitophagy-related gene signature performed good in predicting prognosis (5-year AUC = 0.72, Supplementary Figure 1C).

Characterization of Immune infiltration between high- and low-risk groups

Since mitophagy-related gene signatures exhibited good predictive value in predicting survival, we next evaluated the pathway enrichment and immune infiltration patterns in patients from the high- and low-risk groups. The GSVA analysis showed that several important metastasis-related pathways, including epithelial mesenchymal transition, myogenesis, coagulation and angiogenesis were mainly enriched in high-risk group (Figure 6A). By using CIBERSORT analyses, we observed that memory B cells, activated dendritic cells, activated mast cells, plasma cells, and activated memory CD4⁺ T cells were most notably enriched in low-risk patient group (Figure 6B and C).

Figure 6 Immune characteristics of the high- and low-risk groups. A: Differences in metastasis-related pathways scored by gene set variation analysis; B and C: Comparison of the TIICs between the two risk groups. ^aP < 0.05; ^bP < 0.01; ^cP < 0.001.

Prediction of chemotherapy drug sensitivity

Chemotherapy drugs including oxaliplatin, 5-Fluorouracil and irinotecan are the mainstay treatment for CRC. Considering the distinct long-term survival and signaling pathway enrichment between the risk groups, the GSDC database was used to compare the sensitivity to chemotherapy drugs and several targeted inhibitors between low and high-risk patients. It demonstrated the estimated IC₅₀ of oxaliplatin, 5-Fluorouracil and irinotecan were dramatically higher in the high-risk group. In addition, the IC₅₀ of several other inhibitors including the Wnt inhibitor, KRAS inhibitor, Erlotinib, and Afatinib were notably higher in high-risk group (Figure 7).

Figure 7 Estimated IC₅₀ of different chemotherapy drugs between the high- and low-risk groups. A: Oxaliplatin; B: 5-Fluorouracil; C: Irinotecan; D: Cisplatin; E: Wnt_C59; F: KRAS_Inhibitor; G: Gefitinib; H: Erlotinib; I: Nilotinib; J: Afatinib; K: ABT737; L: Paclitaxel. ^aP < 0.05; ^bP < 0.01; ^cP < 0.001.

The effect of CXCL3 on mitophagy in CRC cells

Paired transcriptome profiling data of cancer and normal tissues were obtained from GSE32323 dataset. We focused on the upregulated mitophagy related gene CXCL3, which showed the highest fold change in cancer tissues (Figure 8A). We then enforced CXCL3 expression in CRC cells (Figure 8B). The result showed that CXCL3 over-expression increased the proliferative ability of cells (Figure 8C) and weakened their sensitivity to oxaliplatin treatment (Figure 8D). We next examined mitochondrial DNA (mtDNA) to observe how ULK1 modulated mitophagy. The decrease in mtDNA was attenuated in the groups with CXCL3 overexpressed under hypoxic conditions (Figure 8E and F).

Figure 8 CXCL3 attenuates mitophagy in colorectal cancer cells. A: Differentially expressed mitophagy related genes included in the signature; B: quantitative reverse transcriptase polymerase chain reaction (qRT-PCR) and Western blot showing enforced expression of CXCL3 in colorectal cancer (CRC) cells; C: CXCL3 overexpression significantly enhanced cell proliferation and cell viability; D: CXCL3 overexpression weakened the sensitivity of CRC cells to oxaliplatin treatment; E and F: The mtDNA of CRC cells with or with CXCL3 overexpression was detected using qRT-PCR. ^aP < 0.05; ^bP < 0.01; ^cP < 0.001.

DISCUSSION

In our study, we integrated and analyzed transcriptomic patterns mediated by mitophagy-related genes. Using the “NMF” method, we grouped patients in the GEO datasets into three clusters with distinct clinical features, among which cluster III showed the most favorable prognosis and the most enriched immune cell infiltration. Here, we demonstrated that cluster III exhibited the enrichment of memory B cells, activated dendritic cells and activated memory CD4⁺ T cells. As previously reported, dendritic cells are a diverse group of specialized antigen-presenting cells with key roles in the initiation and regulation of innate and adaptive immune responses[16,17]. With their potent antigen presenting ability, dendritic cells have been long considered as critical factor in antitumor immunity[18]. CD4⁺ T cells, as a major population of cells in the tumor microenvironment (TME), are required for efficacious anti-tumor immunity[19-21]. CD4⁺ T cell are important in fully equipping cytotoxic T lymphocyte cells to target and eliminate cancer cells and facilitate the improvement in outcomes of all cancer immunotherapy strategies[22-25]. In this part, our research delineated intra-tumor heterogeneity of transcriptomic classifications mediated by mitophagy-related genes, and characterized their TME infiltration characteristics. We believe that the molecular classification mediated by mitophagy-related genes could provide new insight for risk stratification and therapeutic management in patients with CRC.

To further quantify the survival risk of CRC patients and facilitate prognostic predictions, we next identified mitophagy subtype related genes and determined the most useful genes to develop the final prognostic signature. Based on the LASSO regression model, a total of 36 mitophagy-related DEGs were identified. In training and validation sets, patients were successfully divided into two groups with significantly different prognoses. Here, we demonstrated that a high risk score was consistently related to a worse prognosis, and that this model exhibited excellent prognostic prediction abilities. In addition, we developed a nomogram based on mitophagy-related gene signatures to facilitate clinical follow-up and management in colorectal cancer. In both the training and validation sets, we confirmed the strong ability of prediction in our nomogram.

Most of the mitophagy-related genes have been experimentally investigated in colorectal cancer and 18 of them were prognostic for CRC[26-31]. We found that PTEN-induced kinase 1 (PINK1) was one of the most significant risk factors. Previous studies have revealed the restrained role of PINK1 mediated mitophagy in innate immunity[32]. In addition, a pan-cancer analysis showed a correlation between PINK1 expression and immune infiltration, including infiltration of B cells, CD8⁺ T cells, CD4⁺ T cells, macrophages, neutrophils, and dendritic cells, suggesting that PINK1 can work as a biomarker for prognosis and immune responses. However, this finding requires further experimental validation[33].

Though a novel mitophagy-related-gene-based signature has been developed in this study, our study mainly focused on the bioinformatic analyses, which still requires a more robust confirmation. Consequently, more in-depth experimental validation is needed to validate these findings.

CONCLUSION

In this study, we have shown for the first time the biological function of mitophagy-related genes in immune infiltration in CRC, and its correlations with patients’ prognosis and chemotherapy response in CRC. These interesting findings would offer novel insights for therapeutic management of CRC patients.

ARTICLE HIGHLIGHTS

Research background

Mitophagy plays essential role in the initiation and progression of colorectal cancer (CRC). Mitophagy plays essential role in the initiation and progression of CRC.

Research motivation

The effect of mitophagy-related genes in CRC remains largely unknown.

Research objectives

To develop a mitophagy-related gene signature to predict the survival, immune infiltration and chemotherapy response of CRC patients.

Research methods

Non-negative matrix factorization was used to cluster CRC patients from Gene Expression Omnibus database (GSE39582, GSE17536, and GSE37892) based on mitophagy-related gene expression. The CIBERSORT method was applied for the evaluation of the relative infiltration levels of immune cell types. The performance signature in predicting chemotherapeutic sensitivity was evaluated based on the Genomics of Drug Sensitivity in Cancer database.

Research results

Three clusters with different clinicopathological features and prognosis were identified. Higher enrichment of activated B cells and CD4⁺ T cells were observed in cluster III patients with the most favorable prognosis. Next, a gene risk model related to mitophagy was produced and patients in training and validation sets were categorized into low-risk and high-risk groups. Low risk patients showed significantly better prognosis, higher enrichment of immune activating cells and greater response to chemotherapy (oxaliplatin, irinotecan and 5-fluorouracil) compared to high risk patients. Further experiments identified CXCL3 as novel regulator of cell proliferation and mitophagy.

Research conclusions

We revealed the biological roles of mitophagy-related genes in the immune infiltration, and its ability to predict patients’ prognosis and response to chemotherapy in CRC. These new findings would offer meaningful insights for the therapeutic management of CRC patients.

Research perspectives

The developed signature in this study will aid in individualizing treatment and follow-up scheme in CRC patients.