Comprehensive Analysis of Oncogenic, Prognostic, and Immunological Roles of FANCD2 in Hepatocellular Carcinoma: A Potential Predictor for Survival and Immunotherapy

Hepatocellular carcinoma (HCC) is the second most fatal cancer worldwide and imposes a major health and economic burden on society, with a five-year survival rate of < 15%^[1]. Most patients lose the opportunity for radical treatment at the time of diagnosis^[2]. Many unsolved difficulties in cancer diagnosis remain; nevertheless, the effectiveness of traditional therapies and conventional molecularly targeted drugs is increasing^[3]. Therefore, novel therapeutic targets and sensitive tumor biomarkers need to be urgently identified.

Ferroptosis is intensely discussed as an emerging form of programmed cell death (PCD) activated by iron-dependent peroxide accumulation in cells^[4]. Compared with apoptosis, ferroptosis involves various crosslinks between metabolism and environmental stress^[5,6]. Although the physiological role and function of ferroptosis remain ambiguous, emerging evidence has revealed that ferroptosis imputes several neoplastic contexts. Recent studies have considered ferroptosis a hallmark of tumor development by promoting cell death and metabolic regulation^[7-11]. Several recognized genes involved in ferroptosis may pave the way for novel strategies for treating patients with HCC, especially those pertinent to immunotherapy development. Boosting the efficacy of immunotherapy has become controversial^[12,13]; however, the heterogeneity of cancer genomes prevents immunotherapy from benefiting a significant proportion of patients. Therefore, identifying reliable biomarkers that can predict therapeutic responses is crucial.

Members of the Fanconi anemia complementation group (FANC) encode the monoubiquitinated protein FANCD2, which responds to homology-directed DNA repair by localizing in nuclear foci with other proteins. However, FANCD2 in HCC has not been comprehensively analyzed. In our study, several public databases were accessed and analyzed to determine the correlation between FANCD2 expression and the prognosis of digestive system cancers. We also explored the potential association between the expression of FANCD2 and tumor immune cell infiltration levels (TILs), tumor mutational burden (TMB), microsatellite instability (MSI), DNA mutations, and various immune checkpoint genes across multiple cancer types. The TIDE algorithm was used to predict the response of FANCD2 to immune checkpoint inhibitors (ICIs). Subsequently, we identified the biological functions and pathways of FANCD2 using gene set enrichment analysis (GSEA), Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG). We considered that FANCD2 has a pan-cancer carcinogenic effect and may serve as an independent risk factor for HCC with an unfavorable indicator of immunotherapy.

In conclusion, FANCD2 has a promising prognostic and prospective predictive value for digestive system cancers and lays a foundation for future evidence-based precision medicine, assisting in correct and reasonable clinical decision-making for patients with HCC.

The Genotype-Tissue Expression [GTEx (https://www.gtexportal.org/home/datasets)] database with data from The Cancer Genome Atlas [TCGA (https://portal.gdc.cancer.gov/)] was used to examine the dysregulation of FANCD2 expression in 33 types of cancerous and normal tissues. Raw counts of RNA sequencing data (level 3) and clinical follow-up information from TCGA were obtained from the Genomic Data Commons (GDC) data portal website. We conducted all statistical analyses using the R program v4.0.3 (the Statistics Department of the University of Auckland, Auckland, New Zealand) in which all expression data were normalized via log2 (TPM +1) conversion. The distribution of FANCD2 gene expression in tumor and normal tissues was assessed, with the horizontal axis representing different tumor tissues, and the vertical axis indicating the gene expression levels. Different colors correspond to distinct groups (^*P < 0.05, ^**P < 0.01, and ^***P < 0.001). Significance was determined using the Wilcoxon test, and the asterisk represents the degree of difference (^*P).

Clinical covariates including age, sex, tumor stage, histological grade, and clinical outcome were obtained from the aforementioned websites. Immunohistochemical (IHC) staining images for validation of FANCD2 protein expression were extracted from The Human Protein Atlas (HPA) project (https://www.proteinatlas.org/) to verify whether liver hepatocellular carcinoma (LIHC) tissues and their counterparts expressed high levels of the FANCD2 protein.

The correlation between FANCD2 expression and prognosis of patients with digestive system cancers was assessed in terms of overall survival (OS), disease-specific survival (DSS), disease-free survival (DFS), and progression-free survival (PFS) rates through forest plots using the “forestplot” R package version 4.0.3 (The R Foundation for Statistical Computing, 2020). Subsequently, R package “ggsrisk” (log rank P < 0.05) was used to calculate the median value of risk scores to classify the patients into the high-risk and low-risk groups. In addition, we used least absolute shrinkage and selection operator (LASSO)-penalized multivariate Cox regression to visualize the associations. For Kaplan–Meier (K–M) curves, the hazard ratios (HR) and 95% confidence intervals (CI) were generated by log-rank tests and univariate Cox proportional hazards regression using R packages “survival” and “survminer”. Receiver operating characteristic (ROC v0.4) analysis was performed to compare the predictive accuracy of each gene and risk score. Univariate and multivariate Cox regression analyses were performed to identify the appropriate terms to build a nomogram using the “rms” R package (P < 0.05).

Based on TIMER data (https://cistrome.shinyapps.io/timer/), FANCD2 expression profiles and immune infiltrates were analyzed in pan-cancers of the digestive system. FANCD2 expression with an immune score or immune checkpoint-related genes was mapped using Spearman’s correlation analysis. As in the heat maps, the horizontal axis represents cancer types, the vertical axis represents immune scores, and the colors indicate the correlation coefficients. For reliable immune score evaluation, we used the R software v4.0.3 which integrates the two latest algorithms, TIMER and xCell, for statistical analysis (^*P < 0.05, ^**P < 0.01, ^***P < 0.001), and the asterisk represents the degree of importance (^*P).

Differentially expressed gene (DEG) analysis was based on a cut-off of adjusted P < 0.01, and |log2(fold change)| > 2 was used to define DEG between subgroups. KEGG pathway and GO biological process analyses were performed to identify the potential molecular mechanisms of these DEGs. Gene Set Variation Analysis (GSVA) was used to determine the different pathways for analyzing FANCD2 in HCC with respect to biological functions. All statistical analyses were performed using the R package ClusterProfiler. Adjusted P-values were calculated using the Benjamini–Hochberg adjustment.

TMB counts the number of insertions and deletions in repetitive gene sequences, whereas MSI counts mutations per million bases. We downloaded data only available from the TCGA database and described the correlation between quantitative variables with a normal distribution using Spearman's correlation analysis (P < 0.05). The horizontal axis in the figure represents the gene expression distribution, and the ordinate represents the expression distribution of the TMB/MSI score. The P-value, correlation coefficient, and correlation calculation method are labeled at the front of the figure.

The (TIDE) algorithm computationally predicted the ICIs response in HCC^[14], which was conducted using the R package (v4.0.3), ggplot2 (v3.3.3), and ggpubr (0.4.0), according to the gene expression profile. Assessment of two mechanisms of tumor immune evasion, dysfunction of tumor infiltration by cytotoxic T lymphocytes (CTLs) and exclusion of CTLs by immunosuppressive factors, was based on gene expression markers. Patients with higher TIDE scores likely showed signs of antitumor immune escape and exhibited a lower response rate to ICB treatment. Different colors represent different levels of FANCD2 expression. The significance between the two groups of samples was determined using the Wilcoxon test (P < 0.05).

The correlations between different expression levels of FANCD2 mRNA and its gene mutation rates in patients were analyzed using the maftools package in the R software based on the TCGA database using the chi-square test (P < 0.05). Lollipop charts of the mutated FANCD2 gene showed the somatic mutation rate its variant classification type. Oncoplots show the somatic landscape of the HCC cohort. As estimated by MutSigCV, genes were ordered in a sidebar plot by their mutation frequencies according to the −log10-transformed q-value. Waterfall plots show mutation information for each gene in each sample. Color annotations of various cancer types are shown at the bottom. The bar plot above the legend indicates the number of mutations.

Alterations in FANCD2 expression levels were estimated in both cancerous and normal tissues using Wilcoxon tests. To explore the correlation between FANCD2 expression levels and gene-related checkpoints, we conducted a Spearman correlation analysis. In the survival analysis, we calculated HR and its P-value for the diagnosis of cancer, adjusted by univariate Cox regression. Stratified survival times were calculated based on the expression levels of FANCD2 using K–M curves. Time–ROC was used to estimate the predictive efficiency of genes for OS, whereas a nomogram was used for overall recurrence. All the above-mentioned analytical methods were performed using the R software v4.0.3 (https://www.aclbi.com/). P < 0.05 was considered significant.

Recent studies based on TCGA and GTEx databases have revealed significant differences in the mRNA expression of FANCD2 in human common digestive system cancers versus adjacent normal tissues, including LIHC.

By integrating TCGA survival analysis results through single-factor Cox regression analysis, we examined FANCD2’s prognostic value in these cancers. Figure 1A demonstrates the significant association of OS (P = 0.0123), rectum adenocarcinoma (READ) (P = 0.0046), and stomach adenocarcinoma (STAD) (P = 0.0108). We also examined DSS, DFS, and PFS rates using the same methodology (Figure 1B–D). High FANCD2 expression levels may be a risk factor for digestive system cancers and is associated with poor prognosis. Moreover, the survival differences were significant only in HCC (P < 0.05); however, whether high FANCD2 expression levels exhibited oncogenic features and were associated with poor prognosis was confirmed by the following K–M curves.

Figure 1.  Association between the FANCD2 expression levels and OS (A), DSS (B), DFS (C), PFS (D) in different types of digestive system cancers. OS, overall survival; DSS, disease-specific survival; DFS, disease-free survival; PFS, progression-free survival; CHOL, cholangiocarcinoma; COAD, colon adenocarcinoma; ESCA, esophageal cancer; LIHC, liver hepatocellular carcinoma; READ, rectum adenocarcinoma; STAD, stomach adenocarcinoma; PAAD, pancreatic ductal adenocarcinoma.

The effect of FANCD2 on immune cell recruitment in the tumor microenvironment (TME) and its impact on the prognosis of HCC remains uncertain^[15]. Thus, we explored the association between FANCD2 expression and immune infiltration using the TIMER database. As shown in Figure 2A, the expression of FANCD2 was significantly associated with the abundance of infiltrating immune cells in the pan-cancer digestive system. We also explored the subtypes of these immune cells using the online tool, xCell. Among the 38 immune cell subtypes, Th2 CD4+ T cells and common lymphoid progenitors were significantly positively associated with FANCD2 expression in these cancers. FANCD2 expression positively and significantly correlated with TILs in cholangiocarcinoma (CHOL), liver hepatocellular carcinoma (LIHC), and pancreatic ductal adenocarcinoma (PAAD) cells. An immunoinfiltration heat map of the tumor samples was used to visualize the infiltration of these immune cells (Figure 2B).

Figure 2.  FANCD2 expression significantly correlates with the infiltration levels of various immune cells in the TIMER database (A), xCell (B), and checkpoint gene expression (C) of digestive system cancers. CHOL, cholangiocarcinoma; COAD, colonadenocarcinoma; ESCA, esophageal cancer; LIHC, liver hepatocellular carcinoma; READ, rectumadenocarcinoma; STAD, stomach adenocarcinoma; PAAD, pancreatic ductal adenocarcinoma. ^*P < 0.05, ^**P < 0.01, and ^***P < 0.001.

To determine whether FANCD2 expression correlates with the expression of these checkpoint genes, we conducted a correlation analysis using the TCGA database. In LIHC, STAD, COAD, and PAAD, FANCD2 revealed a strong correlation (P < 0.05) with the expression of checkpoint genes including CTLA4 and CD274 (Figure 2C). These findings, particularly for HCC, suggested that FANCD2 levels were positively correlated with CTLA4 expression (P < 0.001, Spearman = 0.366). Thus, FANCD2 could modulate immune checkpoint activity and predict tumor immune responses.

Data available from the TCGA and GTEx databases demonstrated that FANCD2 was upregulated in all HCC cell lines (Figure 3A).

Figure 3.  Stick chart showing the high expression patterns of FANCD2 in all HCC cell lines (A), with immunohistochemical staining for FANCD2 in HCC from HPA (B) (50x). HCC, Hepatocellular carcinoma.

We further conducted a protein level validation of FANCD2 upregulation in HCC by downloading IHC images from the HPA database. FANCD2 protein levels were consistently higher in HCC tissues than in adjacent normal tissues (Figure 3B). Overexpression of FANCD2 was associated with HCC at both the gene and protein expression levels.

In this section, we assess the prognostic value of upregulated FANCD2 in patients with HCC. Patients were categorized into two subgroups according to the expression level of FANCD2. The FANCD2-low group exhibited a higher survival rate than the other groups, as shown in Figure 4A. The K–M survival curves revealed that high expression levels of FANCD2 were significantly associated with unfavorable OS (HR = 1.559, P = 0.0123) (Figure 4B). Subsequently, receiver operating characteristic (ROC) curve analysis was performed to evaluate the effectiveness of FANCD2 in predicting the prognosis of patients with HCC. The curve (AUC) values for predicting the 1-, 3-, and 5-year OS were 0.700, 0.629, and 0.608, respectively (Figure 4C).

Figure 4.  Clinical relevance of FANCD2 in patients with HCC from the TCGA database. For OS outcome, (A) distribution of risk score, survival status, and expression of FANCD2; (B) Kaplan–Meier survival curves for patients stratified by different expression levels of FANCD2; and (C) Time–ROCs for 1-year, 3-year and 5-year survival prediction. HCC, hepatocellular carcinoma; TCGA, The Cancer Genome Atlas; OS, overall survival; ROC, receive operating characteristic curve; AUC, area under the curve.

A nomogram based on discrimination and calibration methods for calculating OS was established to identify the clinical application of the prediction model. Univariate analysis revealed that FANCD2 expression (HR = 1.473, P < 0.001) and pTNM stage (HR = 1.376, P < 0.001) were correlated with OS rates (Figure 5A), indicating that FANCD2 was an independent predictor of HCC progression after adjusting for age, sex, race, and pathological stage in a multivariate Cox regression model (HR = 1.458, P < 0.001; Figure 5B). The C-index results reflected the relatively favorable predictive performance of the nomogram. Calibration plots indicated favorable concordance between the predictors and observers at 1, 3, and 5 years (Figure 5C, D).

Figure 5.  Evaluation of the prognostic value of overall survival in patients with FANCD2 in HCC Forest plot of hazard ratios and P-values of the constituents involved in both univariate (A) and multivariate (B) Cox regression considering clinical information and prognostic FANCD2 in hepatocellular carcinoma (HCC); (C) nomogram to predict 1-, 3- and 5-year survival rates; and (D) calibration curve for the nomogram model in HCC. The dashed diagonal line represents an ideal nomogram.

The prognosis of patients with cancer is greatly affected by immune surveillance and immune responses triggered by immune checkpoints. An increasing number of genes are tightly linked to or recognized as components of immune response checkpoint components^[16]. As ICIs have been approved for the treatment of tumors with the ability to improve their overall survival and objective response rates, they have become one of the first-line regimens for patients with HCC^[17,18]. Therefore, identifying predictive biomarker responses to ICI-based immunotherapy is essential for optimizing therapeutic outcomes. The TIDE score was a more accurate predictor of survival outcomes than the PD-L1 expression level and TMB in patients treated with ICBs^[14]. A recent study has reported its utility in predicting or evaluating the efficacy of ICB therapy^[19].

We observed that the expression of FANCD2 was positively correlated with all types of immune cells and may be involved in the regulation of the tumor immune response; hence, we investigated the relationship between FANCD2 expression and TMB/MSI. The associations between the FANCD2-high group and both biomarkers were not statistically significant (Figure 6A, B). The TIDE score was significantly higher in the FANCD2 high-expression group than in the low-expression group (Wilcoxon rank-sum test, P < 0.001; Figure 6C), which reinforced our hypothesis that patients with HCC with a high FANCD2 expression level are more probable to have poor responsiveness to ICB therapy.

Figure 6.  High FANCD2 expression levels were not significantly associated with TMB (A) or MSI (B) in HCC, and patients treated with ICB had higher TIDE scores (C) than those in the low-expression group.

To further confirm the molecular mechanisms by which FANCD2 is involved in HCC, we divided these samples into two subgroups according to the expression level of FANCD2 to perform both GSEA and GO/KEGG analysis. In contrast to the FANCD2-low group, 481 upregulated genes (such as TFRC, BRCA1, EPCAM, and AFP) and 95 downregulated genes were identified in the FANCD2-high group (such as GPX4, SLC10A1, and CYP3A4) (Supplementary Table S1).

Based on the enrichment analysis of the GO/KEGG pathway, the upregulated genes were significantly related to cell cycle regulation, whereas the downregulated genes were mainly involved in metabolism by CYP450, consequently regulating nuclear and organelle fission, DNA repair, ferroptosis, and fatty acid and xenobiotic metabolism. On this basis, we performed GSVA, annotated pathways of interest. As displayed in Figure 7 the upregulated genes were enriched in various pathways in the Hallmark database including cell cycle, Myc-targets, PI3K/AKT/ mTOR signaling, cellular response to hypoxia, and so on. These results indicate that FANCD2 plays an integral role in many aspects of tumor biology and the microenvironment.

Figure 7.  Using GSVA upregulation of FANCD2 was revealed to be associated with the activation of multiple pathways in HCC. GSVA, Gene Set Variation Analysis; HCC, hepatocellular carcinoma.

The gene mutation landscape of HCC has also been investigated. As shown in Figure 8A, compared with other cancer-promoting genes the FANCD2-caused somatic mutation rate, owing to by missense mutations, was only 0.84%. Compared with other mutated genes of the same type, in HCC FANCD2 mutation rates can counter to its expression level (Figure 8B). Among the two groups of patients with HCC with different FANCD2 mRNA expression levels, most had wild-type FANCD2 with subtle differences in mutation rates (Figure 8C).

Figure 8.  Lollipop charts (A) showed mutated sites of FANCD2, the mutational landscape (B) and box-plot (C) showed only a tiny difference in mutation rates among the two groups stratified by high-/low-expression level of FANCD2 in TCGA HCC cohort. GSVA, gene set variation analysis; HCC, hepatocellular carcinoma.

Cancer is a major public health issue worldwide, owing to its high morbidity and mortality rates. According to the latest cancer statistics, liver cancer has become one of the top three most fatal cancers globally, particularly in China. Early recognition and effective treatment are critical prerequisites for improving the prognosis, as common treatments such as surgical resection, interventional therapy, and adjuvant chemotherapy have limited efficacy to varying degrees^[20].

In diverse diseases, ferroptosis, a form of PCD independent of apoptosis, has been implicated in several pathways^[21]. However, research on the role of ferroptosis in oncology and treatment responses is still in its infancy, particularly for HCC. Previous studies^[22,23] have suggested that FANCD2 may play an oncogenic role in HCC progression and thus act as a potential biomarker for prognosis. Nevertheless, increasing FANCD2 expression leads to the upregulation of the FA pathway, and owing to DNA repair, cancer cells become more resistant to DNA-damaging anticancer agents. The mechanisms modulating ferroptosis sensitivity mainly involve tumor suppression genes, hypoxia-inducible factors, and the degree to which cells are in a mesenchymal state. Cancer cells are treatment-resistant and possess stem cell-like features^[10,24], which may be one of the major causes of targeted therapy failure in most patients with HCC.

To elucidate the detailed molecular mechanisms underlying these unexpected phenomena, we comprehensively examined the expression of FANCD2 macroscopically and microcosmically. In the pan-cancer dataset, a significant upregulation of FANCD2 in in paracancerous tissues compared with normal tissues. The expression levels are related to the abundance of TILs and immune checkpoints in various cancers of the digestive system. High expression levels of FANCD2 are correlated with worse OS, DSS, DFS, and PFS rates in several cancers, particularly HCC. We also analyzed the correlation between FANCD2 expression and HCC. We conducted mRNA-level- and protein-level validation of high-level expression of FANCD2. Most patients with HCC with high expression FANCD2 levels were accompanied with the upregulation of ferroptosis-promoting genes and downregulation of the suppressing genes, suggesting that FANCD2 plays a role in oncogenesis by promoting iron-dependent cell death. FANCD2 showed good predictive accuracy for clinical outcomes in HCC according to univariate and multivariate Cox regression analyses; as such, we identified its value as an independent risk factor. Enrichment analysis demonstrated that FANCD2 participates in various signaling pathways, which may consequently exert a pleiotropic effect on HCC malignancy in many aspects such as DNA repair, ferroptosis, hypoxia, and metabolism.

These findings support our hypothesis that FANCD2 plays a tumor-promoting role in the digestive system and influences TILs in the TME through immune cell recruitment. In addition, FANCD2 acts as an independent risk factor for HCC, has recognized value in predicting tumor aggressiveness and prognosis, and may be a potential biomarker for poor responsiveness to ICB therapy. However, this was not a prospective study, and all data analyzed were obtained from public databases; the lack of in vitro or in vivo data limits the results. Despite several mutated genes in liver cancer, very few mutations in genes that drive the cancer; the mutated genes in patients with liver cancer are too dispersed, and certain complex interaction mechanisms were observed among the individual genes, which has led to less-than-ideal data from clinical drug trials of targeted drugs against driver genes. Consequently, the data on actual clinical cases in this area are becoming increasingly less available. We will conduct clinical diagnostic and therapeutic studies on liver cancer-specific driver target genes to provide theoretical results. Clinical studies are needed to validate its prognostic value, and more in-depth experimental studies are required to reveal the mechanisms of action of FANCD2 in ferroptosis and tumor immunity.

In conjunction with ICIs and immune cell therapy, tumor immunotherapy has become a new frontier in oncology research. TMB, MSI, and DNA damage repair (DDR) have emerged as markers for tumor immunotherapy; however, < 10% of patients with digestive system tumors have benefited from ICI treatment^[17,25,26], indicating that the guiding significance of these markers for immunotherapy is limited, and the prevalence of a low treatment response rate also needs to be further explored. FANCD2, a classic DDR gene, is closely associated with DNA mismatch repair. Almost all FANCD2 genes were wild type, regardless of the mRNA expression level, and the related somatic mutation rate was < 1%. Second, in patients with HCC with high FANCD2 expression levels, the expression levels of most genes belonging to DDR were also significantly increased, possibly resulting in cancer promotion from genomic instability through DDR dysfunction. Therefore, the treatment outcome in this area may be futile. DDR gene mutations are associated with higher TMB and TILs^[27]. In immunotherapy, patients with harmful DDR mutations have higher objective response rates and better clinical outcomes than patients with DDR gene wild-type. In addition, PD-1 antibodies and targeted CTLA-4 blocking are ineffective in patients with non-MSI-H digestive system tumors^[28-32]. Therefore, patients with HCC with high FANCD2 expression may have a negative DDR and lower TMB/MSI, suggesting that these patients have a poor response to ICB.

Solid tumors consist of highly heterogeneous tumor cells and a tumor microenvironment (TME), which contains various signaling molecules and is considered a highly complex tissue. Local immunity is constantly remodeled in an immunomodulatory microenvironment. The tumor immune microenvironment (TIME) is an immune system composed of different cell lines of the immune system and their interactions in the TME niche, which has become pivotal to the occurrence, development, and treatment of cancer. TIME can be divided into invasive inflammation (as “hot tumor”) and invasive rejection (as “cold tumor”). Compared with the surface of “cold tumors”, that of “hot tumors” has more specific molecular markers, which can be recognized and attacked by T cells, thus triggering anti-tumor immune response, such as HCC^[33,34]. The systemic tumor immune environment (STIE) is primarily controlled by the circulating blood and lymphatic vessels via immune regulatory molecules and immune cells, which are crucial for tumor metastasis and immune infiltration. Recent studies have also found that it contains a new tumor antigen (neoantigen) that can inhibit tumor growth^[35]. TIME and STIE are closely related and play complementary role^[36]. Future research may be able to reshape the STIE and TIME by blocking inhibitory receptors, costimulatory receptors, cytokines, and immune checkpoints and combining them with traditional therapeutic strategies to achieve better anti-tumor effects.

We comprehensively provided a deeper understanding of the role of FANCD2 in oncogenesis in digestive system cancers for the first time. We propose that FANCD2 is an independent risk factor for HCC and may be used as a new biomarker for poor immunotherapy and worse prognosis. Although we provided novel insights into the role of the ferroptosis regulator FANCD2 in HCC for prognosis, immune infiltration, and immune therapy, its potential biological functions need to be elucidated in prospective, multicenter studies. Moreover, our study preliminarily identified a strong correlation between FANCD2 and ICIs, and a more ideal management for diagnosis and treatment could be explored in the future.