Identification of DEGs related to ferroptosis and iron metabolism in HCC
A total of 104 ferroptosis- and iron metabolism-related genes were identified to match the mRNA-sequencing data in the TCGA and ICGC databases. Using limma with an absolute log2-fold change (FC)> 1 and an adjusted P value <0.05 to perform differential expression analysis, we identified 24 DEGs (17 upregulated and 7 downregulated) in TCGA (Figure 1A and 1C) and 16 DEGs (13 upregulated and 3 downregulated) in ICGC (Figure 1B and 1D) that were related to ferroptosis and iron metabolism in HCC.
Comprehensive analysis of the ferroptosis- and iron metabolism-related genes closely associated with prognosis in HCC
We performed univariate Cox regression to explore the relationship between the expression of the 24 DEGs obtained from TCGA and prognosis using 371 HCC samples with OS rates and survival status in TCGA. Sixteen DEGs were statistically significant (P < 0.05) and considered to be associated with the prognosis of HCC. Then, LASSO Cox regression was applied to these genes. LASSO is a penalized regression method that adjusts the regression coefficient with L1 penalty to reduce the final weight of most potential indicators to zero, thereby decreasing the number of indicators with a final weight of nonzero[32]. Based on the LASSO regression with 10-fold cross-validation, we screened 7 genes with a repetition frequency greater than 900 times in 1000 substitution samplings (Figure 1E-F). Matching the 7 genes with 16 DEGs in ICGC, we finally determined that 4 genes (ABCB6, FLVCR1, SLC48A1 and SLC7A11) were significantly associated with prognosis in HCC.
Building the prognostic signature based on the four ferroptosis- and iron metabolism-related genes and validating its predictive performance
Based on a multivariate Cox regression of the four genes (ABCB6, FLVCR1, SLC48A1 and SLC7A11), we built a prognostic signature. Prognostic index (PI) = (0.135 * expression level of ABCB6) + (0.167 * expression level of FLVCR1) + (0.051 * expression level of SLC48A1) + (0.083 * expression level of SLC7A11). The optimal cut-off value was determined to be 1.4 using X-tile software and performed to divide 370 patients with HCC in the HCC cohort from TCGA into the high-risk and low-risk groups (Figure S1). The underlying diseases of HCC (including viral hepatitis, alcohol consumption, non-alcoholic fatty liver or hemochromatosis) were determined to not affect the expression profiles of these 4 genes in patients (Figure S2). OS was significantly worse in the high-risk groups than that in the low-risk groups (P<0.001, HR=3.70, 95% CI:2.22-6.25) (Figure 2 A). Figure 2C shows the distribution of risk scores corresponding to gene expression levels. The area under the curve (AUC) in the time-dependent ROC at 0.5, 1, 3 and 5 years reached 0.73, 0.77, 0.71 and 0.64 (Figure 2D), indicating great specificity and sensitivity of the prognostic signature in predicting OS. We then used the 243 HCC samples in the ICGC to validate the predictive performance of the prognostic signature. PI was calculated according to the formula mentioned earlier, and the optimal cutoff value determined by X-tile software for dividing 243 HCC samples into the high-risk group and low-risk group was 21.3. Consistent with the above results, patients with HCC in the high-risk group had a significantly lower OS than those in the low-risk group (P<0.001, HR=2.70, 95% CI: 1.49-5.00) (Figure 2B). The risk score distribution and gene expression are shown in Figure 2E. The AUCs for 0.5-, 1-, 3- and 5-year OS were 0.72, 0.67, 0.73 and 0.62, respectively (Figure 2F).
Construction and validation of the predictive nomogram in the HCC cohort from TCGA
To determine whether the predictive ability of the prognostic signature in predicting OS was independent of other traditional clinical characteristics (including age, AFP, weight, vascular tumor cell, sex, pathological grade and TNM stage), we performed univariate and multivariate Cox regression analyses on these variables using 370 HCC samples with clinical information in TCGA (Table S1). The results determined that TNM stage (HR=2.038) and risk score of the prognostic signature (HR=1.258) were independent predictive factors for predicting OS (Figure 3A). The proportional hazards of the two independent predictive factors was exhibited in Figure S3. Based on the two independent predictive factors, we constructed a predictive nomogram to quantify the prediction results of individual survival probability at 1, 3 and 5 years (Figure 3B). The C index for the nomogram was 0.66, with 1000 cycles of bootstrapping (95% CI: 0.55-0.72), and the calibration curves of the nomogram showed great consistency between the predicted OS rates and actual observations at 1, 3 and 5 years (Figure 3C-E).
We then performed ROC curve analysis to validate the predictive value of the nomogram. The AUCs for 1-, 3- and 5-year OS with the nomogram were 0.644, 0.694 and 0.667, respectively, superior to a single independent predictive factor (Figure 3F-H). To further determine the value of the nomogram in clinical decision making, we performed DCA. DCA is a new reliable evaluation tool that quantifies the clinical value of a nomogram by analyzing the clinical results obtained from the decision based on the nomogram and has important value in determining the diagnosis and adjusting the prognosis strategy[33]. We found that compared to a single independent predictive factor, the nomogram could obtain the optimal net benefit at 1, 3 and 5 years (Figure 3 I-J).
The diagnostic models were established and validated for high specificity and sensitivity
A diagnostic model integrating the four genes was established to distinguish HCC from normal subjects using a stepwise logistic regression method. Diagnostic scores were identified as follows: logit (P = HCC) = - 15.2439 + (-0.0327 × ABCB6 expression level) + (8.0880 × FLVCR1 expression level) + (3.1229 × SLC48A1 expression level) + (0.1703 × SLC7A11 expression level). Applying the diagnostic model, there was 92.00% sensitivity and 98.00% specificity in the HCC cohort from TCGA (containing 50 normal samples and paired 50 HCC samples) (Figure 4A) and 88.07% sensitivity and 92.08% specificity in the HCC cohort from ICGC (containing 202 normal samples and 243 HCC samples) (Figure 4B). ROCs in the HCC cohort from TCGA (AUC=0.980) (Figure 4C) and ICGC (AUC=0.956) (Figure 4D) were also determined to have great value in accurately distinguishing HCC from normal samples. Unsupervised hierarchical clustering of the four genes indicated a superior ability to differentiate HCC from normal samples (Figure 4E and 4F).
Since nodules less than 2 cm in the liver were difficult to distinguish from HCC through radiological or pathological examinations[34], we also constructed a diagnostic model based on the four genes in the training cohort (GSE6764) (containing 35 HCC samples and 17 dysplastic nodule samples) for differentiating nodules from HCC samples and validated it in the test cohort (GSE98620) (containing 49 HCC samples and 24 dysplastic nodule samples). Diagnostic scores were identified as follows: logit (P = HCC) = -13.9106 + (1.3676 × ABCB6 expression level) + (-0.1018 × FLVCR1 expression level) + (-0.2817 × SLC48A1 expression level) + (1.1909 × SLC7A11 expression level). The AUCs for the diagnostic model reached 0.973 in the training cohort, with 97.14% sensitivity and 94.12% specificity (Figure 5A and 5C), and 0.786 in the test cohort, with 79.59% sensitivity and 54.17% specificity (Figure 5B and 5D). Figures 5E and 5F show unsupervised hierarchical clustering of the four genes.
Comparison of the immune microenvironment of patients with HCC between the high-risk and low-risk groups
Since drugs targeting immune checkpoints have been shown to achieve antitumor effects by reversing the immunosuppressive effects of tumors, the expression of immune checkpoints has attracted widespread attention as a biomarker for identifying patients with HCC to receive immunotherapy[35]. The TMB can be used to predict the efficacy of immune checkpoint blockade and has been proven to be a biomarker for identifying patients who can benefit from immunotherapy in several cancer types[36]. In this study, we analyzed the association between risk scores and TMB. Figures 6A and 6B indicate the differences in TMB in somatic cells in patients with HCC between the high- and low-risk groups. Patients in the high-risk group had a higher TMB than patients in the low-risk group (Figure 6C). A higher OS rate was obtained in patients with low risk and low TMB group than that in patients with high risk and high TMB group (P<0.0001) (Figure 6D).
The differences in immune infiltration of 22 immune cell types obtained from 289 patients with HCC from the TCGA database are shown in Figure 7A, which may represent an intrinsic feature that can characterize individual differences. Patients with HCC in the high-risk group had higher ratios of M0 macrophages, follicular helper T cells, memory B cells, and neutrophils than those in the low-risk group (P <0.05) (Figures 7C-F). Figure 7B shows the relationship between the risk score and the expression of immune checkpoints. We found that the expression levels of CD83, B7H3, OX40 and OX40L in the high-risk group were significantly higher than those in the low-risk group (P <0.05) (Figure 7G-J), suggesting that the poor prognosis of high-risk patients was partly due to the immunosuppressive microenvironment. The results above indicated that abnormal immune infiltration and expression differences of immune checkpoints in HCC can be used as prognostic indicators and targets for immunotherapy, with important clinical significance.
Internal and external validation of the expression patterns and prognostic predictive performance of the four ferroptosis- and iron metabolism-related genes
The expression levels of ABCB6, FLVCR1, SLC48A1, and SLC7A11 were significantly higher in the HCC cohort from ICGC than in normal samples (P<0.001) (Figure 8A-D), which was consistent with the predictive analysis of diagnosis and prognosis, demonstrating that the four genes were suitable for constructing diagnostic and prognostic models. For further validation, we detected the expression characteristics of the four genes in the GSE6764 cohort. The four genes presented markedly higher expression in HCC than in dysplastic nodule samples, consistent with the findings above (Figure 8E-H). Meanwhile, we also evaluated the expression levels of the proteins encoded by these four genes in four pairs of human HCC tissues and corresponding non-tumor tissues to validate the clinical relevance of the four genes. The results of immunohistochemistry (IHC) staining confirmed that ABCB6, FLVCR1, SLC48A1 and SLC7A11 were strongly positive in HCC tissues compared with normal tissues (Figure 8I-L). In addition, the expression profiles of the four genes in multiple cell lines are shown in Figure 8M-P.
Since the four genes exhibited high expression in the tumor tissues, we explored the correlation among the genes. The expression of ABCB6 had synergy with the expression of FLVCR1, as well as the expression of ABCB6 and SLC7A11, ABCB6 and SLC48A1, and SLC48A1 and SLC7A11, which also had the same positive correlation (Figure 9A-D). The correlation between the expression of the four genes by HCC cells and the immune infiltrate is shown in Figure 9E-H. K-M curve analysis was performed to validate the predictive value of the four genes in OS. Genes with high expression had lower OS rates than those with low expression (Figure 9I-L). ROCs validated the predictive performance with high sensitivity and specificity (Figure 9M-P).
Inhibition of erastin on the proliferation and tumorigenesis of HCC and its possible molecular mechanism
As an inducer of ferroptosis, erastin was used to evaluate its influence on the development and progression of HCC[37].The chemical formula of erastin was showed in Figure 10A. Performing the CCK-8 assay, we found that erastin treatment inhibited cell proliferation in a dose-dependent manner (Figure 10B-C). And we validated that erastin treatment significantly increased the accumulation of ROS (Figure 10D-G) and iron (Figure 10H-I) in cells. In addition, Ferrostatin-1 and NAC, regarded as ferroptosis inhibitor and ROS inhibitor, were obviously rescued the anti-proliferation effect of erastin in HCC cells(Figure 10J-K) indicating that erastin induced ferroptosis to inhibit the proliferation and progression of HCC.To further evaluate the anti-tumoral effect of erastin in vivo, we constructed subcutaneous HCC xenograft models in male BALB/c nude mice by subcutaneous injection of SMMC-7721 cells. Then we treated tumor-bearing mice with erastin and vehicle, respectively. To access the potential toxicity of elastin to organs, we also performed the same elastin treatment on the non-tumor bearing male BALB/c nude mice. Figures 11A-C indicates that erastin significantly inhibits the rate of tumor volume and weight gain in mice. Importantly, we tested important organs (heart, liver, lung and kidneys) in tumor-bearing and non-tumor-bearing mice treated with erastin and confirmed that erastin treatment is nontoxic (Figure 11D). Compared with vehicle treated tumor-bearing mice, erastin-treated tumor-bearing mice did not undergo significant changes in body weight (Figure 11E). Lower expression levels of Ki67 and N-cadherin were exhibited in tumor tissues under erastin treatment (Figure 11F). Moreover, we observed that the expression of ABCB6, FLVCR1 and SLC7A11 in erastin-treated tumor tissues was significantly lower than that in vehicle-treated tumor tissues, but there was no significant difference in the expression of SLC48A1 between the two groups (Figure 11G).
As it was determined that erastin inhibited the proliferation and progression of HCC, we explored the possible molecular mechanism by which erastin achieves antitumor effects. In the Cancer Therapeutics Response Portal (CTRP) database (http://portals.broadinstitute.org/ctrp/), 52 genes were shown to be regulated by erastin, and their association is exhibited in Figure 12A. By performing Gene Ontology (GO) (Figure 12B) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses (Figure 12C) on these genes, we found that erastin could cause changes in signaling cascades, including Th17 cell differentiation and the IL-17 signaling pathway (P<0.05). This result indicated that the IL-17 signaling pathway is a potential target affected by erastin in this study.