Prognostic prediction and immune infiltration analysis based on ferroptosis and EMT state in hepatocellular carcinoma

Background: Ferroptosis is one of the main mechanisms of sorafenib against hepatocellular carcinoma (HCC). Epithelial-mesenchymal transition (EMT) plays an important role in the heterogeneity, tumor metastasis, immunosuppressive microenvironment, and drug resistance of HCC. However, there are few studies looking into the relationship between ferroptosis and EMT and how they may affect the prognosis of HCC collectively.

Methods: We downloaded gene expression and clinical data of HCC patients from the Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC) databases for prognostic model construction and validation respectively. The Least absolute shrinkage and selection operator (LASSO) Cox regression was used for model construction. The predictive ability of the model was assessed by Kaplan–Meier survival analysis and receiver operating characteristic (ROC) curve. We performed the expression profiles analysis to evaluate the ferroptosis and EMT state. CIBERSORT and single-sample Gene Set Enrichment Analysis (ssGSEA) methods were used for immune infiltration analysis.

Results: A total of thirteen crucial genes were identified for ferroptosis-related and EMT-related prognostic model (FEPM) stratifying patients into two risk groups. The high-FEPM group had shorter overall survivals than the low-FEPM group (p<0.0001 in the TCGA cohort and p<0.05 in the ICGC cohort). The FEPM score was proved to be an independent prognostic risk factor (HR>1, p<0.01). Furthermore, the expression profiles analysis suggested that the high-FEPM group appeared to have a more suppressive ferroptosis status and a more active EMT status than the low- FEPM group. Immune infiltration analysis showed that the myeloid-derived suppressor cells (MDSCs), and regulatory T cells (Tregs) were highly enriched in the high-FEPM group. Finally, a nomogram enrolling FEPM score and TNM stage was constructed showing outstanding predictive capacity for the prognosis of patients in the two cohorts.

Conclusion: In conclusion, we developed a ferroptosis-related and EMT-related prognostic model, which could help predict overall survival for HCC patients. It might provide a new idea for predicting the response to targeted therapies and immunotherapies in HCC patients.

1 Introduction

Primary liver cancer is the sixth most common cancer and the third leading cause of cancer death worldwide. As the most common pathological type, hepatocellular carcinoma (HCC) causes a significant cancer burden worldwide (1). Chronic hepatitis B or C virus infection or excessive alcohol consumption are the main causes of HCC (2). With the progress in prevention and systemic treatment of HCC, the incidence and mortality of liver cancer have gradually decreased in recent years (3). Nevertheless, high postoperative recurrence rate and low 5-year survival rate remain clinical problems (4). In-depth understanding of the molecular mechanism of HCC progression and screening out potential tumor-specific markers are of great importance for early diagnosis, treatment and especially prognostic prediction of HCC patients. Ferroptosis is a novel form of regulatory cell death mediated by iron-dependent lipid peroxidation (5), which is closely related to a variety of human diseases. Since it was first defined in 2012, ferroptosis has attracted great attention of medical researchers especially in the field of cancer treatment. Sorafenib is a first-line therapy for advanced HCC. Sorafenib, as a multi-kinase inhibitor, can not only inhibit the proliferation and angiogenesis of tumor cells by inhibiting RAF kinase and several receptor tyrosine kinases (6), but also exert cytotoxic effects by inducing HCC apoptosis (7). In addition, sorafenib is associated with ferroptosis which is considered as a new treatment mechanism (8). Studies have shown that there are multiple resistance mechanisms of ferroptosis in the treatment of HCC with sorafenib, which may lead to the decreased therapeutic effect of sorafenib (9, 10). Therefore, it is necessary to understand the ferroptosis status of different HCC patients to improve the effect of individualized targeted therapies.

EMT is known for participating in various pathophysiological processes like tissue fibrosis and cancer progression (11). With the deepening of research, it has been realized that EMT is a reversible and dynamic process in which epithelial cells acquire mesenchymal phenotype and behavior through dedifferentiation (12). In addition to increasing cell invasiveness and promoting tumor metastasis, EMT is also involved in cancer stem cells (CSC), tumor immune escape and drug resistance, which is considered as a multifunctional tool for tumor progression (13–15). EMT accounts for the complex cellular heterogeneity of HCC to some extent (12), and has been regarded as the main cause of secondary sorafenib resistance in several studies (16, 17). Targeting EMT is expected to be a new therapy to overcome drug resistance in HCC.

As two important biological processes of cancer, the potential relationship between ferroptosis and EMT has attracted more and more attention. Previous study showed that the activation of ZEB1 during EMT promoted the synthesis of polyunsaturated fatty acid (PUFA), and the metabolism of PUFA relied on GPX4. This made cancer cells vulnerable to ferroptosis (18). Subsequently, researchers used human adrenal cortical carcinoma cell line SW13 as a model to explore the link between ferroptosis and EMT (19). They found that treatment with histone deacetylase (HDAC) inhibitors could induce EMT in SW13 cells and increase their sensitivity to ferroptosis by altering intracellular iron levels. In addition, EMT-driven Discoidin Domain Receptor Tyrosine Kinase 2 (DDR2) upregulation stimulated ferroptosis susceptibility in recurrent breast cancer through the Hippo pathway (20). However, the regulatory effect of ferroptosis on EMT has been inconsistent in some studies. A single-cell RNA sequencing study found that the tendency of ferroptosis positively regulated EMT in lung adenocarcinoma (LUAD) epithelial cells (21). In contrast, ferritinophagy-mediated ferroptosis inhibited EMT in gastric cancer cells and HepG2 cells (22, 23). Exploring the crosstalk between ferroptosis and EMT could help develop novel cancer therapies (24).

Increasing evidence suggested there was a complex relationship between ferroptosis and EMT. In this study, we attempt to develop a prognostic model combining ferroptosis-related genes (FRGs) and EMT-related genes (ERGs). A new ferroptosis-related and EMT-related prognostic model (FEPM) was constructed and of great value to predict the prognosis of HCC patients.

2 Materials and methods

2.1 Data collection and genes acquisition

The level-3 mRNA expression data and corresponding clinicopathological information of HCC patients were downloaded from the Cancer Genome Atlas (TCGA) website¹ and the International Cancer Genome Consortium (ICGC) website². After excluding patients lacking crucial clinical information, the training cohort (TCGA) comprised 365 patients with HCC. According to the specific sample naming rules of TCGA database, transcriptome sequencing results of 365 tumor samples and 50 normal tissue samples were extracted from gene expression data files. The validation cohort (ICGC) included 231 patients with HCC from Japan. Accordingly, we obtained the gene expression data from 231 tumor samples and 199 normal tissue samples. Table 1 showed the clinicopathological characteristics of the two cohorts. Log₂ (TPM+1) was used to normalize gene expression data. 259 FRGs were derived from FerrDb database³. 1011 ERGs which undertake the function of encoding proteins were selected from dbEMT 2.0 database⁴.

TABLE 1

Table 1 Baseline characteristics of HCC patients in the TCGA and ICGC cohorts.

2.2 Screening crucial genes for the model

To begin with, the gene expression data of TCGA dataset was analyzed by the R package “DESeq2” (25). The gene count matrix was prepared as the input file for differentially expression analysis. The comparison was about tumor samples versus normal tissue samples. The thresholds for screening differentially expressed genes (DEGs) were set at |Log2FC| > 1 and adjusted p< 0.05. As the Venn diagram showed, the overlapping genes were identified as differentially expressed FRGs and ERGs. Thereafter, prognostic genes were recognized by univariate Cox regression analysis. The normalized expression values of differentially expressed FRGs and ERGs from tumor samples were analyzed as covariates in the univariate Cox regression. To avoid multicollinearity, FRGs and ERGs with prognostic significance (p< 0.0005) passed through the LASSO regression for variable selection and shrinkage. This analysis process generated the crucial genes participating in model construction and their corresponding coefficients by means of the “glmnet” package (26). Finally, the FEPM score was obtained according to the following formula where Coef_i is the coefficient and Exp_i is the expression value of each crucial gene.

$$FEPM score=\sum _{i=1}^{\mathit{\text{n}}}Coe{f}_i \times Ex{p}_i$$

2.3 Functional analysis

To explore the enriched pathways of different subgroups, Gene Set Enrichment Analysis was performed in GSEA software (version 4.2.2). Significantly enriched pathways were satisfied with the following standards simultaneously: nominal p value< 0.05, false discovery rate (FDR)< 0.25, and normalized enrichment score > 1. The gene set “c2.cp.kegg.v7.5.1.symbols. gmt” was selected as a reference gene set.

The prognostic genes with statistically significance (p< 0.0005) were uploaded into the STRING database⁵, for the purpose of a protein–protein interaction (PPI) network. Cytoscape (version 3.9.0) was utilized to visualize the network. The clustering effect of the prognostic model was evaluated by principal component analysis (PCA).

2.4 Analysis of immune cell infiltration

The enrichment scores of immune cells were calculated according to the single-sample GSEA (ssGSEA) approach. This analysis was implemented via the R package “GSVA”. The phenotype feature list of immune cells was retrieved from the previous pan-cancer analysis publication (27). The immune infiltration data that were generated by the CIBERSORT algorithm can be found on the GDC website⁶. A validated leukocyte gene signature matrix (LM22) was used for the deconvolution in the CIBERSORT analysis (28).

2.5 Statistical analysis

R software (version 4.1.2) and IBM SPSS Statistics (version 21.0) were utilized for statistical analysis. The “survminer” package was used to perform univariate and multivariate Cox regressions. The hazard function h(t) meaning the risk of dying at survival time “t” was used as the response variable of Cox model. FEPM score and the clinical features including gender (female vs. male), age (> 60 vs.<= 60), histological grade (G3-G4 vs. G1-G2), Child-Pugh classification (B-C vs. A), liver fibrosis (Yes vs. No), α-fetoprotein (AFP > 400 vs.<= 400 ug/L), vascular invasion (Yes vs. No), and TNM stage (III-IV vs. I-II) were used as the covariates. The survival analysis for different subgroups was displayed in the Kaplan–Meier curves. The differences between the curves were compared utilizing the log-rank test. The count data were compared utilizing the Chi-square test. The significance of differences in specific gene expressions or immune cell infiltration was assessed by Wilcoxon test. The “rms” package was used for construction and verification of nomogram. The “ggDCA” package was used for decision curve analysis. The “ggplot2”, “ggpubr”, “ggbiplot”, “pheatmap” packages were applied to visualization.

3 Results

3.1 A 13-gene signature showing high accuracy of survival prediction in the TCGA cohort

We made a flow chart to refine the main idea of this study (Figure 1). A total of 1208 ERGs and FRGs were included in this study, among which 336 genes were differentially expressed genes (DEGs) between tumors and normal tissues (Figures 2A, B). Subsequent univariate Cox regression analysis identified 38 candidate prognostic genes. PPI network of these genes showed an intricate interaction between FRGs and ERGs (Figure 2C). The LASSO Cox regression analysis identified 13 crucial genes for ferroptosis-related and EMT-related prognostic model (FEPM) based on the optimal value of λ (Figures 2D, E). The expressions of the 13 genes were as shown in the heat map (Figure 2F). Among them, there were 9 ERGs (PPARGC1A, MMP1, EZH2, STC2, KRT17, SPP1, BSG, MYCN and HOXD9), 3 FRGs (SLC7A11, STMN1 and SRXN1) and one gene (SQSTM1) belonging to ERGs as well as FRGs. The risk score for each patient was calculated based on the formula mentioned above. The TCGA cohort was classified as a high-FEPM group (182 patients) and a low-FEPM group (183 patients), regarding the median risk score as the cutoff value (Figure 3A). As illustrated in Kaplan-Meier curves, patients in the high-FEPM group had significantly shorter overall survivals (p< 0.0001) than those in the low-FEPM group (Figure 3C). Patients with higher risk scores died earlier and survived for a shorter time (Figure 3B). The area under the curve (AUC) reached 0.820 at 1 year, 0.753 at 3 years, and 0.737 at 5 years according to time-dependent ROC curves (Figure 3D). PCA showed that patients could be distinctly separated into two subgroups (Figure 3E).

FIGURE 1

Figure 1 Study flow chart.

FIGURE 2

Figure 2 Screening crucial genes for the model. (A) A volcanic plot of DEGs. (B) A Venn diagram indicating that 336 differentially expressed FRGs and ERGs were identified in the TCGA cohort. (C) A PPI network suggesting the relationship between prognostic FRGs and ERGs. (D, E) LASSO regression was performed to generate the crucial genes. (F) A heatmap showing the expression of the 13 crucial genes in the tumors and normal tissues of the TCGA cohort.

FIGURE 3

Figure 3 Survival analysis of FEPM in the TCGA cohort. (A, B) Distributions of FEPM scores and survival status of HCC patients in the TCGA cohort. (C)Kaplan–Meier curves showed that patients in the high-FEPM group had worse prognosis than those in the low-FEPM group in the TCGA cohort. (D) ROC curves of FEPM for predicting the 1/3/5-years survival in the TCGA cohort. (E) PCA plots of the TCGA cohort based on FEPM.

3.2 Validation of the prognostic model in the ICGC cohort

To avoid overfitting, the ICGC cohort were enrolled in the study to validate the model. Risk scores were figured out according to the same calculation in TCGA cohort. The ICGC cohort was classified as a high-FEPM group (115 patients) and a low-FEPM group (116 patients) (Figure 4A). Consistent with the results of the TCGA cohort, the prognosis of patients in the high-FEPM group was worse than that in the low-FEPM group (p< 0.05) (Figure 4C). More death events occurred in the high-FEPM group (Figure 4B). Because only two patients survived over 5 years, the 5-years AUC was not calculated in the ICGC cohort. The AUC reached 0.664 at 1 year, 0.663 at 2 years, and 0.670 at 3 years according to ROC curves (Figure 4D). PCA revealed distinct clustering of patients in the two subgroups. (Figure 4E).

FIGURE 4

Figure 4 Survival analysis of FEPM in the ICGC cohort. (A, B) Distributions of FEPM scores and survival status of HCC patients in the ICGC cohort. (C)Kaplan–Meier curves showed that patients in the high-FEPM group had worse prognosis than those in the low-FEPM group in the ICGC cohort. (D) ROC curves of FEPM for predicting the 1/2/3-years survival in the ICGC cohort. (E) PCA plots of the ICGC cohort based on FEPM.

3.3 FEPM score is an independent prognostic factor

As shown in Table 2, PPARGC1A was a protective factor for HCC patients (HR< 1, p< 0.001), and the rest genes were risk factors (HR > 1, p< 0.001). The heat map showed different expression levels of these genes in the two subgroups of TCGA cohort (Figure 5A). With the ascending of FEPM score, the expression level of PPARGC1A decreased but the expression levels of other genes increased. Moreover, it is interesting to find an evident correlation between FEPM score and clinicopathological features including TNM stage and tumor grade.

TABLE 2

Table 2 The 13 genes in the LASSO model.

FIGURE 5

Figure 5 (A) A heatmap showed the relationship between FEPM score and clinicopathological features in the TCGA cohort. (B-E) Univariate and multivariate Cox regression analyses proved that FEPM score was an independent prognostic factor in the TCGA cohort and ICGC cohort.

To further demonstrate the relationship between FEPM and clinicopathological features, we contrasted the baseline data between the two subgroups using the chi-square test. Patients in the high-FEPM group exhibited more advanced TNM stages (III-IV, p< 0.001), higher histological grades (G3-G4, p< 0.001), higher alpha-fetoprotein levels (> 400, p = 0.003) and higher vascular invasion possibilities (p< 0.001) in the TCGA cohort. A later TNM stage was also observed in the high-FEPM group in the ICGC cohort (III-IV, p = 0.006) (Table 3).

TABLE 3

Table 3 Relationships between FEPM and clinicopathological features of HCC patients.

Afterwards, FEPM score and clinicopathological features were submitted to univariate Cox regression analysis as predictors. The forest plots showed that FEPM score was a significant prognostic risk factor for HCC patients, either in TCGA cohort (HR = 3.611, 95%CI: 2.741-4.758, p< 0.001) or ICGC cohort (HR = 1.754, 95%CI: 1.242-2.479, p = 0.001) (Figures 5B, D). Finally, the statistically significant predictors were submitted to multivariate Cox regression analysis. As shown in forest plots, FEPM score was an independent prognostic factor for HCC patients (TCGA cohort: p< 0.001; ICGC cohort: p = 0.005) (Figures 5C, E).

3.4 Analysis of ferroptosis and EMT status

To analyze the ferroptosis status, we contrasted the expression profiles of suppressors of ferroptosis (SOFs) in the two subgroups. FTH1, GPX4, HELLS, ATF4, OTUB1, CA9, HSPB1, CD44, HMOX1, SLC7A11, SQSTM1, HSF1, ACSL3, HSPA5, SCD, NQO1 are SOFs that have been thoroughly researched (29, 30). In the TCGA cohort, it was amazing to notice that the expressions of all these SOFs were significantly upregulated in the high-FEPM group (Figure 6A). Except for GPX4 and HSPA5, SOFs also expressed at a higher level in the high-FEPM group of ICGC cohort (Figure 6B). In addition, the expressions of some drivers of ferroptosis (DOFs) were downregulated in the high-FEPM group, including CDO1, ACO1, GOT1, MAP1LC3A and PEBP1 (Supplementary Figures 1A, B). These findings suggested that the high-FEPM group appeared to have a suppressive ferroptosis status.

FIGURE 6

Figure 6 Analysis of Ferroptosis and EMT status. (A, B) Boxplots showing the expression profiles of SOFs in the TCGA and ICGC cohorts. (C, D) Boxplots showing the expression profiles of EMT markers in the TCGA and ICGC cohorts. (E) Pathways enriched in the low-FEPM group in the TCGA cohort. ^* p< 0.05; ^**p< 0.01; ^***p< 0.001; ^****p< 0.0001; ns, not significant.

GSEA analysis was accomplished in the TCGA cohort. As a result, many ferroptosis-related metabolism approaches (including primary bile acid biosynthesis, fatty acid metabolism, Glycine, serine and threonine metabolism, retinol metabolism, and cytochrome P450 activity) were significantly enriched in the low-FEPM group (Figure 6E).

To clarify whether FEPM could identify EMT status, we compared the expressions of a group of widely recognized epithelial markers, mesenchymal markers, and EMT transcription factors (EMT-TFs) between different groups (12, 13, 31). In the TCGA cohort, we found that the expressions of mesenchymal markers (CDH2, VIM, FSP1, CTNNB1, and FN1) and EMT-TFs (TWIST1, TWIST2, ZEB1, and ZEB2) were considerably increased in the high-FEPM group (Figure 6C). In the ICGC cohort, the high-FEPM group also showed great elevation in the expressions of mesenchymal markers (FSP1, CTNNB1, and FN1) (Figure 6D). These findings suggested that the high-FEPM group might have a more active EMT status.

3.5 Immunosuppressive microenvironment in the high-risk group

To find out how the FEPM reflects the tumor immune microenvironment, we computed the enrichment scores of immune cells in TCGA cohort by means of ssGSEA. Then we identified two subtypes including high immune infiltration cluster and low immune infiltration cluster via the method of cluster analysis (Figure 7A). To find out whether the high-FEPM group possessed high or low immune infiltration abundance, we performed the chi-square test. The result showed no significant statistical difference (p > 0.05) (Supplementary Table 1). However, we found a significant upregulation of immunosuppressive cells, such as myeloid-derived suppressor cells (MDSCs), and regulatory T cells (Tregs) (Figures 8A, B). Besides, a group of immune checkpoints including CTLA-4, PD-L1, PD-1, BTLA, LAG-3, TIGIT, TIM-3 and VISTA highly expressed in the high-FEPM group (Figures 8C–J). CIBERSORT analysis showed higher infiltration levels of neutrophils, macrophages M0, memory B cells, follicular helper T cells, regulatory T cells (Tregs), and activated memory CD4 T cells in the high-FEPM group. The infiltration levels of resting mast cells, monocytes, resting NK cells, and resting memory CD4 T cells were higher in the low-FEPM group (Figure 7B).

FIGURE 7

Figure 7 Analysis of tumor immune microenvironment. (A) A heatmap of ssGSEA in the TCGA cohort showing the high immune infiltration cluster (cluster 1) and the low immune infiltration cluster (cluster 2). (B) Comparison between the fractions of immune cells in the high-FEPM and low-FEPM groups of the TCGA cohort via the CIBERSORT method. ^* p< 0.05; ^**p< 0.01; ^***p< 0.001; ^****p< 0.0001; ns, not significant.

FIGURE 8

Figure 8 (A, B) Violin plots showing the upregulation of immunosuppressive cells (MDSCs and Tregs) in the high-FEPM group based on the ssGSEA enrichment scores of TCGA cohort. (C-J) Violin plots showing the upregulation of a group of immune checkpoints in the high-FEPM group in the TCGA cohort. ^**p< 0.01; ^***p< 0.001.

3.6 The nomogram performs well for prognosis prediction in the two cohorts

To improve the prediction ability and facilitate clinical application, we constructed a nomogram and appraised its performance in the TCGA and ICGC cohorts. The nomogram enrolled two independent prognostic factors including FEPM score and TNM stage (Figure 9A). According to the time-dependent ROC curves, AUCs of the nomogram reached 0.826 (1-year), 0.762 (3-years), and 0.733 (5-years) respectively in the TCGA cohort (Figures 9B–D). The AUCs at 1, 2, and 3 years reached 0.803, 0.719 and 0.701 in the ICGC cohort, indicating desirable sensitivity and specificity (Figure 9E). The C-indexes of TCGA cohort and ICGC cohort were respectively 0.734 (95%CI: 0.683-0.785) and 0.713 (95%CI: 0.637-0.789). Calibration curves further showed that the predictive values of the nomogram were highly consistent with the actual values (Figures 10A–F). Decision curve analysis (DCA) showed that the nomogram could provide clinical benefits within a wide range of thresholds (Figures 10G–I).

FIGURE 9

Figure 9 (A) Nomogram based on FEPM score and TNM stage. (B-D) ROC curves for the nomogram, FEPM score, and TNM stage in the TCGA cohort at 1, 3, and 5 years. (E) ROC curves for the nomogram in the ICGC cohort at 1, 2, and 3 years.

FIGURE 10

Figure 10 (A-F) Calibration curves showed the nomogram had great consistency. (G-I) Decision curve analysis showed that the nomogram could provide clinical benefits within a wide range of thresholds.

4 Discussion

In this study, we developed a risk prediction model combining ferroptosis and EMT. Patients in the high-risk group had shorter survival times and more advanced tumor stages. Furthermore, we conducted immune infiltration analysis and established a nomogram. ROC curves, C-indexes, calibration curves and DCA analyses showed that the nomogram had a strong predictive ability and played a certain guiding role in clinical decision-making. For high-risk patients, we can increase the frequency of follow-up to monitor tumor progression and choose more aggressive treatments to benefit patients.

Consistent with the findings of two previous studies (29, 30), we discovered that patients in the high-FEPM group displayed a suppressive ferroptosis state. HSPA5, as a promoter of unfolded protein response (UPR), is beneficial to cell survival under endoplasmic reticulum stress (32). GPX4 is a key regulator of ferroptosis and lipid peroxidation (33). Studies have proved that HSPA5 can bind to GPX4 and inhibit GPX4 protein degradation, thereby inhibiting ferroptosis in pancreatic ductal adenocarcinoma (PDAC) cells and colorectal cancer (CRC) cells (34, 35). In addition, activation of PERK/ATF4/HSPA5 pathway attenuated dihydroartemisinin-induced ferroptosis in glioma cells (36). However, no significant upregulation of HSPA5 and GPX4 was observed in the high-risk group of the ICGC cohort, suggesting that the HSPA5/GPX4 pathway might not work in subgroup stratification of the ICGC cohort. GSEA results showed that biological processes such as fatty acid metabolism, amino acid metabolism, retinol metabolism, bile acid biosynthesis, and cytochrome P450 redox reaction were significantly enriched in the low-risk group. These biological processes were closely related to lipid peroxidation and reactive oxygen species (ROS) production, suggesting that the low-FEPM group presented a more active ferroptosis status.

Based on a set of EMT molecular markers, we found a more active EMT status in the high-risk group. Although mesenchymal markers and EMT-TFs were elevated in the high-risk group, epithelial markers were not significantly downregulated but highly co-expressed. Increasing evidence proved that during cancer progression, most cancer cells developed partial EMT and expressed both epithelial and mesenchymal markers, known as the so-called hybrid E/M phenotype (37). For example, tamoxifen-induced dual recombinase lineage tracing systems have found that partial EMT cells rather than full EMT cells played a major role in the metastasis of breast cancer (38). Cytokeratin 18, despite a cytoskeletal protein of epithelial cells, critically contributed to initiating TGF-β1-induced EMT (39). Besides, the co-expression of cytokeratin and vimentin in some tumor types indicated higher invasive and metastatic potential (40).

We evaluated the tumor immune microenvironment in the TCGA cohort by ssGSEA and CIBERSORT. Recently, three immunogenomic subtypes of HCC have been identified by ssGSEA, which can distinguish the prognosis of different patients (41). In this study, we also identified two subtypes of immune infiltration by ssGSEA method, representing the high and low levels of immune infiltration abundance. To our disappointment, we found that the high-FEPM group and the low-FEPM group had no significant difference in the overall immune infiltration abundance. However, the infiltration degree of immunosuppressive cells MDSCs and Tregs, together with the expression levels of a set of immune checkpoint molecules, were prominently upregulated in the high-FEPM group. In addition, the macrophages M0 was significantly accumulated in the high-FEPM group according to the CIBERSORT result. Studies have showed that macrophages M0 can polarize into M2 phenotype under the induction of some stimuli, suppressing antitumor immunity (42). Macrophages M2 are one category of the tumor associated macrophages (TAMs). GPX4 inhibitors can promote cell death in TAMs and Tregs by inducing ferroptosis, thus reversing the immunosuppressive microenvironment (43). However, the antitumor effects of activated CD8⁺Tcells will be reduced at the same time due to their vulnerability to ferroptosis. It reminds us to utilize ferroptosis inducers cautiously after evaluating the tumor immune microenvironment.

In fact, EMT is closely bound up with tumor immune suppression and immune escape. For example, Snail promoted ovarian cancer progression by upregulating CXCR2 ligands to recruit MDSCs (44). ZEB1 induced the accumulation of MDSCs by upregulating the inflammatory cytokines in breast cancer (45). Besides, EMT-TFs increased the expression of immune checkpoint ligand PD-L1 (46, 47). On the other hand, CD4 ⁺ CD25 ⁺ Tregs promoted HCC invasion by secreting high levels of TGF-β1 to induce EMT (48). TAMs could help tumor cells migrate by secreting multiple proteases (42). The feedback loop between EMT and immune suppression promoted tumor progression (15). Hence, EMT inhibitors may be effective for the improvement of immunosuppressive microenvironment.

Over the past decade, inducing ferroptosis has exhibited inestimable therapeutic potential, which attracts great interests of many researchers. For instance, researchers have found that artesunate could cooperate with sorafenib to induce ferroptosis in HCC via different mechanisms (49). However, the inevitable drug resistance arising from long-term use will greatly limit the efficacy of these drugs. In addition, anti-CTLA-4 and anti-PD-L1/PD-1 antibodies are well-known immunotherapies of HCC. The EMT-induced tumor immunosuppressive microenvironment will greatly reduce the response of patients to these immune checkpoint inhibitors (ICIs) (15). In this study, we identified a high-risk group of patients with suppressive ferroptosis status and active EMT status. They might be lack of response sensitivity to immunotherapy because of their immunosuppressive microenvironment. Although the 13-gene prognostic model for HCC has been developed with a strong ability of prognostic prediction, as a retrospective study, this study indeed has some inevitable selection bias. Further prospective studies of larger sample’s queue are needed prior to its clinical application. In summary, a ferroptosis-related and EMT-related prognostic model was developed and validated in HCC. The implications of ferroptosis and EMT for immunotherapy were discussed emphatically in the study. In the future, the combination of ferroptosis inducers and EMT inhibitors may be a promising treatment for HCC.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

Author contributions

XD and ZL contributed to the study conception, design, and the whole process of article revision. Material preparation, data collection and analysis were performed by JW, SYL, CR and DG. The first draft of the manuscript was written by JW and ZL, and all authors commented on previous versions of the manuscript. The figures were made by ZL, LHL and SXL, HG, DL, and LL critically revised the manuscript. All authors contributed to the article and approved the submitted version.

Funding

All authors declare that they have no competing interests. This study was supported by Joint project of Henan Province Medical Science and Technology Program for Tackling Key Problems (LHGJ20190035) and Key Projects of Institutions of Higher Learning in Henan Province (21A320062).

Conflict of interest

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

Publisher’s note

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

Supplementary material

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

Footnotes

1. ^ https://portal.gdc.cancer.gov/
2. ^ https://dcc.icgc.org/
3. ^ http://www.zhounan.org/ferrdb/legacy/index.html
4. ^ http://dbemt.bioinfo-minzhao.org/index.html
5. ^ http://string-db.org
6. ^ https://gdc.cancer.gov/about-data/publications/panimmune

References

1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: Cancer J Clin (2021) 71(3):209–49. doi: 10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

2. McGlynn KA, Petrick JL, El-Serag HB. Epidemiology of hepatocellular carcinoma. Hepatology (2021) 73 Suppl 1(Suppl 1):4–13. doi: 10.1002/hep.31288

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Zhou M, Wang H, Zeng X, Yin P, Zhu J, Chen W, et al. Mortality, morbidity, and risk factors in China and its provinces, 1990-2017: A systematic analysis for the global burden of disease study 2017. Lancet (2019) 394(10204):1145–58. doi: 10.1016/S0140-6736(19)30427-1

PubMed Abstract | CrossRef Full Text | Google Scholar

4. European Association for the Study of the Liver. EASL clinical practice guidelines: Management of hepatocellular carcinoma. J Hepatol (2018) 69(1):182–236. doi: 10.3389/10.1016/j.jhep.2018.03.019

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Dixon SJ, Lemberg KM, Lamprecht MR, Skouta R, Zaitsev EM, Gleason CE, et al. Ferroptosis: an iron-dependent form of nonapoptotic cell death. Cell (2012) 149(5):1060–72. doi: 10.1016/j.cell.2012.03.042

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Wilhelm SM, Carter C, Tang L, Wilkie D, McNabola A, Rong H, et al. BAY 43-9006 exhibits broad spectrum oral antitumor activity and targets the RAF/MEK/ERK pathway and receptor tyrosine kinases involved in tumor progression and angiogenesis. Cancer Res (2004) 64(19):7099–109. doi: 10.1158/0008-5472.CAN-04-1443

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Liu L, Cao Y, Chen C, Zhang X, McNabola A, Wilkie D, et al. Sorafenib blocks the RAF/MEK/ERK pathway, inhibits tumor angiogenesis, and induces tumor cell apoptosis in hepatocellular carcinoma model PLC/PRF/5. Cancer Res (2006) 66(24):11851–8. doi: 10.1158/0008-5472.CAN-06-1377

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Louandre C, Ezzoukhry Z, Godin C, Barbare JC, Mazière JC, Chauffert B, et al. Iron-dependent cell death of hepatocellular carcinoma cells exposed to sorafenib. Int J Cancer (2013) 133(7):1732–42. doi: 10.1002/ijc.28159

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Sun X, Niu X, Chen R, He W, Chen D, Kang R, et al. Metallothionein-1G facilitates sorafenib resistance through inhibition of ferroptosis. Hepatology (2016) 64(2):488–500. doi: 10.1002/hep.28574

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Liang Y, Zheng T, Song R, Wang J, Yin D, Wang L, et al. Hypoxia-mediated sorafenib resistance can be overcome by EF24 through Von hippel-lindau tumor suppressor-dependent HIF-1α inhibition in hepatocellular carcinoma. Hepatology (2013) 57(5):1847–57. doi: 10.1002/hep.26224

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Nieto MA, Huang RY, Jackson RA, Thiery JP. EMT: 2016. Cell (2016) 166(1):21–45. doi: 10.1016/j.cell.2016.06.028

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Giannelli G, Koudelkova P, Dituri F, Mikulits W. Role of epithelial to mesenchymal transition in hepatocellular carcinoma. J Hepatol (2016) 65(4):798–808. doi: 10.1016/j.jhep.2016.05.007

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Brabletz S, Schuhwerk H, Brabletz T, Stemmler MP. Dynamic EMT: a multi-tool for tumor progression. EMBO J (2021) 40(18):e108647. doi: 10.15252/embj.2021108647

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Shibue T, Weinberg RA. EMT. CSCs, and drug resistance: the mechanistic link and clinical implications. Nat Rev Clin Oncol (2017) 14(10):611–29. doi: 10.1038/nrclinonc.2017.44

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Taki M, Abiko K, Ukita M, Murakami R, Yamanoi K, Yamaguchi K, et al. Tumor immune microenvironment during epithelial-mesenchymal transition. Clin Cancer Res (2021) 27(17):4669–79. doi: 10.1158/1078-0432.CCR-20-4459

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Xia S, Pan Y, Liang Y, Xu J, Cai X. The microenvironmental and metabolic aspects of sorafenib resistance in hepatocellular carcinoma. EBioMedicine (2020) 51:102610. doi: 10.1016/j.ebiom.2019.102610

PubMed Abstract | CrossRef Full Text | Google Scholar

17. van Malenstein H, Dekervel J, Verslype C, Van Cutsem E, Windmolders P, Nevens F, et al. Long-term exposure to sorafenib of liver cancer cells induces resistance with epithelial-to-mesenchymal transition, increased invasion and risk of rebound growth. Cancer letters. (2013) 329(1):74–83. doi: 10.1016/j.canlet.2012.10.021

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Viswanathan VS, Ryan MJ, Dhruv HD, Gill S, Eichhoff OM, Seashore-Ludlow B, et al. Dependency of a therapy-resistant state of cancer cells on a lipid peroxidase pathway. Nature (2017) 547(7664):453–7. doi: 10.1038/nature23007

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Oliveira T, Hermann E, Lin D, Chowanadisai W, Hull E, Montgomery M. HDAC inhibition induces EMT and alterations in cellular iron homeostasis to augment ferroptosis sensitivity in SW13 cells. Redox Biol (2021) 47:102149. doi: 10.1016/j.redox.2021.102149

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Lin CC, Yang WH, Lin YT, Tang X, Chen PH, Ding CC, et al. DDR2 upregulation confers ferroptosis susceptibility of recurrent breast tumors through the hippo pathway. Oncogene (2021) 40(11):2018–34. doi: 10.1038/s41388-021-01676-x

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Yao J, Zhang Y, Li M, Sun Z, Liu T, Zhao M, et al. Single-cell RNA-seq reveals the promoting role of ferroptosis tendency during lung adenocarcinoma EMT progression. Front Cell Dev Biol (2021) 9:822315. doi: 10.3389/fcell.2021.822315

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Li H, Zhou W, Wei H, Li L, Wang X, Li Y, et al. Ferritinophagic flux was a driving force in determination of status of EMT, ferroptosis, and NDRG1 activation in action of mechanism of 2-pyridylhydrazone dithiocarbamate s-acetic acid. J Oncol (2021) 2021:3015710. doi: 10.1155/2021/3015710

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Guan D, Zhou W, Wei H, Wang T, Zheng K, Yang C, et al. Ferritinophagy-mediated ferroptosis and activation of Keap1/Nrf2/HO-1 pathway were conducive to EMT inhibition of gastric cancer cells in action of 2,2’-di-pyridineketone hydrazone dithiocarbamate butyric acid ester. Oxid Med Cell Longevity (2022) 2022:3920664. doi: 10.1155/2022/3920664

CrossRef Full Text | Google Scholar

24. Ebrahimi N, Adelian S, Shakerian S, Afshinpour M, Chaleshtori SR, Rostami N, et al. Crosstalk between ferroptosis and the epithelial-mesenchymal transition: Implications for inflammation and cancer therapy. Cytokine Growth Factor Rev (2022) 64:33–45. doi: 10.1016/j.cytogfr.2022.01.006

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol (2014) 15(12):550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Softw (2010) 33(1):1–22. doi: 10.18637/jss.v033.i01

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep (2017) 18(1):248–62. doi: 10.1016/j.celrep.2016.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods (2015) 12(5):453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Wen K, Yan Y, Shi J, Hu L, Wang W, Liao H, et al. Construction and validation of a combined ferroptosis and hypoxia prognostic signature for hepatocellular carcinoma. Front Mol Biosci (2021) 8:809672. doi: 10.3389/fmolb.2021.809672

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Liu Y, Zhang X, Zhang J, Tan J, Li J, Song Z. Development and validation of a combined ferroptosis and immune prognostic classifier for hepatocellular carcinoma. Front Cell Dev Biol (2020) 8:596679. doi: 10.3389/fcell.2020.596679

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Zeisberg M, Neilson EG. Biomarkers for epithelial-mesenchymal transitions. J Clin Invest (2009) 119(6):1429–37. doi: 10.1172/JCI36183

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Wang J, Lee J, Liem D, Ping P. HSPA5 gene encoding Hsp70 chaperone BiP in the endoplasmic reticulum. Gene (2017) 618:14–23. doi: 10.1016/j.gene.2017.03.005

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Yang WS, SriRamaratnam R, Welsch ME, Shimada K, Skouta R, Viswanathan VS, et al. Regulation of ferroptotic cancer cell death by GPX4. Cell (2014) 156(1-2):317–31. doi: 10.1016/j.cell.2013.12.010

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Zhu S, Zhang Q, Sun X, Zeh HJ 3rd, Lotze MT, Kang R, et al. HSPA5 regulates ferroptotic cell death in cancer cells. Cancer Res (2017) 77(8):2064–77. doi: 10.1158/0008-5472.CAN-16-1979

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Wang R, Hua L, Ma P, Song Y, Min J, Guo Y, et al. HSPA5 repressed ferroptosis to promote colorectal cancer development by maintaining GPX4 stability. Neoplasma (2022) 69(5):1054–69. doi: 10.21203/rs.3.rs-855955/v1

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Chen Y, Mi Y, Zhang X, Ma Q, Song Y, Zhang L, et al. Dihydroartemisinin-induced unfolded protein response feedback attenuates ferroptosis via PERK/ATF4/HSPA5 pathway in glioma cells. J Exp Clin Cancer Res: CR. (2019) 38(1):402. doi: 10.1186/s13046-019-1413-7

CrossRef Full Text | Google Scholar

37. Saitoh M. Involvement of partial EMT in cancer progression. J Biochem (2018) 164(4):257–64. doi: 10.1093/jb/mvy047

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Lüönd F, Sugiyama N, Bill R, Bornes L, Hager C, Tang F, et al. Distinct contributions of partial and full EMT to breast cancer malignancy. Dev Cell (2021) 56(23):3203–21.e11. doi: 10.1016/j.devcel.2021.11.006

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Jung H, Kim B, Moon BI, Oh ES. Cytokeratin 18 is necessary for initiation of TGF-β1-induced epithelial-mesenchymal transition in breast epithelial cells. Mol Cell Biochem (2016) 423(1-2):21–8. doi: 10.1007/s11010-016-2818-7

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Kuburich NA, den Hollander P, Pietz JT, Mani SA. Vimentin and cytokeratin: Good alone, bad together. Semin Cancer Biol (2021) 86(Pt 3):816–26. doi: 10.1016/j.semcancer.2021.12.006

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Zhuang W, Sun H, Zhang S, Zhou Y, Weng W, Wu B, et al. An immunogenomic signature for molecular classification in hepatocellular carcinoma. Mol Ther Nucleic Acids (2021) 25:105–15. doi: 10.1016/j.omtn.2021.06.024

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Pan Y, Yu Y, Wang X, Zhang T. Tumor-associated macrophages in tumor immunity. Front Immunol (2020) 11:583084. doi: 10.3389/fimmu.2020.583084

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Xu H, Ye D, Ren M, Zhang H, Bi F. Ferroptosis in the tumor microenvironment: perspectives for immunotherapy. Trends Mol Med (2021) 27(9):856–67. doi: 10.1016/j.molmed.2021.06.014

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Taki M, Abiko K, Baba T, Hamanishi J, Yamaguchi K, Murakami R, et al. Snail promotes ovarian cancer progression by recruiting myeloid-derived suppressor cells via CXCR2 ligand upregulation. Nat Commun (2018) 9(1):1685. doi: 10.1038/s41467-018-03966-7

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Katsura A, Tamura Y, Hokari S, Harada M, Morikawa M, Sakurai T, et al. ZEB1-regulated inflammatory phenotype in breast cancer cells. Mol Oncol (2017) 11(9):1241–62. doi: 10.1002/1878-0261.12098

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Noman MZ, Janji B, Abdou A, Hasmim M, Terry S, Tan TZ, et al. The immune checkpoint ligand PD-L1 is upregulated in EMT-activated human breast cancer cells by a mechanism involving ZEB-1 and miR-200. Oncoimmunology (2017) 6(1):e1263412. doi: 10.1080/2162402X.2016.1263412

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Chen L, Gibbons DL, Goswami S, Cortez MA, Ahn YH, Byers LA, et al. Metastasis is regulated via microRNA-200/ZEB1 axis control of tumour cell PD-L1 expression and intratumoral immunosuppression. Nat Commun (2014) 5:5241. doi: 10.1038/ncomms6241

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Shi C, Chen Y, Chen Y, Yang Y, Bing W, Qi J. CD4(+) CD25(+) regulatory T cells promote hepatocellular carcinoma invasion via TGF-β1-induced epithelial-mesenchymal transition. OncoTargets Ther (2019) 12:279–89. doi: 10.2147/OTT.S172417

CrossRef Full Text | Google Scholar

49. Li ZJ, Dai HQ, Huang XW, Feng J, Deng JH, Wang ZX, et al. Artesunate synergizes with sorafenib to induce ferroptosis in hepatocellular carcinoma. Acta Pharmacol Sinica (2021) 42(2):301–10. doi: 10.1038/s41401-020-0478-3

CrossRef Full Text | Google Scholar