Identication of Iron Metabolism-Related Genes Signature and Novel Role of FLVCR1 in Hepatocellular Carcinoma

Background: Hepatocellular carcinoma (HCC) is the main and highly malignant histological subtype of liver cancer. We tried to construct a novel signature with iron metabolism-related genes to provide new therapeutic targets and improve the prognosis for HCC patients. Methods: The gene expression data of 70 iron metabolism-related genes and its relevant clinical information were obtained from The Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC) databases. Consensus clustering analysis was performed to determine clusters of HCC patients with different OS. Cox regression and LASSO regression analyses were used to establish a prognostic signature. Receiver operating characteristic (ROC) and Kaplan–Meier analyses were carried out to examine the predicated performance of the signature. Results: Consensus clustering analysis determined two clusters of HCC patients with different OS(p<0.01), TNM stage(p<0.05) and pathological grade(p<0.05). A nine-gene prognostic signature established with iron metabolism-related genes can independently predicate the prognostic of HCC patients. The ROC curves showed a great performance of the signature. In addition, FLVCR1, a hub gene with the highest mutation frequency in our signature, showed the signicantly prognostic value in HCC patients. High FLVCR1 expression was signicantly associated with poor prognosis and aggressive progression in HCC patients. The promoter methylation level of FLVCR1 was lower in HCC samples with aggressive progression status. The FLVCR1 expression was positively correlated with the inltration level of B cell, CD4+ T cell, macrophage, neutrophil and dendritic cell. Conclusion: Our study rst established a signature related to iron metabolism and identied FLVCR1 as a potential therapeutic target. These ndings provided more treatment strategies for HCC patients. processes cell epithelial to mesenchymal transition(25), JAK/STAT and ERK we performed functional analyses of the differentially regulated genes between the two clusters of HCC. the results showed that these genes are signicantly related to cancer-related biological processes and signaling pathways. Combined with these results, it could be determined that iron metabolism-related genes play a critical effect in regulating the malignant process of


Introduction
Hepatocellular carcinoma (HCC) is the most common histological subtype of liver cancer, accounting for > 80% of primary liver cancer (1). It is estimated that HCC is the fourth most common cause of cancerrelated deaths worldwide, and incidence and mortality continue to increase (2,3). Approximately 85% of HCC cases are assessed to occur in developing countries or regions, especially in Eastern Asia and sub-Saharan Africa (4). Chronic hepatitis B virus (HBV) and hepatitis C virus (HCV) infection are the main risk factors for HCC (5). Surgical resection, liver transplantation, tumor ablation, chemotherapy have been proven to be bene cial to the survival of HCC patients (6). However, most patients are diagnosed at an advanced stage, which leads to awfully poor prognosis for HCC patients.
Iron is an indispensable trace element for human health, and it plays a special role in the biological processes of cells, such as DNA synthesis, oxygen transport and cellular respiration (7). In general, iron homeostasis is maintained in healthy cells by balancing iron uptake, utilization, storage and export.
Dysregulation of this balance is related to various diseases, including cancer (8). Numerous studies have indicated that iron has dual properties in the progression of cancer(9, 10). On the one hand, excess iron provides fuel to meet the demands of the rapid growth of tumor cells (11). The aberrant expression of iron metabolism-related genes in cancer typically has the characteristics of upregulation of genes for iron uptake and downregulation of genes for iron out ow (12). Yet, excessive iron can produce reactive oxygen species (ROS) and destroy lipids in cell membranes, thereby triggering ferroptosis (13). Ferroptosis is a novel type of iron-dependent programmed cell death, which seems different from apoptosis, necrosis, and autophagy (14). In recent years, ferroptosis has been regarded as a new strategy of cancer treatment, and a variety of ferroptosis inducers have been applied (15,16).
Since the liver is the chief organ for storage of iron, excess iron is an important risk factor for HCC (17).
Therefore, targeting iron metabolism is one of the promising therapeutic strategies in the treatment of HCC. In this study, we tried to rst construct and validate a reliable prognostic model using iron metabolism-related genes for HCC patients. We suggested that our model can provide new targets and strategies for HCC treatment.

Data Sources
Seventy iron metabolism-related genes were collected from the published literature (18). RNA-seq transcriptome and clinical data sets were downloaded from The Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC)databases. The expression data of FLVCR1 in Wurmbach and Guichard cohorts were obtained from the Oncomine database. Patients with unknown clinicopathological characteristics, survival time and status are excluded.

Consensus clustering analysis
Based on the expression pattern of iron metabolism-related genes, all tumor samples were clustered into several subgroups by using "ConsensusClusterPlus" package of R. The difference of overall survival (OS) of HCC patients in different subgroups were analyzed using the 'survival' package.

Functional analyses
GO (gene ontology) function and KEGG (kyoto encyclopedia of genes and genomes) pathway analyses were performed using clusterpro ler package, and P<0.05 was used as a threshold to screen the major enrichment functions and pathways. Gene set enrichment analyses (GSEA) were conducted to observe the top enriched functions and signal pathways in the Molecular Signatures Database with adj. p < 0.05.

Construction and validation of a prognostic signature
We rst performed univariate analysis on the TCGA cohort to screen for genes associated with overall survival (OS) (p<0.05), and used LASSO regression analysis to further select them. Finally, the genes were identi ed, and their regression coe cient(β) were obtained using multivariate cox regression analysis.
After calculating the risk score of each patient, all patients were divided into two risk groups according to the median risk score. The predicated capability of the prognostic signature was evaluated using receiver operating characteristic (ROC) and Kaplan-Meier analysis in TCGA and ICGC cohorts.

Establishment of a nomogram
Cox regression analyses were conducted to estimate the prognostic value of multi-gene signature and clinicopathologic characteristics. Independent prognostic factors identi ed by multivariate Cox regression analysis were incorporated into the establishment of a prognostic nomogram. A calibration plot was conducted to evaluate predicated performance of the prognostic model.

Cell culture
The human liver cell line L-02 and HCC cell lines SMMC-7721, CSQT2, Huh7 and HepG2 were purchased from the Shanghai Institute of Biosciences Cell Resource Center, Chinese Academy of Sciences. All cells were cultured in Dulbecco's modi ed Eagle's medium (DMEM; Gibco) containing 10% fetal bovine serum, and placed at 37 °C in a cell culture incubator with 5% CO2.

Quantitative real-time polymerase chain reaction (qRT-PCR)
Collected cells were thoroughly ground by adding an appropriate amount of TRIzol (Invitrogen Corporation), and RNA was extracted based on the kit directions. cDNA was nally obtained by a reverse transcription kit (PrimeScript RT reagent kit, TaKaRa Corporation). The speci c primers used were as follows: FLVCR1 (F:5'-GAGCAGAGATGTCACATGGCGAAG-3', R:5'-CTCATGGCTTGGTGCTGTCCTTG-3').
8. Immune in ltration TIMER (https://cistrome.shinyapps.io/timer/) is a website for systematically analyzing immune in ltration across diverse cancer types. In this research, relationships between the in ltration level of tumor immune in ltration cells (TIICs) and FLVCR1 expression or copy number variation were conducted. 9. Statistical analyses SPSS 21.0, Graph Prism7.0 and R software 3.6.1 were used for all statistical analyses. In general, p<0.05 is considered a statistically signi cant difference.

The molecular mechanism of iron metabolism-related genes
To explore the molecular mechanism of 70 iron metabolism-related genes, the gene enrichment analyses were performed. GO analysis demonstrated that iron metabolism related mechanisms such as iron ion homeostasis, transition metal ion homeostasis, cellular iron ion homeostasis, cellular transition metal ion homeostasis, iron ion transport and iron ion transmembrane transport etc. were enriched (Figure1a).
Ferroptosis, mineral absorption, porphyrin and chlorophyll metabolism, ABC transporters and vitamin digestion and absorption were associated with these genes in KEGG analysis (Figure1b).

Association between iron metabolism-related genes and HCC prognosis
Based on the expression pattern of 70 iron metabolism-related genes, we determined the optimal specimen clustering number as two via consensus clustering analysis (Figure2a). The Kaplan-Meier curve revealed that the HCC patients in cluster 2 had a distinctly worse OS than those in cluster 1 (p =0.006) (Figure2b). Then, the relationship analyses of clinicopathological characteristics between the two clusters showed signi cant differences for the stage and grade (p < 0.05) (Figure2c). These ndings indicated that the clustering results were associated with the prognosis and malignancy of HCC.
Then, we further explored the difference between the two clusters. Setting adjust p = 0.05 and |log2FC| =0.5, we screened the differentially expressed genes in cluster 2 compared to cluster 1 (Figure3a).
Subsequently, we annotated the biological functions of these genes using functional enrichment analyses. GO enrichment showed that differential genes were signi cantly enriched in terms associated with cell junction, such as cell-cell junction, collagen-containing extracellular matrix and bicellular tight All patients in TCGA and ICGC cohorts were divided into two risk groups based on the risk score. The survival curves displayed that patients with high risk had shorter OS time than those with low risk (p<0.001). Consistent with the result was found in the ICGC cohort(p<0.001) (Figure5a). The distribution and survival status of patients for the nine-gene prognostic signature were shown in the TCGA and ICGC cohorts (Figure5b). The number of deaths increased as the risk score increased. ROC analysis was performed to assess the predictive capacity of the prognostic signature at 1, 2, and 3 years. The results showed that AUCs for 1-, 2-and 3-year OS reached 0.806, 0.781 and 0.750 in the TCGA cohort and 0.747, 0.772 and 0.813 in the ICGC cohort, respectively (Figure5c).

Independent prognostic and predictive value of the multi-gene prognostic signature
To examine whether the prognostic signature can independently predict the prognosis of HCC patients, the Cox regression analyses were carried out again. After univariate Cox regression analysis, stage and risk score were found to be closely related to OS in HCC patients (Figure6a). After multivariate Cox regression analysis, stage and risk score were identi ed as independent prognostic factors (Figure6c). Consistent results were obtained in the ICGC cohort (Figure6b, d). Then, all independent prognostic factors determined by Cox regression analyses were used to construct a nomogram for predicting the 1-

FLVCR1 gene expression in HCC and its impact on prognosis
PPI network showed that FLVCR1, ABCB6 and ALAS1 was hub genes of the prognostic signature (Figure8a). In addition, we examined the changes in genetic alterations of the nine genes and eventually found that ampli cation was the most common type of mutation (Figure8b). Among them, the FLVCR1 gene had the highest mutation frequency. These results implied that FLVCR1 had a special role in the progression of HCC.
As the results showed, FLVCR1 expression was abnormally upregulated in several cancers, such as BLCA, HCC and STAD (Figure8c). The upregulated expression of FLVCR1 in HCC was veri ed again based on the ICGC, Wurmbach and Guichard cohorts (Figure8d). The RT-qPCR analysis revealed that mRNA expression levels of FLVCR1 were higher in HCC cells than LO2 cells (Figure8e). The relationship analysis between clinicopathological characteristics and FLVCR1 showed signi cant difference for the age(p=0.0015), T(p=0.0058), stage(p<0.001) and grade(p<0.001) (Figure9a). Promoter methylation levels of FLVCR1 were signi cantly lower in primary HCC samples than normal samples(p<0.001). Patients with TP53 mutation, higher tumor grade and stage had higher promoter methylation levels of FLVCR1(Figure9b). In addition, survival curves demonstrated that high FLVCR1 expression was related to poor OS, DFI, DSS and PFI in TCGA cohort (Figure9c). High FLVCR1 expression was closely related to poor OS in ICGC cohort.

Potential role of FLVCR1 in HCC
Many studies reported that TIICs are closely related to prognosis of cancer patients. Therefore, we explored the relationship between FLVCR1 expression and immune in ltration using TIMER. The results indicated that FLVCR1 expression was positively associated with tumor purity and the in ltration level of B cell, CD4+ T cell, macrophage, neutrophil and dendritic cell (Figure10a). In addition, copy number alteration of FLVCR1 signi cantly correlated with in ltration level of TIICs. Elevated arm-level gain of FLVCR1 results in lower in ltration levels of CD8+ cells, macrophage, neutrophil and dendritic cell compared with normal samples (Figure10b). Furthermore, GSEA revealed that FLVCR1 was signi cantly involved in some cancer-related pathways, including cell cycle, fatty acid metabolism, glucose metabolism, MYC targets, E2F targets, WNT signaling pathway and PI3K signaling pathway (Figure10c).

Discussion
Disordered iron metabolism is one of the hallmarks of malignancies, in which iron is required for the growth, metastasis and survival of tumor cells (19,20). Previous studies have shown that iron plays a signi cant role in the development of HCC (21). HCC requires excess iron to meet the demands of rapid proliferation and DNA duplication (22). However, excessive iron can lead to ferroptosis, which has great potential in the treatment of HCC. As a commonly used chemotherapy drug for the treatment of advanced HCC, sorafenib has been found to kill HCC cells by inducing ferroptosis (23). Given the dual role of iron metabolism in HCC, it is valuable to nd new targets in iron metabolism for predicting the prognosis of HCC patients more accurately and providing new treatment strategies.
In our work, based on the expression pattern of iron metabolism-related genes, we identi ed two clusters of HCC subgroups using consensus clustering analysis. The results indicated that there is a signi cant difference in OS, stage and grade between the two clusters, which suggested that the expression pattern of iron metabolism-related genes was closely related to the development of HCC. Previous research found that iron metabolism-related genes could in uence numerous biological processes and signaling pathways, such as cell cycle (24), epithelial to mesenchymal transition (25), JAK/STAT (26) and ERK signaling pathways (27). Hence, we performed functional analyses of the differentially regulated genes between the two clusters of HCC. Similarly, the results showed that these genes are signi cantly related to cancer-related biological processes and signaling pathways. Combined with these results, it could be determined that iron metabolism-related genes play a critical effect in regulating the malignant process of HCC.
Iron metabolism is associated with the prognosis of cancer patients, including HCC patients (21). Therefore, we performed Cox regression analysis and LASSO regression analysis to construct a nine-gene prognostic signature with iron metabolism-related genes, which consists of FBXL5, SLC48A1, BMP6, HAVCR1, ALAS1, CD163, PCBP2, FLVCR1 and ABCB6. The high-and low-risk groups divided based on the median risk score showed signi cant difference in OS, and the signature had a great performance for predicating the prognosis of HCC patients. Then, the results of multivariate Cox regression analysis further indicated that this prognostic signature could independently predict the prognosis of HCC patients. In addition, a distinct difference of OS between low-and high-risk groups was observed in all subgroups strati ed by age, gender, stage and grade. Patients with high risk in all subgroups had shorter OS time than those with low risk. These results indicated that our nine-gene prognostic signature has reliable predictive performance and great clinical value.
In this nine-gene prognostic signature, FLVCR1 was identi ed as a hub gene with the highest mutation frequency, which implied FLVCR1 plays a vital role in the progression of HCC. Feline Leukemia virus subgroup C receptor (FLVCR1) functions as a mammalian cell heme exporter. FLVCR1a and FLVCR1b are two different isoforms of FLVCR1, which are expressed on the plasma membrane and mitochondria, respectively (28). FLVCR1 is involved in a variety of biological processes, including cell proliferation, cell apoptosis, oxidative stress response, cellular differentiation and metabolism (29)(30)(31). In addition, the carcinogenesis of FLVCR1 in several malignancies has gradually been uncovered. For example, Changliang Peng found FLVCR1 can promote the proliferation of synovial sarcoma through regulating autophagy (32). Although, Xianli Wei, reported the expression of FLVCR1 in HCC tumor tissues was upregulated compared with the adjacent non-tumor tissues (33), the role of FLVCR1 in progression of HCC is still unclear. In our study, GSEA showed that FLVCR1 might affect the malignant progression of HCC by regulating some biological processes such as cell cycle, metabolic pathways and cancer-related pathways. This result is consistent with the previous reports. Furthermore, we found that FLVCR1 expression was positively associated with the in ltration level of some TIICs. In the tumor microenvironment, TIICs were educated by tumor cells to adopt an iron-donor phenotype, which promotes tumor proliferation and metastasis. The mechanism by which FLVCR1 affects the progress of HCC may also regulate the iron out ow of TIICs. This requires further to explore in vitro and vivo experiments, but the interaction between iron metabolism genes and the tumor microenvironment is a potential strategy for HCC treatment.
Our study also has several potential limitations. First, this is a retrospective study, and multi-center prospective trials are still required to validate our conclusion. In addition, in vivo and in vitro experiments are needed to further verify the function and mechanism of FLVCR1 in HCC.

Conclusion
In summary, this study is the rst to use large-scale data to construct an iron metabolism-related prognostic signature that predicts the prognosis of HCC patients and reveal the potential role of FLVCR1 in the progression of HCC. These ndings could improve risk management and provide new therapeutic targets for HCC patients.      Survival risk assessment of prognostic signature in TCGA and ICGC cohort. a. Kaplan-Meier survival curves between the high-and low-risk groups. b. The distribution of survival status, risk score and corresponding expression pro les. c. ROC curves at 1, 2, and 3 years. Identi cation of the independent prognostic factors and construction of a nomogram. a, b. Univariate Cox analysis of the risk score and clinicopathological characteristics in TCGA and ICGC cohorts; c, d.

Declarations
multivariate Cox analysis of the risk score and clinicopathological characteristics in TCGA and ICGC cohorts; e. A nomogram constructed with independent prognostic factors; f. Calibration curves for nomogram model at 1,2, and 3 years; g. Heatmap and clinicopathologic characteristics of the two risk groups. *p<0.05, **p<0.01, ***p<0.001 Functional enrichment analyses, immune microenvironment and subgroup analysis for two risk groups in TCGA cohort. a. The survival curves for the low-and high-risk subgroups strati ed by age, gender, grade and stage; b. Top 5 signi cantly altered GO and KEGG pathways in high-or low-risk patients; c.
Correlation between our signature and TIICs.   Prognostic value of FLVCR1 in HCC. a. Correlation between FLVCR1 expression and clinicopathological characteristics in TCGA cohort; b. Correlation between promoter methylation levels of FLVCR1 and clinicopathological characteristics in TCGA cohort; c. FLVCR1 expression signi cantly correlated with overall survival (OS), disease-speci c survival (DSS), disease-free interval (DFI), progression-free interval (PFI).

Figure 10
Potential role of FLVCR1 in HCC microenvironment. a. Correlation between FLVCR1 expression and in ltration level of TIICs in TCGA cohort; b. Correlation between copy number alteration of FLVCR1 and in ltration level of TIICs; c. Several cancer-related pathways associated with FLVCR1 expression displayed by GSEA.