Evaluation of Ferroptosis-related Gene AKR1C1 as a Novel Biomarker Associated with the Immune Microenvironment and Prognosis in Breast Cancer

Background

Breast cancer is the most common cancer in women, followed by nonmelanoma skin cancer,¹ as well as the worldwide top three common malignancies with lung and colon cancer.² About 1.7 million people were diagnosed with breast cancer per year with one in three of them dying from it.³ Despite significant developments of endocrine therapy and anti-HER2 targeting therapy, breast cancer remains one of the top three common causes of death in most countries, with raising mortality in some districts. Thus, the discovery of novel therapeutic targets in breast cancer is urgently needed for better patient prognosis.

Ferroptosis is the latest-discovered, iron-dependent form of regulated cell death, which is characterized by accumulating reactive oxygen species and lipid peroxidation products to lethal levels.⁴ Besides, ferroptosis plays a pivotal role in normal cells and tissues, it has been increasingly recognized that ferroptosis and related genes are also involved in oncogenesis and cancer development, including breast cancer.^5,6 Studies have also shown that some aggressive malignancies, such as triple-negative breast cancer, are intrinsically sensitive to ferroptosis due to the nonapoptotic feature,⁷ which indicates ferroptosis-related therapies as a promising treatment for these malignancies. More recently, ferroptosis and related genes are proven to be a mediator in immune microenvironment and cancer immunotherapy,⁸ which provides a novel insight of pathophysiological function of ferroptosis. Aldo-keto reductases (AKRs) gene family (AKR1C1-3) has been introduced to potentially help to metabolize lipid peroxides by catalyzing the conversion of aldehydes and ketones to alcohols, therefore, inhibiting the ferroptosis.⁹ Meanwhile, two studies validated the role of AKR1C1 in ferroptosis in lung cancer and melanoma.^10,11 However, no study reported the role of AKR1C1 in breast cancer with detailed mechanisms.

In order to reveal the importance of these ferroptosis-related genes in breast cancer, this study was conducted comprehensively, analyzing the clinical and immunological role of ferroptosis-related genes in breast cancer, which could provide a better understanding of the complex functions and mechanisms of ferroptosis-related genes in breast cancer.

Methods

Data Acquisition

Datasets of breast cancer were searched in the GEO database, which included all the ferroptosis-related genes and adequate clinical information. Three transcriptome microarray datasets (GSE22358, GSE9014, and GSE8977) were retrieved from gene expression omnibus (GEO) database for identification of ferroptosis-related genes. Meanwhile, a Venn diagram was used for illustrating co-differentially expressed genes (co-DEGs). Ferroptosis-related genes were retrieved from the reports from Liang et al’s group.¹² Further validations were realized using transcriptome and clinical data of BRCA patients from TCGA database (http://cancergenome.nih.gov/), including 1097 tumor tissues and 114 normal tissues. RNA-seq data were downloaded with the TPM format and received log2 conversions. The expression of genes were compared using one-way ANOVA or Mann–Whitney U-test. GEPIA (http://gepia.cancer-pku.cn/),¹³ UALCAN (http://ualcan.path.uab.edu/analysis.html),¹⁴ and TNMplot (https://www.tnmplot.com/) were also used to access data from TCGA database.¹⁵

Correlation Analysis

Correlation analysis was realized by The LinkedOmics database (http://www.linkedomics.org/login.php).¹⁶ The LinkFinder module was used for screening the co-expressed genes with Pearson's correlation test, and the results were shown as volcano plot and heatmap. Subsequently, gene ontology biological process (GO_BP), gene ontology cellular component (GO_CC), gene ontology molecular function (GO_MF), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were analyzed by the gene set enrichment analysis (GSEA) in the LinkInterpreter module.

Immunological Analysis

Tumor IMmune Estimation Resource (TIMER) database (https://cistrome.shinyapps.io/timer/)¹⁷ and the TISIDB database (http://cis.hku.hk/TISIDB)¹⁸ were used for validation of the abundance of tumor-infiltrating immune cells. In the meantime, the estimation of stromal and immune cells in malignant tumor tissues using expression data (ESTIMATE) algorithm was applied to generate immune score, representing the infiltration degree of immune cells in tumor tissue.

Statistical Analyses

R 3.3.0 and Statistical Package for Social Sciences 23.0 for Windows (BM Corporation, Armonk, NY, USA) were used to perform statistical analyses. One-way analysis of variance (ANOVA) was used, with homogeneous variance, to analyze differences in immune cell components between normal and breast cancer tissues, while Welch’s ANOVA was applied when variance was heterogeneous. Kaplan–Meier curve was used for survival analysis from Kaplan–Meier Plotter (http://kmplot.com/analysis/),¹⁹ where the log rank test was for comparisons between groups. The receiver operating characteristic (ROC) curve was used to investigate the value of genes in predicting patient survival. The area under the curve (AUC) was used for quantifying the prediction value. Meantime, univariate and multivariate Cox regression analyses were used to confirm the associations between patient characteristics and patient prognosis. Factors with P<0.05 were included in the multivariate Cox regression with the LR method. The main tools involved in this study were summarized in Table 1.

Table 1 The Main Tools Involved in This Study

Results

Identification of Candidate Ferroptosis-Related Genes

The data of GSE22358, GSE9014, and GSE8977 from GEO database were obtained to identify differentially expressed genes in breast cancer. One thousand, three hundred and twenty-five upregulated DEGS and 1625 downregulated DEGS were found in GSE22358, while 10,650 and 7356 in GSE9014, and 927 and 775 in GSE8977. As shown in Figure 1A and B, Venn plot showed 65 upregulated co-DEGS and 77 downregulated co-DEGS. Venn plots (Figure 1C and D) between co-DEGS and ferroptosis-related genes showed that two ferroptosis-related genes (AKR1C1 and AKR1C3) were downregulated. The prognostic value of AKR1C1 and AKR1C3 was subsequently investigated using overall survival (OS) data of breast cancer patients from GSE2990, GSE2034 (Figure 1E–H). Although both the higher AKR1C1 and AKR1C3 showed a better OS survival curve in breast cancer, the differences of AKR1C3 in both GSE2990 and GSE2034 were not significant. Therefore, further investigations were focused on AKR1C1.

Figure 1 Identification of candidate ferroptosis-related genes. (A and B) Venn plots of to identify upregulated (A) and downregulated (B) co-DEGs. (C and D) Venn plots of to identify ferroptosis-related co-DEGs. (E–H). Kaplan–Meier curves of AKR1C1 and AKR1C3 on overall survival of GSE2990 and GSE2034.

Validation of the Expression Level and Clinical Significance of AKR1C1

The expression level of AKR1C1 was significantly lower in breast cancer tissues compared with normal control in GSE22358 (Figure 2A), GSE9014 (Figure 2B), and GSE8977 (Figure 2C). To validate the significance and explore further information of AKR1C1, data from the BRCA cohort of TCGA database and normal breast tissues from GTEx were also used. The results showed a consistent decrease of AKR1C1 expression level in breast cancer, regardless of unpaired tissues (Figure 2D) or paired tissues (Figure 2E). In addition, several databases, such as UALCAN (Figure 2F), TNMplot (Figure 2G), and GEPIA (Figure 2H), were enrolled to validate the downregulated levels of AKR1C1 in BRCA patients. Meanwhile, the diagnostic and prognostic significance were also investigated. The baseline data of patient characteristics were shown in Table 2. Patients were divided into two groups according to the expression of AKR1C1. The median age of both low expression group and high expression group was 58 years old (P=0.778). There was no significant difference in TNM stage between the two groups. One hundred and thirty (24.6%) and 130 (24.4%) patients suffered from stage III or stage IV in the low expression group and high expression group (P=0.732). However, the low expression group had a higher positive rate of PR, ER, and HER2 compared with high expression group (P<0.001). The ROC curve showed a strong diagnostic value of AKR1C1 with AUC of 0.948 (95%CI: 0.933–0.963) (Figure 2I), while the survival curve of OS from Kaplan–Meier Plotter confirmed the high level AKR1C1 representing good patient’s prognosis (Figure 2J). Univariate and multivariate Cox regression analyses still showed independent prognostic value of AKR1C1 (Table 3).

Table 2 The Characteristics of BRCA Patients from TCGA Database

Table 3 Univariate and Multivariate Analyses of Risk Factors with and OS in BRCA Patients

Figure 2 Validation of the expression level and clinical significance of AKR1C1. Expression comparisons of AKR1C1 between breast cancer and normal control in GSE22358 (A), GSE9014 (B), GSE8977 (C), unpaired TCGA-BRCA (D), and paired TCGA-BRCA (E). Expression of AKR1C1 in TCGA-BRCA validated by UALCAN (F), TNM plot (G), and GEPIA (H). The ROC curve of AKR1C1 on diagnosing breast cancer (I). The Kaplan–Meier curve of AKR1C1 on overall survival of TCGA-BRCA (J). *P<0.05; **P<0.01;***P<0.001.

AKR1C1 Co-expression Analyses

To explore the knowledge of biological function of AKR1C1 in breast cancer, the LinkedOmics tool was used to obtain the co-expression pattern of AKR1C1 in TCGA-BRCA. 7441 genes showed a positively correlation with AKR1C1, while 4424 genes a negatively correlation (Figure 3A). Figure 3B and C showed the top 50 genes positively and negatively co-expressed with AKR1C1, respectively. The prognostic values of these genes were subsequently investigated (Figure 3D and E). Most positively co-expressed genes (38 out of 50) had a protective hazard ratio, and most of the negatively co-expressed genes (39 of 50) had an unfavorable hazard ratio, which is consistent with the prognostic role of AKR1C1 in breast cancer.

Figure 3 The co-expression genes with AKR1C1 from the LinkedOmics database. (A) Genes significantly associated with AKR1C1 by Pearson test. (B and C) Heat maps of top 50 genes positively and negatively related to AKR1C1A. Red represents positively linked genes and blue represents negatively linked genes. (D and E) Survival map of the top 50 genes positively (D) and negatively (E) associated with AKR1C1A in TCGA-BRCA. (F–I) GO terms for biological process, cellular component, molecular functions, and KEGG pathways of AKR1C1A by GSEA analyses.

Moreover, GSEA analysis was performed to explore the potential biological functions of AKR1C1. The enrichment results of GO terms and KEGG pathway were shown in Figure 3F–I. Termination of DNA-templated transcription, methyltransferase complex and helicase activity were the top enriched terms in biological process, cellular component, and molecular functions. Fanconi anemia pathway, aminoacyl-tRNA biosynthesis and homologous recombination were the top three enriched KEGG pathways.

The Roles of AKR1C1 in the Immune Microenvironment

Due to the revealed association between ferroptosis and immune microenvironment,²⁰ we subsequently attempted to explore the role of AKR1C1 in the immune microenvironment of breast cancer. We further investigated associations between AKR1C1 with immune infiltration cells using ssGSEA algorithm. Dendritic cells (DC), neutrophils, macrophages, B cells, mast cells, CD8 T cells, NK cells, TFH, and eosinophils were found positively correlated with AKR1C1 in breast cancer (Figure 4A). The top three cell type (DC, neutrophils, macrophages) were subsequently validated using the TIMER database and the TISIDB database, and the results were consistent (Figure 4B–G). Moreover, the ESTIMATE algorithm showed a positive correlation between AKR1C1 and the infiltration immune score in breast cancer (Figure 4H). We further investigated the prognostic value of the immune score, and the results showed patients with higher immune score had a better OS (Figure 4I), indicating AKR1C1 may improve the patient prognosis by influencing immune microenvironment in breast cancer. In addition, we also observed the positive associations between AKR1C1 levels and multiple immune checkpoint molecules, including PD-1 (Figure 4J), PD-L1 (Figure 4K), CTLA4 (Figure 4L), B7-H3 (Figure 4M), and VSIR (Figure 4N).

Figure 4 The associations between immune infiltration with AKR1C1. (A) Immune infiltration cells associated with AKR1C1 using ssGESA algorithm. (B–D) Validations of top three cell (DC, neutrophils, macrophages) by the TIMER database. (E–G) Validations of top three cell (DC, neutrophils, macrophages) by the TISIDB database. (H) The association between AKR1C1 with immune score calculated by ESTIMATE algorithm. (I) The Kaplan–Meier curve of immune score on overall survival of TCGA-BRCA. (J–N) The association between AKR1C1 with PD-1, PD-L1, CTLA4, B7-H3, and VSIR.

Further investigations were focused on broadening the recognition of immune role of AKR1C1 other than immune infiltration cells. TISIDB database was used to analyze the associations between AKR1C1 with immunostimulators, immunoinhibitors, chemokines, and receptors. As shown in Figure 5, multiple immunostimulators and immunoinhibitors showed associations with expression of AKR1C1, with IL-6 as the most significant immunostimulatory marker and BTLA as the most significant immunoinhibitory marker. Meanwhile, Figure 6 showed that CXCL2 was the most significant chemokine related to AKR1C1, and CCR7 was the most significant chemokine receptor. These results together suggested that AKR1C1 might play promising roles in the prognosis of patients with breast cancer through regulating the immune response.

Figure 5 Associations of the AKR1C1 expression with immunomodulators from the TISIDB database. (A) Correlations between immunostimulators and AKR1C1 with the four most significant molecules. (B) Correlations between immunoinhibitors and AKR1C1 with the four most significant molecules.

Figure 6 Associations of the AKR1C1 expression with chemokines and receptors from the TISIDB database. (A) Correlations between chemokines and AKR1C1 with the four most significant molecules. (B) Correlations between receptors and AKR1C1 with the four most significant molecules.

Discussion

In this study, we found a ferroptosis-related gene, AKR1C1, plays a significant role in the oncogenesis and development of breast cancer, functioning as a potential tumor suppressor. Further investigations showed that AKR1C1 was significantly associated with the immune microenvironment of breast cancer, indicating that AKR1C1 might make an impact on patient prognosis through altering immune microenvironment in breast cancer. These results suggested that AKR1C1 had potential value for use as a treatment target or a biomarker of patient prognosis in breast cancer.

There are several studies introducing the significance of ferroptosis and regulators in breast cancer. One of the features in ferroptosis is the inactivation of GPX4, which has also been proven in breast cancer.²¹ Surprisingly, the drug-tolerant breast cancer cells seem to be more vulnerable to inducing ferroptosis by GPX4 inhibition.²² Alternatively, lipid reactive oxygen species (ROS) is the main killer in ferroptosis. Among the regulators involved in lipid metabolism, ACSL4 is specifically associated with ferroptosis.²³ Particularly, ACSL4 was reported to be a promoter of ferroptosis by enriching cellular membranes with long polyunsaturated n-6 fatty acids in breast cancer.²⁴ In this study, we found the expression of AKR1C1 was downregulated in tumor tissues, which might indicate the treatment of inducing ferroptosis may be potential effective in breast cancer.

As an important part of tumor microenvironment (TME), immune microenvironment composed of various immune cells is widely recognized as a biomarker of tumor progression and predictors of clinical responses to therapeutic regimens, especially immunotherapy.²⁵ DC cells are essential for antitumor immune responses due to their functions of cross-presenting tumor antigen to lymphocytes. In breast cancer, a high proportion of DC was reported in less aggressive tumors, which were explained by getting involved in T-cell dependent immune clearance, therefore influencing patient prognosis.²⁶ TME neutrophils (TANs) may have antitumor functions and pro-tumorigenic functions at the same time.²⁷ TANs can function against tumor cells by promoting inflammatory, while it also possesses the anti-inflammatory and immunosuppressive abilities, which may be regulated by TME.²⁸ Similarly, tumor-associated macrophage (TAM) also displays multi-aspect functions ranging from antitumor responses and pro-tumorigenic responses.²⁹ Only one previous study by Wenners et al referred the association between AKR1C1 expression and overall survival, however, it failed to investigate the further underlying mechanisms.³⁰ In this study, we found that AKR1C1 was associated with the immune microenvironment of breast cancer, thus, it may make an impact on patient prognosis and tumor development. As far we know, this was the first study reporting the association between AKR1C1 and the immune microenvironment of breast cancer, suggesting that AKR1C1 could be developed as a therapeutic target or an indicator of immunotherapy.

There were several limitations in this study. This work revealed the significance of AKR1C1 in breast cancer and the association between AKR1C1 and the immune microenvironment by bioinformatic analyses. We strongly recommend further experimental and clinical research to comprehensively validate the biological role of AKR1C1 in influencing the immune microenvironment and patient prognosis of breast cancer.

Conclusion

In summary, this is the first study reporting a ferroptosis-related gene, AKR1C1, which may be associated with the immune microenvironment and therefore influence development of tumor and patient prognosis in breast cancer.