The Landscape of Iron metabolism-related Genes for Overall survival prediction in Patients with Hepatocellular Carcinoma

doi:10.21203/rs.3.rs-264542/v1

PDF

Research

The Landscape of Iron metabolism-related Genes for Overall survival prediction in Patients with Hepatocellular Carcinoma

https://doi.org/10.21203/rs.3.rs-264542/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted 04 Mar, 2021

You are reading this older preprint version

Read the latest preprint version →

Hepatocellular carcinoma (HCC) is the seventh most commonly occurring cancer and the second most common cause of cancer-related death worldwide. Despite improvements in early detection and treatment, the morbidity and mortality remain high because of complex molecular mechanisms and cellular heterogeneity in HCC. Immunotherapy therapies have been identified to be an effective treatment strategy for HCC. However, novel model is still needed to predict the survival and clinical immunotherapy response in HCC. Here we established a gene signature using iron metabolism-related genes form the Cancer Genome Atlas (TCGA), and the survival outcomes were validated from International Cancer Genome Consortium(ICGC). Distinct subtypes (high- and low-risk subtypes) were characterized by different clinical outcomes. The high-risk subtype was featured by better survival outcomes, upregulation of immune checkpoints expression, including programmed death-ligand 1 (PD-L1) expression, cytotoxic T-lymphocyte associated protein 4(CTLA-4) expression, T-cell immunoglobulin mucin 3(TIM-3) expression and T cell Ig and ITIM domain (TIGIT) expression, upregulated cell cycle relevant pathways and better response for immunotherapy. Thus our finding suggested that the novel model may be useful as a biomarker for prognostic predication, immunotherapy and cell cycle inhibitors may be efficacious for high-risk subtype of HCC patients.

Translational Medicine

Cancer Biology

Hepatocellular carcinoma

Iron metabolism

Gene signature

Overall survival

Immune status

Due to the complex molecular mechanisms and cellular heterogeneity in HCC,[1–3] it remains high morbidity and mortality in the world. Despite the improvement of treatment, the morbidity and mortality remain high and it is the most rapidly increasing cancer types leading cause of tumor-related death.[4] In the US, 5-year survival rate is only 18%.[5] Surgery is the main treatment in early and middle period.

Unfortunately, the treatment method for advanced and metastatic cases still has a feeble role because of molecular mechanisms and cellular heterogeneity, and metabolic pathway may be involved in the occurrence and progression of HCC, which mediated by iron, nutrient, oxygen, energy stimulation. It is urgently required to construct prediction model to predict the prognosis for patients with HCC. Although several previous studies have revealed different types of signatures for predicting the prognosis of HCC, effective biomarker for prognosis prediction is still lack. Thus, novel prognostic biomarkers are urgently needed to predict the survival of HCC patients and outline an individualized treatment plan for HCC patients.

Disorder of iron metabolism is associated with cell proliferation and unraveled to be involved in most tumors.[6–8] Ferroptosis, as a hot topic of disorder of iron metabolism, regulates cell death via membrane lipid peroxidation.[9, 10] In recently, increased studies have proved the vital role of iron mechanism on the pathogenesis and inhibition of ferroptosis may increase the anticancer activity in HCC.[11–13] Previous studies revealed iron metabolism-related genes in HCC, including NRF2,[13] MT1G, [12]Rb, [11]CISD1,[14]

Inspired by all these works, we constructed an iron metabolic gene signature by performing training cohort from TCGA dataset. Then, the prognostic value of our 13-gene signature was validated in internal validation cohort, external validation cohort. Finally, functional enrichment analysis was performed to explore the underlying mechanisms and provide efficacious therapeutic schedule for patients in high risk group.

2.1 Data Source and grouping

365 primary HCCs candidates and clinical information(age, gender, BMI, vascular tumor cell, grade, TNM stage) were acquired from The cancer genome atlas-liver hepatocellular carcinoma (TCGA‐LIHC, https://cancergenome.nih.gov/ ) dataset. Tumors samples of 243 HCC patients and clinical datasheet(age, gender, priorMalignancy, TNM stage) were identified from International Cancer Genome Consortium (ICGC, https://icgc.org/ ) dataset. 360 HCC patients with complete survival time from TCGA were used as training cohort. By performing R package Caret, 184 HCC patients from TCGA were identified as internal testing cohort with the ratio of 1:1 in a random manner. 243 patients with complete survival time from ICGC dataset were selected as external testing cohort.

2.2 Identification of iron metabolism-related genes

13 iron metabolism-related gene sets (GO_2_IRON_2_SULFUR_CLUSTER_BINDING.gmt, GO_4_IRON_4_SULFUR_CLUSTER_BINDING.gmt,GO_CELLULAR_IRON_ION_HOMEOSTASIS.gmt, GO_HEME_BIOSYNTHETIC_PROCESS.gmt, GO_IRON_COORDINATION_ENTITY_TRANSPORT.gmt, GO_IRON_ION_BINDING.gmt, GO_IRON_ION_HOMEOSTASIS.gmt, GO_IRON_ION_TRANSPORT.gmt, GO_RESPONSE_TO_IRON_ION.gmt, HALLMARK_HEME_METABOLISM.gmt, HEME_BIOSYNTHETIC_PROCESS.gmt, HEME_METABOLIC_PROCESS.gmt, REACTOME_IRON_UPTAKE_AND_TRANSPORT.gmt) were identified from the GSEA (http://software.broadinstitute.org/gsea/index.jsp). we identified 510 iron metabolism-related genes from the 13 gene sets and selected 482 iron metabolism-related genes from TCGA and ICGC gene expression profile.

2.3 Prognostic signature construction

In the training cohort, DEGs associated with iron metabolism were calculated between patients who survived < 1 year and more than 3 years by using the Limma version 3.36.2 R package, and |Log₂ FC|> 1.5 and adjusted P < 0.05 were considered as selection criterion. Univariate Cox proportional hazards regression analysis was performed to assess the relationship between DEGs and the overall survival (OS) to screen prognostic genes, with adjusted P-value < 0.05 as the significance cutoff. An interaction network for the prognostic DEGs was generated by the STRING database (version 11.0) The LASSO Cox regression was applied to minimize overfitting by performing R package glmnet. A multi-gene marker-based prognostic signature was established to assess risk score for HCC patients, the formula is as follows:

Risk score = βgene1 × Egene1 + βgene2 × Egene2 + βgene3 × Egene3 ... βgeneN × EgeneN. Taking the median risk score as a cutoff value, HCC patients in training cohort were divided into low or high risk group.

2.4 The performance of gene signature

Kaplan–Meier (KM) survival curves was performed to compare the survival rate between low or high risk group, and Time-ROC curve analyses the predictive capacity of the gene signature. Next, Univariate and Multivariate Cox regression analysis were used to assess whether the gene signature could be independent of other clinical parameters, including age, gender, BMI, vascular tumor cell, grade, TNM stage. Furthermore, we applied stratification analysis to assess the predictive capacity of the gene signature in different subgroups stratified clinical variables. Finally, the gene signature was validated in external testing cohort, internal testing cohort.

2.5 Identification of Key Co-expression Modules associated with Risk score Using WGCNA

Gene co-expression network was constructed to explore the modules highly associated with sample traits, 9 was selected as the soft thresholding power to produce a weighted network. The enrolled genes were clustered into 9 modules except the gray module, the modules with high correlation coefficient were considered candidates relevant to traits, which were picked out to perform further pathway enrichment analysis using an online-based web tool “Metascape” (http://metascape.org/).

2.6 Gene set enrichment analysis

To identify potential mechanisms underlying our constructed prognostic gene signature to affect the HCC, GSEA was performed to find enriched terms using TCGA-LIHC dataset, and GSEA was performed based on high risk group and low risk group. The FDR < 0.05 was used to sort the pathways enriched in each phenotype.

2.7 Prediction of immunotherapy response

CIBERSORT computational method was implemented to estimate TIC abundance profile in all tumor samples, and samples with p < 0.05 were selected for the following analysis. The difference analysis of immune cells and the activity of immune-related pathways were calculated by performing single-sample gene set enrichment analysis (ssGSEA) in the "gsva" R package.23 ImmuneCellAI (http://bioinfo.life.hust.edu.cn/web/ImmuCellAI/) was used to forecast the immunotherapy response of high and low risk group.

2.8 Building and validating a predictive nomogram for HCC Survival Prediction

Independent prognostic factors identified by multivariate Cox regression analysis were selected to build a nomogram, The calibration curve plotted to observe the prediction probabilities against the observed rates. Time-ROC analysis was performed to compare the prediction accuracy between the nomogram and single independent prognostic factor.

2.9 Human tissue samples

Twelve HCC and adjacent tissue samples were collected from patients -during surgery at The First Affiliated Hospital of Xiamen University; including.

2.10 Quantitative real-time RT-PCR (qRT-PCR) analyses

Total RNA was extracted from the human tissue specimens using TRIzol reagent (Invitrogen; Thermo Fisher Scientific, Inc.). cDNA was synthesized using ReverTra AceH qPCR RT Kit (TOYOBO) and qRT-PCRs were performed using a LightCycler480II (Roche) with SYBR Green technology (Takara). The relative expression levels of the target mRNA were calculated using the 2–ΔΔCt method (ΔCt = Cttarget gene –Ctinternal control). The primer sequences are shown in TableS1. qRT-PCR was performed for 30 cycles consisting of denaturation (at 94°C for 30 sec), annealing (at 56°C for 30 sec) and elongation (at 72°C for 1 min).

3.1 Study Process and Patients' Information

The entire work of construction process of prediction model for HCC patients was presented in Fig. 1A, and the clinical information from training cohort, internal testing cohort and external testing cohort was summarized in Table 1.

Table 1

Summary of patient demographics and clinical characteristics in training cohort and testing cohort.
Characteristic	Training cohort	Internal testing cohort	External testing cohort
Gender
male	187(67%)	89(63%)	182(75%)
female	89(33%)	55(37%)	61(25%)
Age
<60	129(46%)	59(40%)	45(19%)
>=60 BMI <25 >=25 Vascular tumor cell Yes No Grade 1+2 3+4	147(56%) 144(53%) 132(47%) 94(34%) 182(66%) 161(64%) 165(36%)	85(60%) 69(49%) 75(51%) 44(30%) 100(70%) 82(59%) 62(41%)	198(81%) - - - - - -
Stage
stage I-II	223(74%)	120(75%)	146(60%)
stage III-IV	53(26%)	24(25%)	97(40%)
Vital
living	244(68%)	125(68%)	199(82%)
Dead	126(32%)	59(32%)	44(18%)
Risk group
low risk	185(50%)	96(52%)	123(51%)
high risk	185(50%)	88(48%)	120(49%)

3.2 Establishment of the four-gene-based prognostic gene signature

A total of 37 DEGs were obtained when compared patients who survived < 1 year and more than 3 years using training set in the TCGA database (Fig. 1B). 370 patients with complete survival time from training cohort were included for subsequent univariate Cox proportional hazards regression analysis, 25 genes were significantly associated with OS(Fig. 2C,D). Then the LASSO analysis was performed to narrow the gene range. 13 genes were identified and subsequently used to construct a gene signature. Risk-formula: Risk score = expression of CYP3A5*-0.00144471228429986 + expression of CCNB1*0.00932548817179557 + expression of ABCB6*0.115364054791514 + expression of FLVCR1*0.00981797872412028 + expression of OSBP2*0.0960318647350786 + expression of G6PD*0.00325860881257371 + expression of RAP1GAP*0.00568586514806877 + expression of SLC7A11*0.0217165375535805 + expression of PPAT*0.196619205112227 + expression of CYP2C9*-0.000726941457117346 + expression of PLOD2*0.010247507631253 + expression of IGSF3*0.000979749106599136 + expression of RRM2*0.00176828844166581. The interaction network among these genes was presented in Fig. 2E and the correlation between these genes is presented in Fig. 1F. We then calculated the risk score for each patient in training cohort and median risk score was set as cut-off for risk score, 370 patients in training cohort were divided into high risk (185)and low risk (185)group, 184 patients from internal testing cohort were divided into high risk (86)and low risk (98)group, 243 patients from external testing cohort were divided into high risk (185) and low risk (58)group.

3.3 Diagnostic Values of iron metabolism-related signature

In training cohort, Kaplan–Meier curve, time-dependent ROC and PCA analysis were performed to assess the integrated effects of iron metabolism-related signature on the prognosis(Fig. 2). Patients in high-risk group shown significantly poorer OS than patients in the low-risk group(Fig. 2B, P < 0.001), the expression level of 11 prognostic genes increased with risk score increased, including OSBP2, SLC7A11, IGSF3, ABCB6, FLVCR1, PPAT, RAP1GAP, PLOD2, G6PD, CCNB1 and RRM2. While CYP3A5 and CYP2C9 increased with risk score decreasing(Fig. 2A). The area under the ROC curve (AUC) for 1-year, 2-year, and 3-year OS were 0.828, 0.760, 0.739(Fig. 2C). PCA analysis indicated the patients in different risk groups were distributed in two and three directions(Fig. 2D,E), Therefore, the signature have demonstrated a good performance.

3.4 Further Validation of the iron metabolism-related signature using Independent internal testing cohort and external testing cohort

In order to verify the reliability of the 13-Gene Signature, similar procedures were performed in internal testing cohort(Fig. 3), external testing cohort(Fig. 4). Patients in the high-risk group indicated significantly poorer OS than patients in the low-risk group in internal testing cohort (Fig. 3B)and external testing cohort (Fig. 4B)(all P < 0.001). The AUC for 1-year, 2-year, and 3-year OS were 0.807, 0.817, 0.774 in internal testing cohort(Fig. 3C); 0.734, 0.743, 0.743 in external testing cohort(Fig. 4C). PCA analysis indicated displayed a different distribution pattern of high and low risk according to iron metabolism-related signature in internal testing cohort(Fig. 3D,E) and external testing cohort(Fig. 4D,E).

3.5 Independent prognostic role of the gene signature from the Other Clinicopathological Features

Iron metabolism-related signature could be considered as an independent factor in predicting OS. Univariate and Multivariate Cox regression analysis were conducted to assess independent predictive values. In training cohort, Univariate Cox regression analysis indicated that BMI, TNM stage and signature had a prognostic value(Fig. 2F), Multivariate Cox regression analysis revealed that BMI, TNM stage and signature were independent prognostic factor correlation with worse OS(Fig. 2G). In internal testing cohort, Univariate Cox analysis suggested that BMI, TNM stage and signature were associated with OS(Fig. 3F), Multivariate Cox regression analysis showed that BMI, TNM stage and signature were independent prognostic factor associated with OS(Fig. 3G). In external testing cohort, Univariate Cox analysis suggested that gender, TNM stage and signature had a prognostic value(Fig. 4F), Multivariate Cox regression analysis indicated that gender, TNM stage and signature were independent prognostic factor associated with OS(Fig. 4G). These results confirmed that gene signature is an independent predictor of prognosis in HCC patients.

3.6 Validation of expression level of the four genes and risk score in different subgroups

As showed in Fig. 5A-M, the expression level of ABCB6, CCNB1, FLVCR1, G6PD, IGSF3, OSBP2, PLOD2, PPAT, RAP1GAP, RRM2 and SLC7A11 was significantly higher in high risk group compared with low risk group in training cohort, The expression level of CYP3A5 and CYP2C9 was significantly lower in high risk group compared with high risk group in training cohort. Moreover, the risk score in tumor was obviously higher in than normal liver(Fig. 5N), the AUC was 0.932 (95%CI: 0.919–0.945)(Fig. 5O). In addition, the heatmap showed the expression of the nine genes in high- and low-risk patients in the TCGA dataset. It showed that high- and low-risk group associated with age, sex, survival status, grade and stage(Fig. 6A).

3.7 Subgroup Analysis of OS by clinical variables

We performed survival analysis stratified by TNM stage in training cohort(Fig. 6C,D), internal testing cohort(Fig. 6F,G) and external testing cohort(Fig. 6I,J). Risk score grouped based on 13-gene signature, the log-rank test indicated that HCC patients in high risk group still had obviously worser OS than patients in low risk group in different subgroups in stage I-II and stage III-IV. Furthermore, as the TNM stage and tumor grade increased, the risk score increased in training cohort(Fig. 6B), internal testing cohort (Fig. 6E) and external testing cohort(Fig. 6H).

3.8 Comparing the performance of the gene signature with other gene signatures

In recent years, constructing prognostic gene signature for survival prediction has been noted in many studies. Thus, we further compared the gene signature with other gene signatures, the AUC of a seven gene signature(KIF18B, CEP55, CIT, MCM7, CDC45, EZH2, and MCM5) is 0.70, 0.69, 0.67 at 1, 3 and 5-year,[15] the AUC of a three gene signature(UPB1, SOCS2 and RTN3) is 0.72, 0.70, 0.60 at 1, 3 and 5‐year,[16] the AUC of a four gene signature(PBK, CBX2, CLSPN, and CPEB3) is 0.71 at 1-year, 0.66 at 3 ‐year,[17] the AUC of a three eight signature(DCAF13, FAM163A, GPR18, LRP10, PVRIG, S100A9, SGCB, and TNNI3K) is 0.71, 0.69 at 3 and 5‐yea.[18] Peng‐Fei Chen et al[19] have revealed a four-gene signature(KPNA2, CDC20, SPP1, and TOP2A), however, the study lacks of systematic analysis of the four-gene signature, including Univariate and Multivariate Cox regression analysis, time-ROC analysis. Gao‐Min Liu et al[20] have reported a four-gene signature(ACAT1, GOT2, PTDSS2, and UCK2), the AUC of gene signature is 0.70, 0.71, 0.68 at 1-year, 3-year, 5-year, and the AUC of Nomogram is 0.75, 0.75, 0.68 at 1-year, 3-year, 5-year. Our gene signature revealed that the AUC of gene signature is 0.78, 0.68, 0.67 at 1,3,5-year, and the AUC of our Nomogram is 0.78, 0.75, 0.73 at 1,3,5-year. Thus, it indicates that our model has a high sensitivity and specificity for predict the survival, with strong prediction ability of short-term.

3.9 Functional Annotation for Genes of Interest

We identified 936 up-upregulated DEGs and 162 down-regulated DEGs(corrected P-value < 0.05 and absolute log fold change (FC) > 1), Go enrichment analysis revealed that DEGs were mainly involved in mitotic nuclear division, nuclear division, chromosome segregation, nuclear chromosome segregation(Fig. 7A). KEGG enrichment analysis revealed that DEGs were mainly involved in Cell cycle, ECM − receptor interaction, Glycolysis / Gluconeogenesis and so on(Fig. 7B).

Moreover, we conducted GSEA between high and low risk group compared with the median level of risk score to identify the significant pathways (NOM P-value < 0.05), the genes in high risk group were mainly enriched in cell cycle related pathways, including KEGG_CELL_CYCLE, KEGG_OOCYTE_MEIOSIS, KEGG_PURINE_METABOLISM, KEGG_PYRIMIDINE_METABOLISM, KEGG_RNA_DEGRADATION. As to low risk group, the genes were enriched in metabolic pathways, such as KEGG_DRUG_METABOLISM_CYTOCHROME_P450,KEGG_FATTY_ACID_METABOLISM, KEGG_GLYCINE_SERINE_AND_THREONINE_METABOLISM, KEGG_PRIMARY_BILE_ACID_BIOSYNTHESIS,KEGG_TRYPTOPHAN_METABOLISM(Fig. 7C).

To furtherly investigate the potential mechanisms of Risk score in HCC, WGCNA method was performed to cluster genes that highly correlated with the high risk and low risk group. 9 was selected as the soft thresholding power to produce a weighted network(Fig. 8A). We identified 9 modules (Fig. 8B,C) in the study(excluding a gray module). As showed in Fig. 8D, we found genes in blue module was highly positively associated with Risk score. Go enrichment analysis was performed to explore the potential mechanisms of genes in blue module using online-based web tool “Metascape” (http://metascape.org/), as showed in Fig. 8E and Fig. 8F, genes in the blue module were significantly enriched in metabolic process and cell cycle-related process.

3.10 Predictive immunotherapy response of identified HCC subtypes

We employed CIBERSORT algorithm to identify tumor infiltrating immune subsets proportions in HCC, and 21 kinds of immune cell profiles in HCC samples were constructed and showed in Fig. 9A,B. The difference analysis revealed that almost all types of immune cells (excluding B. cells. naive )were significantly different between high- and low-risk group(Fig. 9C). To further explore the correlation between the risk score and immune status, the enrichment scores of ssGSEA associated with functions or pathways were quantified, Cytolytic_activity, MHC_class_I, Type_I_IFN_Reponse and Type_II_IFN_Reponse were significantly different between the low risk and high risk group(Fig. 9D). Finally, ImmuneCellAI ( http://bioinfo.life.hust.edu.cn/web/ImmuCellAI/) was performed to predict immunotherapy response, Fig. 9E and Fig. 9F suggested that patients in the high-risk group might present with a better response for immunotherapy.

3.11 Difference of immune checkpoints expression between 2 HCC subtypes

We performed a correlation analysis between risk score and immune checkpoints, the result revealed risk score positively correlated with immune checkpoints, including PD-L1(CD274)(Fig. 10A), CTLA-4 (Fig. 10B) and TIM-3 (HAVCR2) (Fig. 10C) and TIGIT(Fig. 10D). Expression of PD-L1(Fig. 10E), CTLA-4(Fig. 10F), TIM-3(HAVCR2)(Fig. 10G) and TIGIT(Fig. 10H) were significantly higher in high-risk group versus low-risk group in TCGA.

3.12 Building and validating a Nomogram to Predict OS

Univariate and Multivariate Cox regression analysis have revealed that TNM stage and riskScore were independent prognostic factor in training cohort and testing cohorts. Thus TNM stage and signature were used to build a nomogram to predict 1-year, 3-year, and 5-year OS in the TCGA cohort(Fig. 11A). The calibration curve analysis confirmed the reliability of the nomogram, which demonstrated that 1-, 3-, and 5-year survival prediction probabilities were closely related to the observed survival probability(Fig. 11B,C,D), and the C-index for the nomogram is 0.746342(0.708523–0.784161), which was more large than Signature(0.726819(0.680894–0.772744)) and Stage(0.664469(0.589597–0.73934))(Table 2). Time- ROC curve analysis was applied to evaluate the prediction value of nomogram, all of AUCs of the nomogram was presented in Fig. 11E in training cohort(AUC of 1-year, 3-year, 5-year), in Fig. 11F in internal testing cohort (AUC of 1-year, 3-year, 5-year is 0.808, 0.813, 0.858), in Fig. 11G in external testing cohort(AUC of 1-year, 3-year, 5-year is 0.844, 0.781, 0.7830) .

Table 2

The C-index of Signature, Stage and Combined model.
Biomarker	C-index	95%CI Lower	95%CI Higher	P-value
Signature	0.726819	0.680894	0.772744	3.67E-22
Stage	0.664469	0.589597	0.73934	1.67E-05
Combined model	0.746342	0.708523	0.784161	4.65E-12

3.13 Expression level of genes of risk score in clinical samples

Finally, we compared expression level of genes in risk score between normal and cancer tissues using gene expression profile in TCGA and ICGC, Six genes were significantly different between normal and cancer tissues in TCGA (Fig. 12A-F) and ICGC(Fig. 12G-L). Mrna expression level of G6PD, CCNB1, FLVCR1, ABCB6, RRM2 were significant higher in cancer tissues than normal tissues. However, Mrna expression level of CYP2C9 was significant lower in cancer tissue than normal tissues. In addition, The six genes were significant associated with overall survival (OS), disease-free survival (DFS), disease specific survival (DSS) and progression free survival (PFS) by performing Kaplan-Meier(https://kmplot.com/analysis/). Furthermore, we performed qRT-PCR to detect the expression level in HCC tissues and adjacent samples(Fig. 12M-R). It revealed the same result, Mrna expression level of G6PD, CCNB1, FLVCR1, ABCB6, RRM2 were significant higher in cancer tissue than normal tissue, and Mrna expression level of CYP2C9 was significant lower in cancer tissue than normal tissues.

In the study, we downloaded RNA sequencing, clinical information from TCGA and ICGC dataset. Next, we identified DEGs between patients with > 3 year OS and < 1 year OS by integrated bioinformatics analysis of iron metabolic genes in TCGA dataset. Then, we constructed and validated an iron metabolic signature in TCGA and ICGC dataset, which showed a strong prognostic performance for predicting the survival of HCC patients. In addition, the gene signature was an independent prognostic factor by conducting Multivariate Cox regression analysis. Meantime, WGCNA analysis, GSEA analysis, Go and KEGG enrichment analysis were performed to reveal molecular mechanisms, and cell cycle therapy which may be benefit to high-risk patients. And immunotherapy response prediction was analyzed and revealed immunotherapy may be significant for high risk group. Furthermore, two independent prognostic factors were combined to build Nomogram for predicting HCC prognosis. Finally, we analyzed the expression of the six genes using qRT-PCR in twelve HCC samples.

Iron plays a crucial role in biological processes in mammalian cells by functioning as iron-and haem-containing enzymes, including the cellular stress response, mitochondrial metabolism, chromatin remodeling, DNA repair, mitochondrial metabolism, DNA replication.

Meantime, iron may promote tumour initiation via free radicals. Recently, proteins involved in iron metabolism are reported to contribute to cancer development, including iron efflux pumps, oxidases and reductases, iron regulators as well as siderophore-binding proteins.2

Ferroptosis is iron-dependent form of regulating cell death, which present different characteristics with apoptosis, necrosis.[9] Reactive oxygen species (ROS) may appear in mitochondria during iron-involving oxidative phosphorylation, and excessive ROS could induce oxidative stress response, which could damage large molecular substances such as nucleic acids, lipids and proteins.[21] Recently, an increasing number of studies has indicated that various kinds of diseases are associated with ferroptosis, including periventricular leukomalacia, Huntington's disease and renal functional damage.[22–24] In addition, many kinds of cancers have revealed the close relationships with ferroptosis, including liver cancer, osteosarcoma, prostate adenocarcinoma, cervical carcinoma and diffuse large B‐cell lymphoma.[25, 26]Ferroptosis is also confirmed to be associated with anti-tumor efficacy of chemotherapeutic drugs and immunotherapy. For advanced liver cancer, Sorafenib, a multi‐kinase inhibitor, is the first choice for treatment for 40% of newly diagnosed HCC patients.[27] which could induce and activate ferroptosis to play an anti-cancer role in HCC patients.

The prognostic model contained 13 genes (ABCB6, CCNB1, FLVCR1, G6PD, IGSF3, OSBP2, PLOD2, PPAT, RAP1GAP, RRM2, SLC7A11, CYP3A5 and CYP2C9), Six genes were significantly different between normal and cancer tissues by analyzing TCGA and ICGC dataset, qRT-PCR demonstrated the same trend, including G6PD, CCNB1, FLVCR1, ABCB6, RRM2 and CYP2C9. And survival analysis revealed the six genes were significantly associated OS,DFS, DSS, PFS using Kaplan-Meier. Thus we focused on the six genes. G6PD is a housekeeping enzyme expressed in most cells. Various studies suggested that G6PD deficiency may be protective against HCC in vitro and in vivo and patients with G6PD deficiency showed a 55% risk reduction for HCC compared with patients with normal G6PD activity.[28] Aberrant activation of G6PD is linked to tumorigenesis and malignancy in rapidly growing cancer cells. Previous studies reported that silencing of PPP enzyme G6PD prevents erastin-induced ferroptosis, which reveals G6PD act a protective role by inducing ferroptosis.[9] Now it’s known, RRM2 is one of the genes with high expression and poor prognosis in malignant tumors. During cell division, RR is essential for DNA synthesis. It encodes enzymes that catalyze the conversion of ribonucleotides to deoxyribonucleic acids. The activity of RRM2 subunit reductase requires two Fe3 + ions to form tyrosil radical on tyr122, and tyrosil radical plays a key role in the enzyme activity. [29]Animal experiments revealed that iron deficiency decreased the expression of RRM2 and reduced the protein abundance of RRM2.[30]And RRM2 may inhibit ferroptosis.[31] The expression level of CCNB1 was abnormally elevated in HCC and high expression is associated with better prognosis, as an oncogene, CCNB could inhibit the apoptosis of HCC cells and promote cell invasion.[32] FLVCR1 is a cell surface heme exporter, essential for erythropoiesis and systemic iron homeostasis.[33] The defect of FLVCR1 gene may lead to porphyrin gene overload, excess porphyrins could cause oxidative damage to the cell membrane components and eventually produce cell death.[34]Francesca Vinchi et al investigated the role of Flvcr1a in liver function in mice, the result revealed that FLVCR1 regulates heme synthesis and degradation as well as cytochromes P450 expression and activity.[35] Further study demonstrated that FLVCR1 mediates heme export from macrophages that ingest senescent red cells and regulates hepatic iron. Iron metabolism is disordered in tumors compared with normal tissues and the iron intake gene, FLVCR1, were the most relevant iron metabolism genes in Hepatocellular Carcinoma.[36] ABCB6 encode a mitochondrial transporter. It localizes to both the outer mitochondrial membrane and the plasma membrane. Expression of ABCB6 has been associated with iron homeostasis, [37] and ABCB6 may play a role in heme metabolism. Previous study has found ABCB6 can rescue the rapid growth phenotype of ATM1 deleted yeast cells, which means that ABCB6 acts as a backup system.[38] In HCC, prognostic analysis has indicated that high expression of ABCB6 was associated with worse prognosis. Because ABCB6 may inhibit the ferroptosis of liver cancer cell, and ferroptosis is considered to be a potential anticancer therapy.[39] Increased expression of cytochrome P450 CYP2C9, together with elevated level of its products epoxyeicosatrienoic acids, is associated with aggressiveness in cancer, previous study revealed that high activity CYP2C9 genotype is associated with increased risk of colorectal cancer, [40] while a low activity CYP2C9 genotype (CYP2C9*2) is associated with increased risk of colorectal adenoma [41] and lung cancer[42] and low expression is associated with poor prognosis in HCC.[43] which is consistent with our study.

As a notable immunotherapy, immune checkpoint blockade have emerged as powerful tools in the management of multiple cancer types, and immune checkpoint drugs have already been FDA approved in renal cell carcinoma, melanoma, non-small cell lung cancer, Hodgkin lymphoma and urothelial bladder cancer, including pembrolizumab, ipilimumab and nivolumab.[44] Several trials investigating the efficacy of immune checkpoint blockades in HCC are in progress (NCT02837029, NCT02947165, NCT02572687, NCT02740985, NCT01853618, NCT02668770, NCT02239900, NCT02423343, NCT01658878). The treatment of anti-PD-1 antibody nivolumab in HCC demonstrated a significant and stable response.[45] Thus it is promising for HCC by applying immunotherapy. In our study, we found risk score has a strongly positive correlation with immune checkpoints, and the high-risk subtype had significantly higher immune checkpoints expression, including PD-L1,PD-L2, CTLA-4, and TIM-3. Moreover, ImmuneCellAI indicated that high-risk group might present with a better response for immunotherapy. Thus high-risk patients may have a better response for immunotherapy.

Cell cycle relevant pathways were significantly upregulated in high risk group, which suggested cell cycle inhibitors may be more suitable for high-risk patients. CDK4 and CDK6 inhibitors have been reported to enhance the response of immune checkpoint blockade therapies in vitro.[46, 47] and CDK4 and CDK6 inhibitors could T cell exclusion signature status to get a better response in vitro.[47] Thus combination of cell cycle inhibitors and immunotherapy may be synergetic to enhance antitumous effect for patients in high risk group for HCC.

Several limitations should be acknowledged in our study. Firstly, a multicenter and prospective cohorts are needed to validate our results. Secondly, effect of combination of cell cycle inhibitors and immunotherapy is also needed to elucidate in vivo and in vitro experiment.

Our finding suggested that the novel model may be useful as a biomarker for prognostic predication, immunotherapy and cell cycle inhibitors may be efficacious for high-risk subtype of HCC patients.

HCC Hepatocellular carcinoma

TCGA the Cancer Genome Atlas

ICGC International Cancer Genome Consortium

PD-L1 programmed death-ligand 1

CTLA-4 cytotoxic T-lymphocyte associated protein 4

TIM-3 T-cell immunoglobulin mucin 3

TIGIT T cell Ig and ITIM domain

OS overall survival

KM Kaplan–Meier

ssGSEA gene set enrichment analysis

AUC area under the ROC curve

DFS disease-free survival

DSS disease specific survival

PFS progression free survival

Acknowledgements

The authors thank all individuals who participated in this study and donated samples.

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.

Author Contributions

Zhipeng Zhu came up with the design and conception. Zhipeng Zhu, Huang Cao, Heqing Huang, Hongliang Zhan, Zhengsheng Liu, Kaihong Lu, Ting Zhang, Chaohao Miao performed data analysis and visualization. The original writing of the draft was by Zhipeng Zhu. The original writing of the editing was by Zhun Wu.

Funding

No funding

Availability of data and materials

The original contributions presented in the study are publicly available. This data can be found here: (https://portal.gdc.cancer.gov/）and (https://dcc.icgc.org/).

Ethics approval and consent to participate

The human tissue samples were approved by the ethics committee of the First Affiliated Hospital of Xiamen University. Access to the de-identified linked dataset was obtained from the TCGA and GEO databases in accordance with the database policy. For analyses of de-identified data from the TCGA and GEO databases.

Consent for publication

Not applicable.

Hoshida Y, Nijman SM, Kobayashi M, Chan JA, Brunet JP, Chiang DY, Villanueva A, Newell P, Ikeda K, Hashimoto M et al: Integrative transcriptome analysis reveals common molecular subclasses of human hepatocellular carcinoma. Cancer Res 2009, 69(18):7385-7392.
Nault JC, Villanueva A: Intratumor molecular and phenotypic diversity in hepatocellular carcinoma. Clin Cancer Res 2015, 21(8):1786-1788.
Torrecilla S, Sia D, Harrington AN, Zhang Z, Cabellos L, Cornella H, Moeini A, Camprecios G, Leow WQ, Fiel MI et al: Trunk mutational events present minimal intra- and inter-tumoral heterogeneity in hepatocellular carcinoma. J Hepatol 2017, 67(6):1222-1231.
El-Serag HB, Rudolph KL: Hepatocellular carcinoma: epidemiology and molecular carcinogenesis. Gastroenterology 2007, 132(7):2557-2576.
Cronin KA, Lake AJ, Scott S, Sherman RL, Noone AM, Howlader N, Henley SJ, Anderson RN, Firth AU, Ma J et al: Annual Report to the Nation on the Status of Cancer, part I: National cancer statistics. Cancer 2018, 124(13):2785-2800.
Huang X: Iron overload and its association with cancer risk in humans: evidence for iron as a carcinogenic metal. Mutat Res 2003, 533(1-2):153-171.
Manz DH, Blanchette NL, Paul BT, Torti FM, Torti SV: Iron and cancer: recent insights. Ann N Y Acad Sci 2016, 1368(1):149-161.
Zhang S, Chang W, Wu H, Wang YH, Gong YW, Zhao YL, Liu SH, Wang HZ, Svatek RS, Rodriguez R et al: Pan-cancer analysis of iron metabolic landscape across the Cancer Genome Atlas. J Cell Physiol 2020, 235(2):1013-1024.
Dixon SJ, Lemberg KM, Lamprecht MR, Skouta R, Zaitsev EM, Gleason CE, Patel DN, Bauer AJ, Cantley AM, Yang WS et al: Ferroptosis: an iron-dependent form of nonapoptotic cell death. Cell 2012, 149(5):1060-1072.
Stockwell BR, Friedmann Angeli JP, Bayir H, Bush AI, Conrad M, Dixon SJ, Fulda S, Gascon S, Hatzios SK, Kagan VE et al: Ferroptosis: A Regulated Cell Death Nexus Linking Metabolism, Redox Biology, and Disease. Cell 2017, 171(2):273-285.
Louandre C, Marcq I, Bouhlal H, Lachaier E, Godin C, Saidak Z, Francois C, Chatelain D, Debuysscher V, Barbare JC et al: The retinoblastoma (Rb) protein regulates ferroptosis induced by sorafenib in human hepatocellular carcinoma cells. Cancer Lett 2015, 356(2 Pt B):971-977.
Sun X, Niu X, Chen R, He W, Chen D, Kang R, Tang D: Metallothionein-1G facilitates sorafenib resistance through inhibition of ferroptosis. Hepatology 2016, 64(2):488-500.
Sun X, Ou Z, Chen R, Niu X, Chen D, Kang R, Tang D: Activation of the p62-Keap1-NRF2 pathway protects against ferroptosis in hepatocellular carcinoma cells. Hepatology 2016, 63(1):173-184.
Yuan H, Li X, Zhang X, Kang R, Tang D: CISD1 inhibits ferroptosis by protection against mitochondrial lipid peroxidation. Biochem Biophys Res Commun 2016, 478(2):838-844.
Xiang XH, Yang L, Zhang X, Ma XH, Miao RC, Gu JX, Fu YN, Yao Q, Zhang JY, Liu C et al: Seven-senescence-associated gene signature predicts overall survival for Asian patients with hepatocellular carcinoma. World J Gastroenterol 2019, 25(14):1715-1728.
Li B, Feng W, Luo O, Xu T, Cao Y, Wu H, Yu D, Ding Y: Development and Validation of a Three-gene Prognostic Signature for Patients with Hepatocellular Carcinoma. Sci Rep 2017, 7(1):5517.
Yan Y, Lu Y, Mao K, Zhang M, Liu H, Zhou Q, Lin J, Zhang J, Wang J, Xiao Z: Identification and validation of a prognostic four-genes signature for hepatocellular carcinoma: integrated ceRNA network analysis. Hepatol Int 2019, 13(5):618-630.
Qiao GJ, Chen L, Wu JC, Li ZR: Identification of an eight-gene signature for survival prediction for patients with hepatocellular carcinoma based on integrated bioinformatics analysis. PeerJ 2019, 7:e6548.
Chen PF, Li QH, Zeng LR, Yang XY, Peng PL, He JH, Fan B: A 4-gene prognostic signature predicting survival in hepatocellular carcinoma. J Cell Biochem 2019, 120(6):9117-9124.
Liu GM, Xie WX, Zhang CY, Xu JW: Identification of a four-gene metabolic signature predicting overall survival for hepatocellular carcinoma. J Cell Physiol 2020, 235(2):1624-1636.
Ray PD, Huang BW, Tsuji Y: Reactive oxygen species (ROS) homeostasis and redox regulation in cellular signaling. Cell Signal 2012, 24(5):981-990.
Friedmann Angeli JP, Schneider M, Proneth B, Tyurina YY, Tyurin VA, Hammond VJ, Herbach N, Aichler M, Walch A, Eggenhofer E et al: Inactivation of the ferroptosis regulator Gpx4 triggers acute renal failure in mice. Nat Cell Biol 2014, 16(12):1180-1191.
Linkermann A, Skouta R, Himmerkus N, Mulay SR, Dewitz C, De Zen F, Prokai A, Zuchtriegel G, Krombach F, Welz PS et al: Synchronized renal tubular cell death involves ferroptosis. Proc Natl Acad Sci U S A 2014, 111(47):16836-16841.
Skouta R, Dixon SJ, Wang J, Dunn DE, Orman M, Shimada K, Rosenberg PA, Lo DC, Weinberg JM, Linkermann A et al: Ferrostatins inhibit oxidative lipid damage and cell death in diverse disease models. J Am Chem Soc 2014, 136(12):4551-4556.
Sun X, Ou Z, Xie M, Kang R, Fan Y, Niu X, Wang H, Cao L, Tang D: HSPB1 as a novel regulator of ferroptotic cancer cell death. Oncogene 2015, 34(45):5617-5625.
Yang WS, SriRamaratnam R, Welsch ME, Shimada K, Skouta R, Viswanathan VS, Cheah JH, Clemons PA, Shamji AF, Clish CB et al: Regulation of ferroptotic cancer cell death by GPX4. Cell 2014, 156(1-2):317-331.
Louandre C, Ezzoukhry Z, Godin C, Barbare JC, Maziere JC, Chauffert B, Galmiche A: Iron-dependent cell death of hepatocellular carcinoma cells exposed to sorafenib. Int J Cancer 2013, 133(7):1732-1742.
Dore MP, Vidili G, Marras G, Assy S, Pes GM: Inverse Association between Glucose6Phosphate Dehydrogenase Deficiency and Hepatocellular Carcinoma. Asian Pac J Cancer Prev 2018, 19(4):1069-1073.
Asperti M, Cantamessa L, Ghidinelli S, Gryzik M, Denardo A, Giacomini A, Longhi G, Fanzani A, Arosio P, Poli M: The Antitumor Didox Acts as an Iron Chelator in Hepatocellular Carcinoma Cells. Pharmaceuticals (Basel) 2019, 12(3).
Key J, Sen NE, Arsovic A, Kramer S, Hulse R, Khan NN, Meierhofer D, Gispert S, Koepf G, Auburger G: Systematic Surveys of Iron Homeostasis Mechanisms Reveal Ferritin Superfamily and Nucleotide Surveillance Regulation to be Modified by PINK1 Absence. Cells 2020, 9(10).
Zhang X, Du L, Qiao Y, Zhang X, Zheng W, Wu Q, Chen Y, Zhu G, Liu Y, Bian Z et al: Ferroptosis is governed by differential regulation of transcription in liver cancer. Redox Biol 2019, 24:101211.
Jin J, Xu H, Li W, Xu X, Liu H, Wei F: LINC00346 Acts as a Competing Endogenous RNA Regulating Development of Hepatocellular Carcinoma via Modulating CDK1/CCNB1 Axis. Front Bioeng Biotechnol 2020, 8:54.
Khan AA, Quigley JG: Heme and FLVCR-related transporter families SLC48 and SLC49. Mol Aspects Med 2013, 34(2-3):669-682.
Krishnamurthy P, Xie T, Schuetz JD: The role of transporters in cellular heme and porphyrin homeostasis. Pharmacol Ther 2007, 114(3):345-358.
Vinchi F, Ingoglia G, Chiabrando D, Mercurio S, Turco E, Silengo L, Altruda F, Tolosano E: Heme exporter FLVCR1a regulates heme synthesis and degradation and controls activity of cytochromes P450. Gastroenterology 2014, 146(5):1325-1338.
Shen Y, Li X, Zhao B, Xue Y, Wang S, Chen X, Yang J, Lv H, Shang P: Iron metabolism gene expression and prognostic features of hepatocellular carcinoma. J Cell Biochem 2018, 119(11):9178-9204.
Paterson JK, Shukla S, Black CM, Tachiwada T, Garfield S, Wincovitch S, Ernst DN, Agadir A, Li X, Ambudkar SV et al: Human ABCB6 localizes to both the outer mitochondrial membrane and the plasma membrane. Biochemistry 2007, 46(33):9443-9452.
Zutz A, Gompf S, Schagger H, Tampe R: Mitochondrial ABC proteins in health and disease. Biochim Biophys Acta 2009, 1787(6):681-690.
Zhang J, Zhang X, Li J, Song Z: Systematic analysis of the ABC transporter family in hepatocellular carcinoma reveals the importance of ABCB6 in regulating ferroptosis. Life Sci 2020, 257:118131.
Martinez C, Garcia-Martin E, Ladero JM, Sastre J, Garcia-Gamito F, Diaz-Rubio M, Agundez JA: Association of CYP2C9 genotypes leading to high enzyme activity and colorectal cancer risk. Carcinogenesis 2001, 22(8):1323-1326.
Barry EL, Poole EM, Baron JA, Makar KW, Mott LA, Sandler RS, Ahnen DJ, Bresalier RS, McKeown-Eyssen GE, Ulrich CM: CYP2C9 variants increase risk of colorectal adenoma recurrence and modify associations with smoking but not aspirin treatment. Cancer Causes Control 2013, 24(1):47-54.
Elfaki I, Mir R, Almutairi FM, Duhier FMA: Cytochrome P450: Polymorphisms and Roles in Cancer, Diabetes and Atherosclerosis. Asian Pac J Cancer Prev 2018, 19(8):2057-2070.
Wang X, Yu T, Liao X, Yang C, Han C, Zhu G, Huang K, Yu L, Qin W, Su H et al: The prognostic value of CYP2C subfamily genes in hepatocellular carcinoma. Cancer Med 2018, 7(4):966-980.
Shah UA, Nandikolla AG, Rajdev L: Immunotherapeutic Approaches to Biliary Cancer. Curr Treat Options Oncol 2017, 18(7):44.
El-Khoueiry AB, Sangro B, Yau T, Crocenzi TS, Kudo M, Hsu C, Kim TY, Choo SP, Trojan J, Welling THR et al: Nivolumab in patients with advanced hepatocellular carcinoma (CheckMate 040): an open-label, non-comparative, phase 1/2 dose escalation and expansion trial. Lancet 2017, 389(10088):2492-2502.
Goel S, DeCristo MJ, Watt AC, BrinJones H, Sceneay J, Li BB, Khan N, Ubellacker JM, Xie S, Metzger-Filho O et al: CDK4/6 inhibition triggers anti-tumour immunity. Nature 2017, 548(7668):471-475.
Jerby-Arnon L, Shah P, Cuoco MS, Rodman C, Su MJ, Melms JC, Leeson R, Kanodia A, Mei S, Lin JR et al: A Cancer Cell Program Promotes T Cell Exclusion and Resistance to Checkpoint Blockade. Cell 2018, 175(4):984-997 e924.

PDF

Version 1

posted 04 Mar, 2021

You are reading this older preprint version

Read the latest preprint version →

The Landscape of Iron metabolism-related Genes for Overall survival prediction in Patients with Hepatocellular Carcinoma

Status:

Version 1

Abstract

Figures

1. Introduction

2. Materials And Methods

2.1 Data Source and grouping

2.2 Identification of iron metabolism-related genes

2.3 Prognostic signature construction

2.4 The performance of gene signature

2.5 Identification of Key Co-expression Modules associated with Risk score Using WGCNA

2.6 Gene set enrichment analysis

2.7 Prediction of immunotherapy response

2.8 Building and validating a predictive nomogram for HCC Survival Prediction

2.9 Human tissue samples

2.10 Quantitative real-time RT-PCR (qRT-PCR) analyses

3. Result

3.1 Study Process and Patients' Information

3.2 Establishment of the four-gene-based prognostic gene signature

3.3 Diagnostic Values of iron metabolism-related signature

3.4 Further Validation of the iron metabolism-related signature using Independent internal testing cohort and external testing cohort

3.5 Independent prognostic role of the gene signature from the Other Clinicopathological Features

3.6 Validation of expression level of the four genes and risk score in different subgroups

3.7 Subgroup Analysis of OS by clinical variables

3.8 Comparing the performance of the gene signature with other gene signatures

3.9 Functional Annotation for Genes of Interest

3.10 Predictive immunotherapy response of identified HCC subtypes

3.11 Difference of immune checkpoints expression between 2 HCC subtypes

3.12 Building and validating a Nomogram to Predict OS

3.13 Expression level of genes of risk score in clinical samples

4. Discussion

5. Conclusions

Abbreviations

Declarations

Acknowledgements

Conflict of Interest

Author Contributions

Funding

Availability of data and materials

Ethics approval and consent to participate

Consent for publication

References

Status:

Version 1