ORIGINAL RESEARCH article
Identification and Validation of a Prognostic Signature for Prostate Cancer Based on Ferroptosis-Related Genes
- 1Department of Urology, Shanghai Tenth People’s Hospital, School of Medicine in Tongji University, Shanghai, China
- 2Department of Urology, The Second Hospital of Hebei Medical University, Shijiazhuang, China
- 3Department of Orthopedics, Jingan District Zhabei Central Hospital, Shanghai, China
Ferroptosis, an iron-dependent form of selective cell death, is involved in the development of many cancers. However, ferroptosis related genes (FRGs) in prostate cancer (PCa) are not been well studied. In this study, we collected the mRNA expression profiles and clinical information of PCa patients from TCGA and MSKCC databases. The univariate, LASSO, and multivariate Cox regression analyses were performed to construct a prognostic signature. Seven FRGs, AKR1C3, ALOXE3, ATP5MC3, CARS1, MT1G, PTGS2, and TFRC, were included to establish a risk model, which was validated in the MSKCC dataset. The results showed that the high-risk group was apparently correlated with copy number alteration load, tumor burden mutation, immune cell infiltration, mRNAsi, immunotherapy, and bicalutamide response. Moreover, we found that TFRC overexpression induced the proliferation and invasion of PCa cell lines in vitro. These results demonstrate that this risk model can accurately predict prognosis, suggesting that FRGs are promising prognostic biomarkers and potential drug targets in PCa patients.
Prostate cancer (PCa) is one of the world’s most common malignancies in men, and is the second-leading cause of cancer-related deaths in Western countries (1). Prostate-specific antigen (PSA) has been used as the standard PCa detection test since the 1990s. However, previous studies found no significant differences in mortality between PSA-screened patients and those who are not screened (2–4). PSA is also known to be a major predictor of PCa prognosis. Nearly 27-53% of patients who have undergone radical prostatectomy and radiotherapy experience PSA recurrence (5). Such biochemical recurrence (BCR) can contribute to the development of advanced castration-resistant prostate cancer (CRPC) stage, leading to an increased risk in distant metastases, prostate cancer-specific mortality and overall mortality (6, 7). Therefore, it is of great significance to identify novel prognostic biomarkers for PCa.
Ferroptosis, a newly identified form of regulated cell death characterized by iron accumulation and lipid peroxidation, is distinct from other forms of regulated cell death (necroptosis, apoptosis, or autophagic cell death) (8). Recently, emerging evidence suggested that ferroptosis is related to cancer initiation, progression, or drug sensitivity (9–11). Ferroptosis inducers can be used to resolve drug resistance and prevent cancer progression or metastasis. For example, Sun et al. reported that metallothionein-1G (MT-1G) knockdown enhanced sorafenib sensitivity in hepatocellular carcinoma via promoting ferroptosis (10). It has been shown that suppression of cysteine dioxygenase 1 (CDO1) increases cellular glutathione (GSH) levels, inhibits reactive oxygen species (ROS) generation, and decreases lipid peroxidation in erastin-treated gastric cancer cells (12).
The role of ferroptosis in PCa has drawn attention in recent years. Butler et al’s study revealed that DECR1 knockdown in PCa cells inhibited tumor cell proliferation and migration via cellular polyunsaturated fatty acids (PUFAs) accumulation, thus enhancing mitochondrial oxidative stress and lipid peroxidation, and promoting ferroptosis (13). Pannexin 2 (PANX2) has been proposed to be a new marker in PCa, that promotes proliferation and invasion of PCa cells through regulating ferroptosis (14). While preliminary evidence has identified several markers that are correlated with PCa ferroptosis, the association between other ferroptosis-related genes (FRGs) and PCa prognosis remains largely unknown.
In our research, we analyzed the mRNA expression profiles of 40 ferroptosis-related genes and clinical data of PCa patients from the TCGA database. We then evaluated their differential expression in different risk PCa samples and investigated the enriched pathways and biological roles. Next, a FRGs-based prognostic model were constructed in the TCGA dataset and validated it in another dataset. Moreover, we chose the hub gene transferrin receptor (TFRC) for further validation in vitro.
Materials and Methods
The RNA-seq (FPKM value) data of 499 PCa and 52 normal prostate tissues with related clinical data were obtained from the TCGA website (https://portal.gdc.cancer.gov/repository). The MSKCC data set downloaded from the GEO dataset (https://www.ncbi.nlm.nih/geo/query/), containing genomic profile of 218 prostate tumors was used as the validation cohort
Gene Signature Building
A total of 40 genes from Stockwell et al’s research were analyzed in the current study (15) (shown in Supplementary Table S1). Univariate Cox was performed to select genes related to the recurrence-free survival (RFS), using P values less than 0.05 as the cut-off. The LASSO method was used to minimize the risk of overfitting. Finally, multiple stepwise Cox regression was used to establish the risk model. The risk score of hub genes was established as (exprgene1 × coefficientgene1) + (exprgene2 × coefficientgene2) + ⋯ + (exprgene7 × coefficientgene7).
Clustering, Genetic Alterations, and Functional Enrichment Analysis
Consensus clustering based on ferroptosis genes was achieved using the “Consensus ClusterPlus” package (16). The “clusterProfiler” R package (17) was used to perform GO and KEGG analyses based on differentially expressed genes (DEGs) (|log2FC| ≥ 1, FDR < 0.05) between different risk groups. The infiltrating score of 28 immune cells was calculated with ssGSEA in the “GSVA” R package (18). The annotated gene set file is shown in Supplementary Table S2. Genetic alterations of selective genes in TCGA patients were acquired from the cBioPortal (19, 20). The relationships between mRNA expression level and Gleason score, lymph nodal metastasis, methylation levels were obtained from UALCAN (21).
Copy Number Variation (CNV) Load, Tumor Mutation Burden (TMB), Neoantigens, Tumor Stemness, and Clonal Score
The TCGA data was divided into two risk groups based on the median risk value. The CNV load at the focal and arm levels were then calculated based on the GISTIC_2.0 results, which were freely available from Broad Firehose (https://gdac.broadinstitute.org/). Each patient’s tumor mutation burden (TMB) was measured as the total number of non-synonymous mutations per megabase. Tumor neoantigens, which can be recognized by neoantigen-specific T cell receptors (TCRs) and participate in T-cell-mediated antitumor immune response, were analyzed using the TCIA (https://tcia.at/). We also analyzed the clonality of PCa patients from the TCIA dataset (Clonality is the critical character of tumors). Considering that tumor stemness is also one of the basic characteristics of the tumor, the mRNAsi data was analyzed as reported by Robertson (22).
Prediction of Immunotherapy and Chemotherapy Responses
The Tumor Immune Dysfunction and Exclusion (TIDE) and subclass mapping (23–25) were performed to assess the clinical response of immune therapy in PCa patients. Chemotherapeutic response to three commonly used drugs for PCa patients in the TCGA and MSKCC datasets was calculated using the Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org) via performing ‘pRRophetic’ R package (26).
Cell Culture, TFRC Stably Over-Expressed Cell Lines Establishment and Reagents
The HEK293T cell line, and prostate cancer cells C4-2, PC-3, LNCaP, 22RV1 were purchased from American Type Culture Collection (ATCC, Manassas, VA). HEK293T cells were cultured in a humidified atmosphere (37°C with 5% CO2) in DMEM whereas prostate cancer cells were grown in RPMI1640. 22RV1 and LNCaP cells were treated with 5-Azacytidine (Sigma, #A2385). The TFRC overexpressing plasmid was designed and synthesized by Genomeditech (Shanghai, China). The TFRC overexpression plasmids, psPAX2 and pMD2.G were mixed and added into HEK293T cells to package lentivirus. After incubation for 48 h, lentivirus supernatants were collected, filtered, and used to infect cells.
Western Blot Analysis
The cells were washed thrice with cold PBS buffer and lysed with RIPA buffer on ice. We loaded the same amounts (40 µg) of different protein samples onto the SDS-PAGE gel and then transferred the protein to PVDF membranes. The membranes were blocked using skim milk (5%) for 1 h and incubated with primary antibodies against TFRC (ABclonal, #A5865) or β-Actin (sc-517582) at 4°C overnight. The next day, we incubated the PVDF membranes with HRP-conjugated secondary antibodies (mouse or rabbit) at room temperature for 1h. The blots were visualized using an ECL system (Thermo Fisher Scientific) after washing thrice using TBST.
Cell Invasion Assay
We used transwell chambers (8 μm pore size, Corning, MA, USA) to detect the migration and invasion of 22RV1 and LNCaP after TFRC overexpression. The upper chamber was pre-coated with Matrigel (Corning, USA) for invasion and was left uncoated for migration. Next, 105 cells were seeded in serum-free media in the upper chamber, while 600 μL of culture media containing was added in the lower chamber. After 12 h (migration) or 24 h (invasion), we used a cotton swab to remove the cells remaining in the upper surfaces of the chambers and fixed the cells on the lower surface in 100% methanol for 15 min. We then stained them using 0.1% crystal violet solution for 20 min. The cells were counted and averaged across five randomly chosen fields using a microscope.
All data analyses were performed using the R platform (v.4.0.2, https://cran.r-project.org/) or GraphPad Prism 8.0. The comparison of mRNA expression between PCa tissues and adjacent nontumorous samples, CNV, TMB, mRNAsi, clonal score, chemotherapy response and immunotherapy response across different risk groups were calculated using the Wilcoxon rank-sum test. The RFS of different groups was measured using Kaplan-Meier log-rank test method. The DEGs between different risk groups (high and low) were determined utilizing the “limma” R package (27). The results of MTT, invasion and migration results were analyzed using Student’s t test. P < 0.05 was considered as statistically significant.
Expression and Correlation of FRGs in the TCGA Cohort
Most FRGs (32/40, 80%) (Figures 1A, C), (24/40,60%) (Figures 1B, E) in the TCGA dataset were abnormally expressed in tumor samples compared with adjacent nontumorous tissues, and paired nontumorous tissues. The interaction network among these 40 genes demonstrated that TP53, GCLC, GCLM, GPX4 and other genes had more interactions (Figure 1D). Moreover, the correlations between these genes were very high (Figure 1F), for example, ACSL4 and NFE2L2, HMGCR and NCOA4, HSPB1 and CRYAB.
Figure 1 The landscape of ferroptosis-related genes (FRG) in prostate cancer. (A) Heatmap of 40 FRGs between 499 prostate cancer and normal tissues. (B) Heatmap of 40 FRGs between 52 tumor and adjacent normal pairs. (C) mRNA expression levels FRGs levels in total prostate cancer and normal tissues. (D) The PPI network acquired from the STRING database among the FRGs. (E) mRNA expression levels FRGs levels in prostate cancer and paired tissues. (F) The correlation among FRG in prostate cancer. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001. ns, not significant.
Clustering, Construction and Validation of FRGs Prognostic Model
The TCGA cohort could be divided into two clusters according to the mRNA expression levels of the 40 FRGs (Figure 2A). Patients in cluster 1 had significantly poorer RFS than those in cluster 2 (Figure 2B), suggesting that FRGs were associated with differences between PCa patients. Therefore, identifying prognostic FRGs was necessary. The univariate Cox regression method was first used to identify FRGs that were associated with prognosis; 17 genes were found to be correlated with PCa RFS (Table 1). The LASSO method was then used for selecting the hub genes (Figure 2C) to minimize the risk of overfitting. To further identify the FRGs with the greatest prognostic value, we conducted multiple stepwise Cox regression and chose seven hub FRGs to construct the prognostic model of PCa patients (Figure 2D). These seven hub genes had greater than 1% genetic alterations, for example TFRC (4%) and ALOXE3 (4%) (Figure 2E). The risk scores of each PCa patient were measured using the following method: Risk score = (0.213*ExpAKR1C3) + (0.224*ExpALOXE3) + (0.183*ExpATP5MC3) + (0.182*ExpCARS1) + (-0.346*ExpMT1G) + (-0.193*ExpPTGS2) + (0.299*ExpTFRC). We then divided the PCa patients of the TCGA cohorts into two groups (high and low risk) according to the median value of the risk score.
Figure 2 RFS of PRAD patients in cluster 1/2 subgroups and risk signature with 7 FRGs. (A) Consensus Clustering matrix for k =2. (B) Kaplan-Meier curves of DFS. (C) The Cross-Validation fit curve calculated by lasso regression method. (D) The coefficients of seven hub FRGs estimated by multivariate Cox regression. (E) Genetic alterations of seven hub FRGs. All the survival analysis are based on the recurrence free survival data.
The K-M plot demonstrated that the high-risk group had unfavorable RFS compared with the low-risk group (P < 0.0001, Figure 3A). Moreover, the area under the receiver operating characteristic curve (AUC) for 1-year, 3-year and 5-year RFS were 0.741, 0.729 and 0.736 (Figure 3B), respectively, which were good classification results. The RFS status of PRAD showed that more patients relapsed in the high-risk group (Figure 3C). The heatmap and bar plot showed that AKR1C3, ALOXE3, ATP5MC3, CARS1, TFRC in the high-risk group were increased, while MT1G and PTGS2 were decreased (Figures 3D, E). Furthermore, univariate and multivariate Cox regression analysis showed that the risk score was an independent factor RFS for PCa patients in the TCGA cohort (HR = 1.11, 95% CI:1.08-1.15, p value < 0.001, Table 2). In the MSKCC dataset, the KM plot results showed that the high risk group had unfavorable RFS, consistent with our observation in the TCGA dataset (HR = 1.76, 95% CI: 1.43-2.17, p < 0.001, Figure 4A and Table 2). While the associated survival statuses and mRNA expression of these seven genes were somewhat different between risk groups (Figures 4B–E), the AUC for the 1-year, 3-year, and 5-year RFS, was greater than 0.7, suggesting that the risk model can be useful in predicting prognosis. Moreover, the risk score of the AUC in the two datasets was marginally higher than that of T status (TCGA:0.731>0.648; MSKCC:0.736>0.468; Figure S1F), but this not the case for the Gleason score.
Figure 3 Risk score based on seven hub FRGs in TCGA PRAD cohort. (A) Survival analysis according to risk score. (B) ROC analysis. (C) Survival status of patients. (D, E) Heatmap and barplot of the seven hub genes between high and low risk group. All the survival analysis are based on the recurrence free survival data. ****P < 0.0001.
Table 2 Univariate Cox regression and multivariate Cox regression analysis of riskScore and clinical information in TCGA and MSKCC cohorts.
Figure 4 Risk score based on seven hub FRGs in MSKCC PRAD cohort. (A) Survival analysis according to risk score. (B) ROC analysis. (C) Survival status of patients. (D, E) Heatmap and barplot of the seven hub genes between high and low risk group. All the survival analysis are based on the recurrence free survival data. *P < 0.05, ***P < 0.001, ****P < 0.0001. ns, not significant.
Functional Analyses in TCGA and MSKCC Cohorts
The DEGs between the high risk and low risk groups were analyzed in the TCGA-PCa cohort to further clarify the biological functions and pathways of FRGs. The set comprised 664 genes (adjusted P < 0.05, |log2FC| > 1), including 419 upregulated and 245 downregulated genes (Figure 5A). GO and KEGG enrichment was then performed using the ClusterProfiler R package. Interestingly, most of the GO terms were enriched for immune-related functions (Figure 5B), such as humoral immune response, regulation of humoral immune response, B cell mediated immunity, immunoglobulin complex, antigen binding and growth factor activity. Meanwhile, the KEGG terms were closely associated with several metabolism processes (Figure 5C), such as ascorbate and aldarate metabolism, steroid hormone biosynthesis, and especially drug metabolism and xenobiotics metabolism by cytochrome P450, which are essential pathways for the metabolism of many drugs. The statuses of the 28 immune infiltrating cells were then calculated using the ssGSEA method since the risk score was strongly associated with the immune response. The activated CD4 T cell, CD 56dim natural killer cell, mast cell, memory B cell, neutrophil, regulatory T cell and type 17 T helper cell were significantly different in both the low- and high-risk groups in the TCGA dataset (P <0.05, Figure 5D). The CD 56dim natural killer cell, central memory CD8 T cell, natural killer cell, natural killer T cell, and type 17 T helper cell were significantly different between groups with low or high risk in the MSKCC dataset (P <0.05, Figure 5D). The risk scores were associated with regulatory T cell, type 17 T helper cell, CD 56 bright natural killer cells and neutrophils (Figures 5E, G), indicating that FRGs can regulate the progress of PCa via the immune process pathway.
Figure 5 Potential biological pathways affected by FRGs. (A) The different expression genes (DEGs) between the high-risk and low-risk groups. (B) The Gene Ontology (GO) enrichment of DEGs. (C) The KEGG enrichment of high and low risk groups. Comparison and correlations of the ssGSEA scores between different risk groups and risk scores in TCGA (D, E) and MSKCC (F, G) dataset. *P < 0.05, **P < 0.01, ***P < 0.001. ns, not significant.
The Distinctions of Gene TMB, CNV, Cancer Stemness Index, and Sensitivity to Immunotherapy/Chemotherapy in the High- and Low-Risk Groups
The GO and KEGG enrichment results showed that the risk score was closely associated with the immune process and drug metabolism pathways. The CNV, TMB, cancer stemness index and clonal score were included to determine whether the risk score was correlated with the immune response and chemotherapy. The high-risk group of the CNV status in the TCGA cohort had both high amplification (Pamp = 3.3e-11, Pamp = 4.4e-14) and deletion (Pdel = 2.3e-09, Pdel = 6.7e-13) in the arm and focal levels (Figures 6A, B). Moreover, the high-risk group in TCGA dataset had high TMB (P = 1.1e-11 Figure 6C), neoantigens burden (P = 1.1e-05, Figure 6D), mRNAsi (P < 0.001, Figure 6E), and clonal score (P < 0.001, Figure 6F). Notably, TMB and neoantigens were potential immunotherapy response biomarkers (28). Subsequently, the immunotherapy response and predictive responses of the three common chemotherapy drugs, cisplatin, docetaxel and bicalutamide were analyzed using the TIDE and GDSC database. The results revealed that the TIDE score was higher in the high-risk group (P = 0.032, TCGA dataset, Figure 6G), while there was no significant difference in the MSKCC dataset (P = 0.81, Figure 6I). Furthermore, patients in the high-risk group appeared to benefit more from anti-PD-1/PD-L1 immunotherapy than did the low-risk group (P = 0.0299, P = 0.2657, Figure S1E). The high-risk group had a high estimated IC50 for bicalutamide (P = 8.5e-11 TCGA; P = 1.5e-05 MSKCC) and low estimated IC50 for docetaxel in TCGA (P = 0.025), and a high estimated IC50 for docetaxel in the MSKCC (P = 0.076) (Figures 6H, J).
Figure 6 The correlations of different risk group with copy number alterations, tumor mutation burden, neo-antigens, tumor stemness, clonal status and immune-/chemotherapy response. (A) Arm-level copy number amplification and deletion. (B) Focal-level copy number amplification and deletion; (C) Tumor mutant burden difference. (D) Neo-antigens. (E) Tumor stemness difference represented by the mRNAsi. (F) Clonal status. (G, H) Immunotherapy response based on TIDE website and estimated IC50 indicates the efficiency of chemotherapy in TCGA. (I, J) Immunotherapy response and estimated IC50 in MSKCC dataset. All the survival analysis are based on the recurrence free survival data.
TFRC Overexpression Facilitates Proliferation, Migration and Invasion in PCa Cell Lines
The hub genes could have more effect on the biological functions since the risk model was strongly associated with the PCa RFS. Therefore, we chose the high genetic alteration gene TFRC to confirm our hypothesis. We first measured the baseline protein levels of TFRC in PCa cells PC-3, LNCaP, 22RV1, and C4-2 using western blots. TFRC expression was relatively lower in LNCaP and 22RV1 cells compared to other cell lines (Figure 7A). Therefore, we used these two cell lines to overexpress TFRC in later experiments. Interestingly, the mRNA expression of TFRC decreased in PCa tissues (Figures 1C, E) compared with the normal samples, but increased in tissues with high Gleason scores and lymph nodes metastases (Supplementary S1A, 1B). Several investigators have reported that methylated CpG sites can block transcription initiation by inhibiting transcription factor binding, such as promoters and distal regulatory regions, then changing the mRNA expression of host genes (29–31). We explored the methylation levels of TFRC and found that TFRC had high methylation levels in PCa tissues (although the P value was not significant, Supplementary S1C), which were negatively associated with the mRNA expression of TFRC (Supplementary S1D). Therefore, we treated the PCa cells with 5-zacytidine [a widely-used DNA methylation inhibitor (32)] and found that the protein level of TFRC was increased after cells were incubated with 0.5 μM 5-zacytidine for 24 h (Figure 7B). We then overexpressed TFRC in 22RV1 and LNCaP cells (Figures 7C, D). The effect of TFRC on PCa cells was evaluated using MTT, colony formation and transwell assays. As shown in Figures 7E, F, upregulation of TFRC induced proliferation in PCa cells. Similarly, the transwell assays indicated that TFRC promotes the migrative and invasive activities of LNCaP and 22RV1 cell lines (Figure 7G).
Figure 7 TFRC expression influenced by 5-zacytidine and overexpressing TFRC promotes proliferation, migration and invasion in Pca cells 22RV1 and LNCaP. (A) The protein levels of TFRC in Pca cells (C4-2, LNCaP, 22RV1 and PC-3) were detected by western blots. (B) 5-azacytidine inhibited DNA methylation and increased the expression of TFRC in LNCaP and 22RV1 cells. (C, D). The efficiency of TFRC overexpressed lentivirus were confirmed using fluorescence microscope and western blots. (E, F). The MTT and colony formation assays indicated that TFRC promotes cells proliferation in LNCaP and 22RV1. (G, H). Transwell assays suggested that TFRC promotes cells migration and invasion in LNCaP and 22RV1. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.
Selective induction of cancer cell death is the most effective therapeutic method for tumors (33). Studies have reported that ferroptosis, a common selective induction cell death, plays pivotal role in tumorigenesis (34, 35). However, its importance in prostate cancer has not yet to be elucidated. In this study, we collected the RNA-seq data to investigate variations in mRNA expression profiles of FRGs in prostate cancer. We used univariate Cox, lasso regression, and multivariate Cox analyses to identify the signature of seven RFGs.
It has been widely reported that the seven hub genes, AKR1C3, ALOXE3, ATP5MC3, CARS1, MT1G, PTGS2, and TFRC are involved in the pathogenesis of several diseases. AKR1C3, a crucial androgenic enzyme, can reprogram AR signaling in advanced prostate cancer (36) and is involved in the production of aromatase substrates in breast cancer (37). ALOXE3, epidermal LOX type 3, converts fatty acid substrates to specific epoxyalcohol derivatives using R-hydroperoxides (38), is involved in late epidermal differentiation (39) and ichthyosis (40). ATP5MC3, also referred to as ATP5G3 (41), encodes a mitochondrial ATP synthase subunit, and is associated with overall survival in clear cell renal carcinoma patients (42). CARS1, cysteinyl-TRNA synthetase 1, is associated with tRNA function and contributes to the development of inflammatory myofibroblastic tumors (43) and kidney cancer (44). CARS1 knockdown has been proven to suppress ferroptosis induced by cysteine deprivation and promote the transsulfuration pathway (45). MT1G, a small-molecular weight protein that has high affinity for zinc ions, can inhibit proliferation by increasing the stability and transcriptional activity of p53 (46). PTGS2, also known as cyclooxygenase 2, is associated with prostanoid biosynthesis (47). Several studies have demonstrated that PTGS2 induces cancer stem cell-like activity and promotes the proliferation, angiogenesis and metastasis of cancer cells (48–50). TFRC is a cell surface receptor necessary for cellular iron (51). Few of these genes have been investigated in PCa (except PTGS2 and MT1G). In this study, all seven genes were shown to correlate with RFS in patients with PCa.
The ferroptosis-related risk model was constructed in the TCGA PCa cohort using the mRNA expression profiles of these seven hub genes. The KM plot and ROC curve showed that this risk score could easily predict long versus short RFS. These findings were also validated in the MSKCC external validation dataset. The PCa patients were divided into two groups according to the risk score to further explore the detailed FRG mechanisms. There were 664 DEGs in these two groups. Interestingly, these DEGs were enriched in immune-related GO terms and metabolism KEGG terms (Figures 5B, C). Subsequently, the immune infiltration status of 28 immune cells was calculated using the ssGSEA method, and the results showed that several immune cells, such as central memory CD8 T cell, CD56bright and natural killer cell were associated with the risk score (Figures 5D–G), suggesting that this risk model is appropriate for immunotherapy. The immunotherapy-related signatures, such as CNV load, TMB, mRNAsi, clonal score, and TIDE score were included for further investigation. These results indicated that the high-risk group was correlated with the CNV load, TMB, neoantigens, mRNAsi, and clonal score (Figures 6A–F), implying that the high-risk score group has a better clinical response to immunotherapy. The TIDE score in the high-risk group was also increased, indicating better outcomes for immune therapy (Figures 6G, I). Since most KEGG terms between the low and high-risk groups were enriched in metabolism (especially in drug metabolism) (Figure 5C), we determined the estimated IC50 for three commonly used drugs. The results showed that there was an increase of bicalutamide in the high-risk group (Figures 6H–J), implying that this group could be bicalutamide resistant. In other words, FRGs could have effects on RFS via regulating drug responses. And the high-risk patients can improve their situations through different immunotherapies, providing a potential strategy for the individual treatment of PCa patients.
Finally, we chose TFRC, which had a high genetic alteration rate, for further validation. Since the mRNA level of TFRC was lower in PCa patients compared with normal tissues, we overexpressed TFRC to determine its effects in PCa cells. LNCaP and 22Rv1 cell lines were selected based on the baseline protein expression levels (Figures 7A, C, D). MTT, clonal formation and invasion assays showed that TFRC overexpression can significantly increase the proliferation (Figure 7F) and invasion (Figure 7G) of PCa cells. We also added 5-Aza to inhibit DNA methylation of TFRC, which demonstrated that the hypermethylation status of the TFRC gene induces a relatively low PRAD tissue expressions, but increases the trend of the Gleason score. This finding suggests that TFRC can play an oncogenic role in PCa, but the detailed mechanism needs to be further explored.
In summary, our study systematically evaluated the expression of FRGs and their potential prognostic value in PCa. Moreover, a risk model of seven FRGs was established in TCGA and validated using the MSKCC dataset. We found that the high-risk group was correlated with elevated CNV load, TMB, neoantigens, mRNAsi, clonal score, and immunotherapy response, but also with a high estimated IC50 for bicalutamide. Finally, we assessed the effects of overexpressed TFRC on proliferation and invasion in PCa cell lines.
This study has several limitations. First, our risk model was established and validated based on public databases, more prospective real-world data are needed to confirm its clinical significance. Second, we only initially explored the effect of TFRC overexpression on biological functions, the specific mechanism in vivo and in vitro should be investigated for further.
We established a ferroptosis-related risk model, that was strongly associated with aberrant CNV load, TMB, mRNAsi, neoantigens, and clonal score. In addition, we validated TFRC to be an oncogenic factor in prostate cell lines. Our research provides new insights into personalized therapies for PCa patients.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
YX and LG: conceptualization, supervision, methodology, and writing -review & editing. HL and LG: data curation, formal analysis, and writing-original draft. JL, T-sZ, and TX: validation, visualization, and investigation. All authors contributed to the article and approved the submitted version.
This work was supported by the National Natural Science Foundation of China (grant NO.81971371).
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.
We thank all the people who contributed to this study.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2021.623313/full#supplementary-material
Supplementary Figure 1 | The associations between TFRC mRNA expression levels and Gleason score (A), nodal metastasis (B) and methylation levels (C, D) in TCGA PCa patients. Subclass mapping (SubMap) analysis of patients respond to anti–PD-1/PD-L1 treatment (E). Multiple ROC (Receiver Operating Characteristic) in TCGA and MSKCC (F). All the survival analysis are based on the recurrence free survival data.
1. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global Cancer Statistics 2018: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin (2018) 68:394–424. doi: 10.3322/caac.21492
2. Andriole GL, Crawford ED, Grubb RL 3rd, Buys SS, Chia D, Church TR, et al. Prostate Cancer Screening in the Randomized Prostate, Lung, Colorectal, and Ovarian Cancer Screening Trial: Mortality Results After 13 Years of Follow-Up. J Natl Cancer Inst (2012) 104:125–32. doi: 10.1016/j.yuro.2012.06.019
3. Martin RM, Donovan JL, Turner EL, Metcalfe C, Young GJ, Walsh EI, et al. Effect of a Low-Intensity PSA-Based Screening Intervention on Prostate Cancer Mortality: The CAP Randomized Clinical Trial. JAMA (2018) 319:883–95. doi: 10.1001/jama.2018.0154
4. Schroder FH, Hugosson J, Roobol MJ, Tammela TL, Zappa M, Nelen V, et al. Screening and Prostate Cancer Mortality: Results of the European Randomised Study of Screening for Prostate Cancer (ERSPC) at 13 Years of Follow-Up. Lancet (2014) 384:2027–35. doi: 10.1016/S0140-6736(14)60525-0
5. Yuan P, Ling L, Fan Q, Gao X, Sun T, Miao J, et al. A Four-Gene Signature Associated With Clinical Features Can Better Predict Prognosis in Prostate Cancer. Cancer Med (2020) 9(21):8202–15. doi: 10.1002/cam4.3453
6. Van den Broeck T, van den Bergh RCN, Arfi N, Gross T, Moris L, Briers E, et al. Prognostic Value of Biochemical Recurrence Following Treatment With Curative Intent for Prostate Cancer: A Systematic Review. Eur Urol (2019) 75:967–87. doi: 10.1016/j.eururo.2018.10.011
7. Meng J, Lu X, Zhou Y, Zhang M, Gao L, Gao S, et al. Characterization of the Prognostic Values and Response to Immunotherapy/Chemotherapy of Kruppel-Like Factors in Prostate Cancer. J Cell Mol Med (2020) 24:5797–810. doi: 10.1111/jcmm.15242
8. Kim EH, Shin D, Lee J, Jung AR, Roh JL. CISD2 Inhibition Overcomes Resistance to Sulfasalazine-Induced Ferroptotic Cell Death in Head and Neck Cancer. Cancer Lett (2018) 432:180–90. doi: 10.1016/j.canlet.2018.06.018
12. Hao S, Yu J, He W, Huang Q, Zhao Y, Liang B, et al. Cysteine Dioxygenase 1 Mediates Erastin-Induced Ferroptosis in Human Gastric Cancer Cells. Neoplasia (2017) 19:1022–32. doi: 10.1016/j.neo.2017.10.005
13. Nassar ZD, Mah CY, Dehairs J, Burvenich IJ, Irani S, Centenera MM, et al. Human DECR1 Is an Androgen-Repressed Survival Factor That Regulates PUFA Oxidation to Protect Prostate Tumor Cells From Ferroptosis. Elife (2020) 9:e54166. doi: 10.7554/eLife.54166
14. Liao D, Yang G, Yang Y, Tang X, Huang H, Shao J, et al. Identification of Pannexin 2 as a Novel Marker Correlating With Ferroptosis and Malignant Phenotypes of Prostate Cancer Cells. Onco Targets Ther (2020) 13:4411–21. doi: 10.2147/OTT.S249752
15. Stockwell BR, Friedmann Angeli JP, Bayir H, Bush AI, Conrad M, Dixon SJ, et al. Ferroptosis: A Regulated Cell Death Nexus Linking Metabolism. Redox Biol Disease Cell (2017) 171:273–85. doi: 10.1016/j.cell.2017.09.021
19. Cerami E, Gao J, Dogrusoz U, Gross BE, Sumer SO, Aksoy BA, et al. The cBio Cancer Genomics Portal: An Open Platform for Exploring Multidimensional Cancer Genomics Data. Cancer Discov (2012) 2:401–4. doi: 10.1158/2159-8290.CD-12-0095
20. Gao J, Aksoy BA, Dogrusoz U, Dresdner G, Gross B, Sumer SO, et al. Integrative Analysis of Complex Cancer Genomics and Clinical Profiles Using the Cbioportal. Sci Signal (2013) 6:pl1. doi: 10.1126/scisignal.2004088
21. Chandrashekar DS, Bashel B, Balasubramanya SAH, Creighton CJ, Ponce-Rodriguez I, Chakravarthi B, et al. Ualcan: A Portal for Facilitating Tumor Subgroup Gene Expression and Survival Analyses. Neoplasia (2017) 19:649–58. doi: 10.1016/j.neo.2017.05.002
22. Robertson AG, Kim J, Al-Ahmadie H, Bellmunt J, Guo G, Cherniack AD, et al. Comprehensive Molecular Characterization of Muscle-Invasive Bladder Cancer. Cell (2017) 171:540–56.e525. doi: 10.1016/j.cell.2017.09.007
25. Lu X, Jiang L, Zhang L, Zhu Y, Hu W, Wang J, et al. Immune Signature-Based Subtypes of Cervical Squamous Cell Carcinoma Tightly Associated With Human Papillomavirus Type 16 Expression. Mol Features Clin Outcome Neoplasia (2019) 21:591–601. doi: 10.1016/j.neo.2019.04.003
26. Geeleher P, Cox N, Huang RS. Prrophetic: An R Package for Prediction of Clinical Chemotherapeutic Response From Tumor Gene Expression Levels. PloS One (2014) 9:e107468. doi: 10.1371/journal.pone.0107468
27. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma Powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res (2015) 43:e47. doi: 10.1093/nar/gkv007
28. Fancello L, Gandini S, Pelicci PG, Mazzarella L. Tumor Mutational Burden Quantification From Targeted Gene Panels: Major Advancements and Challenges. J Immunother Cancer (2019) 7:183. doi: 10.1186/s40425-019-0647-4
32. Christman JK. 5-Azacytidine and 5-Aza-2’-Deoxycytidine as Inhibitors of DNA Methylation: Mechanistic Studies and Their Implications for Cancer Therapy. Oncogene (2002) 21:5483–95. doi: 10.1038/sj.onc.1205699
33. Liu HJ, Hu HM, Li GZ, Zhang Y, Wu F, Liu X, et al. Ferroptosis-Related Gene Signature Predicts Glioma Cell Death and Glioma Patient Progression. Front Cell Dev Biol (2020) 8:538. doi: 10.3389/fcell.2020.00538
36. Liu C, Yang JC, Armstrong CM, Lou W, Liu L, Qiu X, et al. AKR1C3 Promotes Ar-V7 Protein Stabilization and Confers Resistance to AR-Targeted Therapies in Advanced Prostate Cancer. Mol Cancer Ther (2019) 18:1875–86. doi: 10.1158/1535-7163.MCT-18-1322
37. Penning TM. AKR1C3 (Type 5 17beta-Hydroxysteroid Dehydrogenase/Prostaglandin F Synthase): Roles in Malignancy and Endocrine Disorders. Mol Cell Endocrinol (2019) 489:82–91. doi: 10.1016/j.mce.2018.07.002
38. Krieg P, Rosenberger S, de Juanes S, Latzko S, Hou J, Dick A, et al. Aloxe3 Knockout Mice Reveal a Function of Epidermal Lipoxygenase-3 as Hepoxilin Synthase and Its Pivotal Role in Barrier Formation. J Invest Dermatol (2013) 133:172–80. doi: 10.1038/jid.2012.250
40. Eckl KM, de Juanes S, Kurtenbach J, Natebus M, Lugassy J, Oji V, et al. Molecular Analysis of 250 Patients With Autosomal Recessive Congenital Ichthyosis: Evidence for Mutation Hotspots in ALOXE3 and Allelic Heterogeneity in ALOX12B. J Invest Dermatol (2009) 129:1421–8. doi: 10.1038/jid.2008.409
42. Bruggemann M, Gromes A, Poss M, Schmidt D, Klumper N, Tolkach Y, et al. Systematic Analysis of the Expression of the Mitochondrial ATP Synthase (Complex V) Subunits in Clear Cell Renal Cell Carcinoma. Transl Oncol (2017) 10:661–8. doi: 10.1016/j.tranon.2017.06.002
43. Debelenko LV, Arthur DC, Pack SD, Helman LJ, Schrump DS, Tsokos M. Identification of CARS-ALK Fusion in Primary and Metastatic Lesions of an Inflammatory Myofibroblastic Tumor. Lab Invest (2003) 83:1255–65. doi: 10.1097/01.LAB.0000088856.49388.EA
44. Wu G, Wang Q, Xu Y, Li Q, Cheng L. A New Survival Model Based on Ferroptosis-Related Genes for Prognostic Prediction in Clear Cell Renal Cell Carcinoma. Aging (Albany N Y) (2020) 12:14933–48. doi: 10.18632/aging.103553
45. Hayano M, Yang WS, Corn CK, Pagano NC, Stockwell BR. Loss of cysteinyl-tRNA Synthetase (CARS) Induces the Transsulfuration Pathway and Inhibits Ferroptosis Induced by Cystine Deprivation. Cell Death Differ (2016) 23:270–8. doi: 10.1038/cdd.2015.93
47. Bjarnason I, Scarpignato C, Holmgren E, Olszewski M, Rainsford KD, Lanas A. Mechanisms of Damage to the Gastrointestinal Tract From Nonsteroidal Anti-Inflammatory Drugs. Gastroenterology (2018) 154:500–14. doi: 10.1053/j.gastro.2017.10.049
48. Altorki NK, Keresztes RS, Port JL, Libby DM, Korst RJ, Flieder DB, et al. Celecoxib, a Selective Cyclo-Oxygenase-2 Inhibitor, Enhances the Response to Preoperative Paclitaxel and Carboplatin in Early-Stage Non-Small-Cell Lung Cancer. J Clin Oncol (2003) 21:2645–50. doi: 10.1200/JCO.2003.07.127
50. Beeghly-Fadiel A, Wilson AJ, Keene S, El Ramahi M, Xu S, Marnett LJ, et al. Differential Cyclooxygenase Expression Levels and Survival Associations in Type I and Type II Ovarian Tumors. J Ovarian Res (2018) 11:17. doi: 10.1186/s13048-018-0389-9
51. Luo T, Gao J, Lin N, Wang J. Effects of Two Kinds of Iron Nanoparticles as Reactive Oxygen Species Inducer and Scavenger on the Transcriptomic Profiles of Two Human Leukemia Cells With Different Stemness. Nanomaterials (Basel) (2020) 10(10):1951. doi: 10.3390/nano10101951
Keywords: ferroptosis, recurrence free survival, TFRC, gene signature, prostate cancer
Citation: Liu H, Gao L, Xie T, Li J, Zhai T-s and Xu Y (2021) Identification and Validation of a Prognostic Signature for Prostate Cancer Based on Ferroptosis-Related Genes. Front. Oncol. 11:623313. doi: 10.3389/fonc.2021.623313
Received: 29 October 2020; Accepted: 28 June 2021;
Published: 15 July 2021.
Edited by:Zhengfei Zhu, Fudan University, China
Reviewed by:Wei-Ming Li, Kaohsiung Medical University, Taiwan
Giuseppina Comito, Università degli Studi di Firenze, Italy
Copyright © 2021 Liu, Gao, Xie, Li, Zhai and Xu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
†These authors have contributed equally to this work