Immune infiltration and a ferroptosis-related gene signature for predicting the prognosis of patients with cholangiocarcinoma

: Cholangiocarcinoma (CHOL) is a digestive tract tumor with high malignancy and poor prognosis and is extremely challenging to treat. At present, induced cell death holds great promise in tumor therapy. Ferroptosis is a recently proposed pattern of programmed cell death, and numerous studies have shown that it is intimately involved in tumors. However, the roles of differentially expressed ferroptosis-related genes (DEFRGs) in CHOL have not been investigated. Our study was based on The Cancer Genome Atlas (TCGA) database, and DEFRGs were obtained to construct a prognostic riskScore model of CHOL by univariate and multivariate Cox regression analyses. Subsequently, the model was evaluated by nomogram construction, survival analysis, receiver operating characteristic (ROC) analysis, and exploration of the immune microenvironment. The mRNA and protein expression levels of each gene in the model were validated by the Gene Expression Omnibus (GEO) database, quantitative real-time PCR (qRT-PCR) and immunohistochemistry (IHC) staining. Our study found that the construction of a nomogram confirmed the predictive value of the model for overall survival (OS), and it was confirmed to have high diagnostic value by ROC analysis. Our experimental results were almost consistent with our bioinformatics results. In conclu-sion, we found that the prognostic model showed extremely high diagnostic and prognostic value and could predict the possibility of immunotherapy, thus providing a new direction for individualized treatment of patients with CHOL.


Introduction
Cholangiocarcinoma (CHOL) is a highly malignant solid tumor originating from the epithelial cells of the biliary tree. According to anatomic location, it is classified into intrahepatic and extrahepatic cholangiocarcinoma, and the latter is further divided into hilar and distal cholangiocarcinoma [1]. As the second most common primary liver tumor, it has currently an increasing trend in morbidity and mortality, with significant regional variation in incidence, accounting for 80% in Asia and South America [2,3]. CHOL has an insidious early onset, and the late clinical manifestations are mostly nonspecific, mainly related to biliary obstruction caused by the tumor. Most patients are often diagnosed in the middle or advanced stage or when the disease is accompanied by distant metastasis. Although there are many treatment options, the overall prognosis is still poor [4,5]. In addition, current multigenic prognostic models have been widely used to assess the risk of tumors, but a well-established prognostic biomarker for CHOL is still lacking, which makes it difficult for hepatobiliary surgeons to risk-stratify CHOL, affecting the choice of treatment regimen. Therefore, we used bioinformatics approaches to mine The Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO) database and developed new prognostic models, which are essential to improve the prognosis.
Ferroptosis is an iron-dependent regulatory cell death process driven by lipid peroxidation that is different from apoptosis and necrosis. Its characteristic morphologic features are mainly plasma membrane blistering, mitochondrial cristae reduction or disappearance, mitochondrial shrinkage and mitochondrial membrane condensation [6,7]. Since the concept of ferroptosis was first proposed in 2012, many studies have shown that ferroptosis is strongly correlated with systemic lupus erythematosus, inflammatory bowel disease, and tumors [8,9]. Among these, ferroptosis has attracted extensive attention as a potential treatment modality for tumors, and we can find a large number of studies on ferroptosis in relation to tumors such as hepatocellular carcinoma [10], gastric cancer and lung adenocarcinoma [11,12]. Currently, some ferroptosis inducers composed of small molecules or nanomaterials have been developed, which not only help to overcome drug resistance but also prevent the metastasis of tumors. However, the roles of differentially expressed ferroptosis-related genes (DEFRGs) in CHOL have not been deeply studied and should be further explored.
In our study, we thoroughly explored the roles of DEFRGs in CHOL and developed a novel prognostic model based on DEFRGs. First, RNA-sequencing (RNA-seq) data and relevant clinical data were extracted from the TCGA database. Next, the FRGs were extracted from FerrDb, the two databases were integrated, and a prognostic model of CHOL composed of four genes was constructed and subsequently validated in the GEO database (GSE26566).
University from July to November 2021. Paraffin sections of 80 pairs (N=20) of CHOL tissues and normal tissues were retrospectively collected from the Department of Pathology of the First Affiliated Hospital of Chongqing Medical University (collected from December 2019 to June 2021), and all sections were confirmed as CHOL or paracancerous tissues by two pathologists. This study was approved by the Ethics Committee of the First Affiliated Hospital of Chongqing Medical University, and all patients signed informed consent forms.

Quantitative real time PCR (qRT-PCR)
Total RNA was extracted from 10 pairs of fresh CHOL and paracancerous tissues with TRIzol reagent (Takara Biotechnology Co., Ltd., Dalian, China). Total RNA was reverse-transcribed into cDNA using a PrimeScript TM RT kit (Takara Biotechnology Co., Ltd.). Then, cDNA was amplified with specific primers, and the primer sequences involved in this study are shown in Table 1. The 2 -ΔΔCt method was used to calculate the relative expression levels of genes of interest using β-actin as a standard.

Immunohistochemistry (IHC) staining
First, 80 pairs of CHOL and paracancerous paraffin sections were sequentially deparaffinized in xylene, rehydrated in different graded ethanol solutions, and processed for antigen retrieval and peroxidase activity blockade. Finally, we further elucidated the roles of FRGs in CHOL by quantitative real-time PCR (qRT-PCR) and immunohistochemical (IHC) staining of clinical samples. Studies have shown that this model holds great promise in evaluating CHOL prognosis and guiding immunotherapy. The flow chart of our research is shown in Figure 1.

Patients
We randomly collected 10 pairs of fresh CHOL specimens and adjacent normal tissue specimens from the First Affiliated Hospital of Chongqing Medical Am J Transl Res 2022;14(2):1204-1219 Subsequently, the tissue sections were incubated with primary antibodies against STMN1 (YT3465; 1:200 dilution), ALOX5 (GB111330; 1:1000 dilution), BNIP3 (GB111204; 1:1000 dilution) and EGFR (GB111504; 1:1000 dilution). Then, sections were incubated with HRPlabeled goat anti-rabbit antibody (GB23303; 1:200 dilution) for 50 min in a 37°C incubator, stained with DAB and counterstained with hematoxylin. A semiquantitative method was used to evaluate the results of IHC staining by two pathologists blinded to the sample information. The specimens were classified into a high expression group and a low expression group according to the degree of staining, and the images were captured using a Leica microscope imaging system (100 × and 400 × magnification; Leica Microsystems GM-bh, Wetzlar, Germany).

Data collection from public databases
In our study, the RNA-seq datasets and corresponding clinical information of 36 CHOL samples and 9 normal samples were derived from the TCGA database (https://portal.gdc.cancer. gov/repository/), the data cutoff was October 10, 2021, and differentially expressed genes between tumor and normal samples were selected based on gene attributes by the "limma" package. The cutoff value was set at |log 2 fold change| >1 and FDR <0.05. In addition, we extracted FRGs from the FerrDb (http: //www.zhounan.org/ferrdb/) database, which is a database containing ferroptosis-related regulators and ferroptosis-related diseases, and then intersected the data with the above FRGs to obtain DEFRGs. The GSE26566 dataset was acquired from the GEO (http://www. ncbi.nlm.nih.gov/geo) database as a validation set, derived from the GPL6104 platform. This dataset includes 104 freshly-frozen tumors and 59 matched non-cancerous livers from Australia, Europe and the United States, and the experiments performed by Copenhagen University.

Establishment of DEFRGs prognostic model
First, univariate Cox regression analysis was applied to identify the DEFRGs with potential prognostic value, a prognostic model was constructed, and the correlation coefficient was obtained by multivariate Cox regression. Then, the riskScore of CHOL patients was determined by the following formula: RiskScore = e (β1 × expression of STMN1) + (β2 × expression of EGFR) + (β3 × expression of ALOX5) + (β4 × expression of BNIP3) -0.09 (β values for the correlation coefficient and a constant value of 0.09). Then, according to the median value of the riskScore, patients were divided into highand low-risk groups for further analysis.

Construction of the nomogram
Univariate and multivariate Cox regression analyses were used to explore whether the riskScore and clinical values could be used as independent prognostic factors for CHOL patients. Subsequently, the specificity and sensitivity of the prognostic model were verified by receiver operating characteristic (ROC) curve analysis. In addition, based on clinicopathological values and riskScore, we constructed a nomogram by using the "regplot" R package to evaluate the 1-, 3-, and 5-year overall survival of CHOL patients.

Bioinformatics analysis
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of the above DEFRGs were performed using the "clusterprofiler" package and then visualized using "ggplot". Then, we used the single-sample gene set enrichment analysis (ssGSEA) function from the "gsva" package in the R language to calculate individual enrichment scores for each paired sample and gene set [7,13], allowing quantification of the extent of immune cell infiltration in tumors. The 24 immune cells involved in this study were refer- Table 1. Sequences of the primers utilized for quantitative real-time PCR enced to a previous study by our group [14]. Furthermore, we explored immune-related pathways involved in prognostic models.

Statistical analysis
In our study, the statistical analyses were completed with R 4.1.5, SPSS (version 25), and plots of the results were performed using GraphPad Prism software (version 8.0) and R 4.1.5. Student's t test was performed to analyze the difference in gene expression between the two groups of continuous variables, and Pearson's correlation analysis was used. Univariate and multivariate Cox regressions were used to analyze the hazard ratio (HR) and 95% confidence interval (CI) and assess the prognostic value. P value <0.05 was considered significant.

Screening of DEFRGs in CHOL patients
The RNA expression data and corresponding clinical information of 36 CHOL and 9 adjacent normal samples were extracted from TCGA database, hierarchical clustering was used to detect the sample group data (Figure 2A), a total of 2197 differential genes were screened out by "limma" R package, and 259 FRGs were extracted from FerrDb. The intersection of the above differential genes was taken, and a Venn diagram was used for visualization ( Figure 2B). A total of 65 DEFRGs were obtained, and these up-and downregulated genes are presented in a volcano plot ( Figure 2C).

Construction of the ferroptosis-related prognostic model
Univariate Cox regression analysis was applied to evaluate the potential prognostic value of the above DEFRGs. As a result, four DEFRGs related to the overall survival of CHOL, namely, STMN1, ALOX5, EGFR, and BNIP3, were selected for the construction of the prognostic model ( Figure 2D), and multivariate Cox regression analysis was used to calculate the correlation coefficient ( Figure 2E). The riskScore formula is as follows: RiskScore = e (0.009000742 × expression of STMN1) + (-0.137555962 × expression of EGFR) + (-0.023056411 × expression of ALOX5) + (0.191038641 × expression of BNIP3) -0.09 . In addition, the patients were divided into a high-risk group (n=18) and a low-risk group (n=18) according to the median riskScore ( Figure 3A). As shown in the figure survStat ( Figure 3B), patients in the high-risk group had a significantly shorter survival time, and Kaplan-Meier survival analysis further confirmed our conclusion (P<0.001, Figure 3C). In addition, the ROC curve indicated that the riskScore model had higher diagnostic value (AUC=0.991, Figure  3D).

Verification of the prognostic model
First, we established a nomogram by combining the clinical data (sex and stage) of patients and the riskScore to further predict the 1-, 3-, and 5-year overall survival of CHOL patients ( Figure 4C).
In addition, of the 45 samples in the TCGA database, STMN1 and ALOX5 were significantly upregulated compared with normal tissues (all P<0.05), but the expression of BNIP3 was significantly downregulated compared with normal tissues, and there was no significant difference in the expression of EGFR between CHOL and adjacent normal tissues (P<0.05, Figure 5A). In addition, to increase the reliability of genes in the model, we downloaded data on 169 cases of CHOL and adjacent sample data containing gene expression from the GSE26566 cohort as the validation set. Unlike the results obtained from TCGA, EGFR was highly expressed in normal samples (P<0.05, Figure 5B). Moreover, we examined the mRNA and protein expression levels corresponding to the four genes between CHOL tissues and adjacent normal tissues by qRT-PCR (n=10, P<0.05, Figure 5C) and IHC (n=20, Figure 6). We found that the results were almost consistent with those in TCGA and GEO, and the expression level of EGFR was low in CHOL tissues, the results indicated that the difference was more significant after expanding the sample size.
Finally, the survival curves of the above genes were plotted using the Gene Expression Profiling Interactive Analysis (GEPIA) database; patients with low expression of STMN1 and  BNIP3 had a long survival time, and those with high expression of ALOX5 versus EGFR had a long survival time ( Figure 5D-G).

Functional enrichment analysis
We performed functional enrichment analysis of 65 DEFRGs to explore the pathways and biological functions involved in CHOL. The GO results showed that these DEFRGs were mainly related to the response to oxidative stress ( Figure 7A). KEGG analysis indicated that these DEFRGs were mainly enriched in ferroptosis, glutathione metabolism, 2-oxocarboxylic acid metabolism and central carbon metabolism in cancer ( Figure 7B).

Correlation between immune infiltration and prognosis of CHOL
To explore the infiltration of immune cells in CHOL, we analyzed 45 CHOL and paracancer-ous samples in TCGA, which were divided into a high invasion group and a low invasion group according to the cluster analysis of the degree of immune cell infiltration, and the immunoprofile of CHOL was visualized using ssGSEA ( Figure 7C). Based on this analysis, we found significant differences in the infiltration of immune cells in CHOL. Next, the association between the four DEFRGs and the immune infiltrate microenvironment was further analyzed ( Figure 8A)  cells. Tgd, cytotoxic cells and neutrophils were more significantly enriched in the high-risk group of the prognostic model.
Next, we validated the correlation between immune cells and the prognostic model by immune cell markers, which showed that BNIP3 was negatively correlated with SPON2 and IL-21R, markers of NK CD56dim cells, and strongly positively correlated with CYP4F3 and VNN3, markers of neutrophils ( Figure 8C).
ALOX5 was positively correlated with the molecular markers LAMP3 and OAS3 of aDCs and negatively correlated with the immune marker RORC of Th17 cells ( Figure 8D). STMN1 was positively correlated with the biomarkers BIRC5, CDC25C, HELLS, NEIL3 and CENPF of Th2 cells ( Figure 8E), which were consistent with our previous conclusion.
In addition, enrichment analysis of immune function in CHOL showed that the riskScore was negatively correlated with the IFN-γ signature. Then we explored the pathways related to CHOL by GSEA. The riskScore was positively correlated with tyrosine metabolism, fatty acid metabolism, and cytochrome P450 drug metabolism, but negatively correlated with complexity and coordination cascades ( Figure  8F).

Discussion
Cell death has long been considered a typical feature of malignant tumors, as it is very important for maintaining body balance and preventing excessive tumor proliferation. Cancerous cells always show high metabolic levels, antioxidant modification, and increased iron intake to maintain unlimited proliferation [15]. The proposal of ferroptosis alters people's understanding of cell death and emphasizes the important role of excessive iron metabolism in carcinogenesis [16]. In recent years, this unique cell death mode has become a new research hotspot. A considerable number of scholars have explored the mechanism of ferroptosis in tumors, such as hepatocellular carcinoma [17], glioma [18], and lung adenocarcinoma [19]. However, the mechanism of ferroptosis in CHOL is not clear. In our study, we found that these DEFRGs were mainly related to oxidative stress, ferroptosis, and tumor metabolism. In addition, a nomogram based on clinicopathological values and riskScore was constructed to help clinicians make better decisions and treatment plans. The ROC curve indicated that the model also has high specificity and sensitivity and can be used for the early diagnosis of CHOL. In addition, single-sample immune cell enrichment analysis showed that ferroptosis was closely related to immunity in CHOL. These results suggest that the model has great research value in CHOL.
At present, although great progress has been made in the treatment of CHOL, its prognosis is still poor due to its complex pathogenesis. Therefore, understanding the specific pathogenesis of CHOL is very important for the development of targeted therapy. The GO and KEGG results suggested that the DEFRGs were significantly enriched in ferroptosis, reactive oxygen species (ROS), cell oxidative stress, iron ion reaction, carbon metabolism, and other related pathways. These biological processes are crucial for the implementation of ferroptosis [20,21]. A large number of studies have shown that ferroptosis is a ROS-and iron-dependent cell death mode [22]. Therefore, the close relationship between ferroptosis and oxidative stress can be understood as iron-dependent lipid peroxidation accumulation. Iron chelators and fatsoluble antioxidants can inhibit ferroptosis [23]. We further analyzed the gene enrichment pathway in CHOL by GSEA, which showed that it was significantly enriched in amino acid, fatty acid and other metabolic pathways, and the riskScore was positively correlated with most pathways ( Figure 8F). Moreover, studies have shown that fatty acids can be used as an energy source for the progression of CHOL [24], which further confirms our conclusion.
It is worth noting that the prognostic model we constructed is closely related to the ferroptosis of tumors. The mitotic receptor BNIP3 is a key mitochondrial autophagy receptor in cells and is necessary to prevent increases in ROS levels and cell death. It is regulated by upstream hypoxia inducible factor (HIF-1) and is highly expressed in normal tissues; it can inhibit the growth of many types of tumors [25,26]. However, when cells become cancerous, the interaction between inflammation and tumors increases the production of ROS in tumor cells [25,27], and the production of ROS may also be carcinogenic, leading to oxidative DNA damage and genetic instability [28]. At this time, the BNIP3 gene is a risk factor in tumors, which may be the result of the joint action of many factors. The lipase family encoded by the ALOX5 gene plays an important role in leukotriene biosynthesis and participates in a variety of processes, such as inflammation and tumors [25].
In the study of stroke, some scholars proposed that N-acetylcysteine (NAC) can prevent ferroptosis by neutralizing the toxic lipids produced by ALOX5 [29]. However, the role of ALOX5 in CHOL needs to be further explored. STMN1 is the main regulator of microtubule dynamics and plays an important role in the regulation of the cell cycle [30], which is closely related to the proliferation and division of tumor cells [31].
In recent years, there has been a close relationship between immunity and tumor progression, and this is also a research hotspot. A large number of studies have shown that the tumor immune microenvironment has important value in predicting tumor prognosis and evaluating therapeutic effects [32]. Different types of immune cell infiltration may create a favorable or unfavorable environment for tumor cell proliferation and metastasis [33]. Interestingly, this study found a correlation between the expression level of genes in the prognostic model and immune infiltrating cells.
Among them, BNIP3 has a strong correlation with most immune cells, especially NK CD56dim cells and neutrophils. First, for NK CD56dim cells, increasing evidence supports that NK CD56dim cells play a cytotoxic role mainly by mediating the antitumor response [34,35]. We also found that the riskScore was negatively correlated with IFN-γ in CHOL (Figure 8F). Notably, NK CD56dim cells secrete IFN-γ to activate antitumor adaptive immunity in head and neck cancer [36], further confirming our conclusion. Second, it is well known that the presence of neutrophils in tumors is usually related to the poor prognosis of patients because they participate in the proliferation, invasion, and immunosuppression of various tumors and induce blood vessels by releasing ROS, reactive nitrogen species (RNS) or protease to promote tumorigenesis [37][38][39]. This is consistent with the aforementioned poor prognosis of patients with high expression of BNIP3 in CHOL ( Figure 5D).
In addition, we found that the expression of ALOX5 was positively correlated with aDC and its molecular markers LAMP3 and OAS3. From the cluster analysis of immune cells ( Figure  8B), aDCs were positively correlated with NK CD56dim cells. Therefore, aDCs may also play an antitumor role in tumors. LAMP3 is induced by hypoxia and is upregulated in esophageal cancer [40]. Another molecular marker, OAS3, is an important part of interferon-induced antiviral enzymes, and studies have shown that the high expression of OAS3 in breast cancer is closely related to prognosis [41]. The expression of ALOX5 is negatively correlated with Th17 cells. Some studies have shown that in colorectal cancer, patients with high infiltration of Th17 and high expression of its immune marker RORC have a poor prognosis [42]. This seems to be one of the reasons why ALOX5 is highly expressed in CHOL, but the prognosis of patients is better ( Figure 5E).
In addition, the oncogene STMN1 is positively correlated with Th2 cells in CHOL. Studies have shown that patients with high tumor malignancy have a low Th1/Th2 ratio; generally speaking, tumor cells produce immunosuppressive cytokines and prostaglandins, which relatively enhance the immune effect of Th2 cells and significantly reduce the antitumor ability [42][43][44]. Interestingly, we know that BIRC5, CDC25C, HELLS, NEIL3, and CENPF are common immune markers of Th2 cells. Some scholars believe that BIRC5 plays an important role in cell division and can inhibit apoptosis [45]. In hepatocellular carcinoma, the upregulation of BIRC5 expression promotes the proliferation of tumor cells [46]. CDC25C plays a significant role in regulating the cell cycle and is highly expressed in prostate cancer [47]. The high expression of HELLS and NEIL3 is related to the shortening of survival time associated with liver cancer and other tumors [48,49]. The upregulated expression of CENPF is closely related to the poor prognosis and progression of liver cancer, lung adenocarcinoma and other tumors [50,51]. These findings seem to explain the poor prognosis of patients with high expression of STMN1 in CHOL ( Figure 5F).
Our study has some limitations. First, the amount of patient information available in public databases is insufficient, and clinical staging, treatment options and other valuable data could not be obtained, therefore, a large amount of clinicopathological information needs to be collected later for prospective studies. Moreover, the mechanism by which the DEFRGs regulated CHOL in this prognostic model needs to be explored in further experimental studies.

Conclusions
In conclusion, this study is the first to construct a polygenic prognostic model of ferroptosis in CHOL, showing that ALOX5, BNIP3, EGFR and STMN1 act together to strongly contribute to CHOL cell apoptosis, ferroptosis and the immune microenvironment and confirming the diagnostic and prognostic value of this model in CHOL. As relevant studies are scarce, this topic also provides ideas and directions for research related to ferroptosis in CHOL, and provides a reference for finding new treatment modalities.