This study included two independent PCa cohorts. Data of the PCa patient cohort were downloaded from The Cancer Genome Atlas (TCGA) database and used as the training set (2021.03.01), which included 483 patients with comprehensive transcriptome (FPKM standardized data) and clinical information¹. Data of the MSK-IMPACT Clinical Sequencing Cohort (MSKCC) PCa cohort were obtained from the Cbioportal database² and used as the validation set, which included 138 patients with complete expression profile (normalized log2 mRNA expression data) and clinicopathological information. In addition, several Gene Expression Omnibus (GEO) datasets were selected to verify our results, including GSE54460, GSE70769, GSE104749, GSE88808, GSE70768, GSE69223, GSE68555, GSE55945, GSE46602, GSE38241, GSE35988, GSE32571, GSE32448, GSE6919, and GSE3325³. The 84 ferroptosis-related genes were taken from a previous study (Liang et al., 2020) and the ‘‘WP_FERROPTOSIS’’ gene set, which was obtained from the Molecular Signatures Database⁴. The clinicopathological features of the two cohorts are summarized in Table 1.

The risk prognostic model was constructed based on the TCGA cohort. The “Limma” package was used to obtain differentially expressed ferroptosis-related genes between PCa tissues and normal tissues, with a false discovery rate (FDR) of <0.05 as the boundary. Univariate Cox regression analysis was employed to identify the ferroptosis-related genes associated with the BCR of PCa. The genes contained in the intersection of the above two analyses were used as the core genes to construct the protein interaction network based on the STRING database⁵ and the expression correlation network based on the expressions of these core genes. Lasso regression was employed to avoid overfitting of the final prediction model. The calculation formula of the model is as follows: Riskscore=∑m=1nCoefm×Expm, where Coef_m is the risk coefficient and Exp_m is the relative mRNA expression of each ferroptosis-related gene. The median risk score was used as a cutoff to distinguish patients with high- and low-risk PCa.

The TCGA and MSKCC cohorts were used to evaluate and verify the prognosis model separately, and the same statistical methods were applied in both cohorts. Log-rank and Kaplan–Meier (K–M) tests were used to visualize the difference in the biochemical relapse-free survival (bRFS) between the two risk groups. The sensitivity and specificity of survival prediction were tested by the receiver operating characteristic (ROC) analysis. An area under the ROC curve (AUC) was calculated as an index of the prediction accuracy, and “vegan” and “stats” packages were used for the principal component analysis (PCA) to explore the distribution among groups. Univariate and multivariate Cox regression analyses were used to investigate whether the risk score can be used as an independent prognostic factor in PCa. The correlation between the risk score and other clinicopathological features was shown by a heat map.

Hallmark and Gene Ontology (GO) gene sets were downloaded for the gene set enrichment analysis (GSEA) to determine which gene sets were significantly different between the high- and low-risk groups. A total of 1,000 gene set permutations were performed to finally obtain the normalized enrichment score, normalized p-value, and FDR. The ferroptosis potential index (FPI) was calculated according to the methods of Liu et al. (2020c) by subtracting the enrichment score of the negative- from the positive-core machine components through single-sample GSEA (ssGSEA) by using the “GSVA” package.

Single-sample GSEA was also used to quantify the immune infiltration level (Rooney et al., 2015). The annotated gene set file was derived from the study of Liang et al. (2020). Finally, the enrichment levels of 16 immune cells and 13 immune-related pathways in each PCa sample were quantified, and the results were expressed by box plots. Moreover, we predicted the TME of PCa by using the “ESTIMATE” package to calculate the immune/stromal/ESTIM/tumor purity score (Yoshihara et al., 2013).

The response of each patient with PCa in both cohorts to bicalutamide and docetaxel was obtained from the Genomics of Drug Sensitivity in Cancer⁶ using the “pRRophetic” R package (Geeleher et al., 2014). Half-maximal inhibitory concentration (IC50) was used to measure the response of tumor cells to drugs.

Human prostate and PCa cell lines (PWR1E and DU145) were obtained from the American Type Culture Collection (Manassas, VA, United States) and the Shanghai Institute of Biochemistry and Cell Biology, Chinese Academy of Sciences (Shanghai, China), which were cultured in keratinocyte serum-free medium containing 50 μg/ml bovine pituitary extract and 5 ng/ml epidermal growth factor for PWR1E and RPMI-1640 medium supplemented with 2 mM L-Glutamine, 10% fetal bovine serum, 1% penicillin/streptomycin for DU145. These cells were grown at 37°C with 5% CO₂ in a humidified incubator. The lentiviral vector carrying AIFM2 shRNA or NFS1 shRNA and the corresponding control shRNA were synthesized and purchased from GenePharma and GenScript (China), respectively.

Cell viability was assessed using the MTT assay to evaluate the effects of different concentrations of erastin on cell viability. PWR1E or DU145 cell lines at 2 × 10³ cells/well were incubated in 96-well plates for 24 h. An MTT solution (5 mg/ml) was cultured for 4 h. The cells were then treated with 150 μl of dimethylsulfoxide. Cell viability was determined by enzyme-linked immunosorbent assay at 570 nm.

Various cell biological methods were used to evaluate ferroptosis. Lipid peroxidation was assessed using the BODIPY^TM 581/591 C11 (D3861, Thermo Fisher Scientific), and oxidation of the polyunsaturated butadienyl portion of the dye results in the shifting of the fluorescence emission peak from ∼590 (red) to ∼510 nm (green). PWR1E or DU145 cell lines were incubated with 10 μM erastin for 24 h. Then 1 μg/ml of Hoechst 33342 and 1 μm of BODIPY^TM 581/591 C11 were added to the culture medium for living cell imaging at the last hour of incubation. Ferroptosis-related indexes, including the levels of Fe²⁺ release, malondialdehyde (MDA), reactive oxygen species (ROS), and GSH, were determined using the Iron Assay Kit (ab83366, Abcam, United Kingdom), lipid peroxidation (MDA) assay (ab118970, Abcam), DCF ROS/RNS Assay (ab238535, Abcam), and GSH/GSSG Ratio Detection Assay Kit II (ab205811, Abcam), respectively.

Total RNA was extracted from the cells using Trizol reagent (Invitrogen, Waltham, MA, United States). Approximately 500 ng of RNA was used for the reverse transcription reaction with PrimeScript RT Master Mix [Takara Biotechnology (Dalian) Co., Ltd., Japan]. Then, with GAPDH as the internal control, real-time quantitative polymerase chain reaction (RT-qPCR) was performed using Premix Ex Taq^TM II [Takara Biotechnology (Dalian) Co., Ltd.] with the Roche Light Cycler 480 Real-Time PCR system. The sequence of primers is as follows: AIFM2 forward: TTACAAGCCAGAGACTGACCAA, reverse: ACAAGGCCTGTCACTGAAGAG; NFS1 forward: CA CTCCCGGACACATGCTTAT, reverse: TGTCTGGGTGGTGAT CAAGTG; GAPDH forward: ATCATCAGCAATGCCTCCTG, reverse: ATGGACTGTGGTCATGAGTC.

Protein samples (50 μg) were separated by sodium dodecyl sulfate–polyacrylamide gel electrophoresis in 4–12% gel and transferred to nitrocellulose membranes for reaction with antibodies against target genes. Secondary antibodies, i.e., horseradish peroxidase-conjugated rabbit anti-goat and rabbit anti-mouse IgG, were detected by using SuperSignal Chemiluminescent substrate (Pierce Biotechnology, Rockford, IL, United States).

Tumors and adjacent normal tissues of 52 patients with PCa who underwent radical prostatectomy at the Beijing Hospital, Xiangya Hospital, and the Second Xiangya Hospital were examined by immunohistochemistry (IHC). Sections were incubated with anti-AIFM2 and anti-NFS1 (1:100 dilution). The sections were scored for staining intensity according to the following scale: negative, no staining; weak, weak staining (light yellow); moderate, moderate staining (yellowish brown); and strong, strong staining (brown). The experiments were reviewed and approved by the Ethics Committee of the Beijing Hospital (2018BJYYEC-085-03). The patients provided their written informed consent to participate in this study.

Prostate cancer cells in the logarithmic growth stage were put into six-well plates. The culture continued for 1–3 weeks, during which the solution was changed every 3 days, and the cell status was observed every day until the number of cells in most single clones was greater than 50. Then the cells were cleaned with phosphate-buffered saline buffer, made in methanol approximately 30 min, and stained with crystal violet dye at a dose of 1%. Thereafter, the number of colonies was calculated.

All animal experiments were approved by the Institutional Ethics Committee of Xiangya Hospital, Central South University. A subcutaneous xenograft model was established by subcutaneously injecting male nude mice with 1 × 10⁶ DU145 cells on the right side. The tumor volume was measured with calipers and repeatedly measured every 7 days (length × width²)/2. At 28 days following implantation, the mice were euthanized by cervical dislocation. Then, the xenografts were removed, fixed, weighed, photographed, and preserved.

All data sorting and analyses were completed by the R 3.6.1 software. For continuous variables with normal distribution and homogeneity of variance, an independent sample t-test was used; otherwise, Wilcoxon rank-sum test was selected. Pearson correlation coefficient test was used to analyze the correlation. A value of p < 0.05 was considered significant. Based on the results of the multivariate Cox proportional hazards analysis, a nomogram was developed to predict 1-, 3-, and 5-year bRFS rates. To evaluate the prediction performance of the nomogram, the consistency index (C index), ROC, and calibration curve were used to evaluate the nomogram.

Figure 1 shows the flow diagram of this study. Among 84 ferroptosis-related genes, 31 were underexpressed, and 27 were highly expressed in tumor tissues compared with the corresponding normal tissues in the TCGA cohort (FDR < 0.05) (Figure 2A and Supplementary Table 1). To identify prognostic genes associated with the bRFS, the univariate Cox regression analysis suggested that 22 ferroptosis-related genes correlated with the prognosis of PCa (Supplementary Table 2). Only four genes serve as protective factors, while the other 18 genes were risk factors (Figure 2B). Taking the intersection of the two lists of genes, 17 ferroptosis-related genes were identified as candidate genes for the subsequent construction of a prognostic model (Figure 2C). The expression of these genes in tumor tissues and adjacent normal tissues is shown in Figure 2D.

These 17 bRFS-related genes were uploaded to STRING to construct a protein–protein interaction (PPI) network (Figure 3A). The correlation analysis also showed that these genes were strongly correlated at the transcriptional level (Figure 3B). To screen the colinearity of these 17 genes, the LASSO Cox regression analysis was performed to determine the real bRFS-affecting factors and finally identified a prognostic panel of 9 ferroptosis-related genes (Figures 3C,D). The following formula was used to calculate the risk assessment score: risk score = AIFM2 × (0.19867) + AKR1C1 × (0.05148) + AKR1C2 × (0.04941) + CBS × (0.02658) + FANCD2 × (0.00202) + FTH1 × (0.00027) + G6PD × (0.00189) + NFS1 × (0.13708) + SLC1A5 × (-0.00602).

Based on the signature calculation formula, all patients in the TCGA cohort were divided into high- and low-risk groups according to their median risk score. The K-M curve showed that the high-risk group had a higher probability of BCR (p < 0.001) (Figure 4A). When evaluating bRFS prediction, the 1-, 3-, and 5-year AUCs of the developed gene signature were 0.680, 0.738, and 0.767, respectively (Figure 4B). PCA demonstrated that patients in various risk groups showed different two-dimensional spatial distributions (p < 0.001) (Figure 4C). When all patients were ranked according to their risk scores, the proportion of patients with BCR in the high-risk group was significantly higher than that in the low-risk group (Figure 4D). Univariate and multivariate Cox regression analyses were implemented to evaluate whether the model was an independent predictor among other clinical factors such as age, prostate-specific antigen (PSA), Gleason score (GS), and staging. We found that the risk score remained independently associated with bRFS not only at the univariate but also at the multivariate analysis when combined with all the clinical features (p < 0.05) (Figures 4E,F and Supplementary Table 3).

This study evaluated the MSKCC to verify the predictive robustness of this model. To avoid false positives caused by information bias, the same statistical methods were used. Surprisingly, the validation results were highly consistent with that in the TCGA cohort. Survival analysis also demonstrated that a high-risk score was associated with the poor bRFS of patients with PCa (Figure 5A). The AUC values for the prognostic model for bRFS were 0.766 at 1 year, 0.729 at 3 years, and 0.726 at 5 years (Figure 5B). The PCA results also suggested that different risk subgroups showed significant discrete tendency directly in the two-dimensional plane (Figure 5C). The BCR rate in the high-risk group was significantly higher than that in the low-risk group (Figure 5D). The risk score was still an independent prognostic factor in the MSKCC (Figures 5E,F and Supplementary Table 3).

To test the predictive power of the prognostic risk model for clinical features, we correlated the risk score with the clinical features (i.e., age, BCR state, PSA, GS, WHO ISUP classification, and T-staging). The risk grouping system showed a significant correlation with the BCR state, GS, WHO ISUP classification, and T-staging but not with age and PSA in both cohorts. The high-risk group showed a higher percentage of BCR and higher WHO ISUP classification, GS, and T-staging than the low-risk group. The heat map of gene expression showed that risk factors AIFM2, AKR1C1, AKR1C, CBS, FANCD2, FTH1, G6PD, and NFS1 were highly expressed in the high-risk group, while only SLC1A5, as a protective factor, was highly expressed in the low-risk group in both cohorts (Figure 6 and Supplementary Table 4).

To explore the mechanism of this risk model, GSEA was used to analyze the potential biological processes and pathways. GSEA results showed that “steroid hormone secretion,” “regulation of endocrine processes,” and “positive regulation of humoral immune response” were significantly enriched in the high-risk group (Figures 7A–C), while the low-risk group was significantly enriched in “androgen response” (Figure 7D). In addition, the levels of Ki67 expression were compared between the high- and low-risk groups, which revealed that the level of Ki67 expression was higher in the high-risk group (Figure 7E). The FPI model constructed by Liu et al. (2020c) showed that a high FPI predicts poor prognosis in various tumors, and it was associated with many important metastasis and immune-related pathways. Therefore, the relationship between the risk score of our prognostic model and FPI was explored. Similarly, the FPI of the high-risk group was significantly higher than that of the low-risk group (Figure 7F). Notably, all the presented results are consistent in the TCGA and MSKCC cohorts (Supplementary Table 5).

Given the above GSEA results and results of previous studies that iron death may be associated with immune cells, immune function, and immune microenvironment, we compared the differences in immune-related parameters between high- and low-risk groups in terms of the infiltration of 16 immune cells and activity of 13 immune-related pathways. In this study, most immune cell and immune function scores were higher in the high-risk group than in the low-risk group (Figure 8). In particular, aDCs, Tfh, TIL, APC co-inhibition, checkpoint, cytolytic activity, T-cell co-inhibition, and T-cell co-stimulation scores were significantly enriched in the high-risk group (Supplementary Table 6). We used the “ESTIMATE” package to evaluate TME and obtained similar results. With the increase in the risk score, the immune/stromal/ESTIM scores as a predictor of TME also increased, and the purity of the tumor decreased (Figure 9A and Supplementary Table 7).

Alongside radical surgery and radiotherapy, drug therapy, including endocrine therapy and chemotherapy, is an important treatment of PCa. In the GSEA, “androgen response” was significantly enriched in the low-risk group, so we explored whether patients in the high- and low-risk groups responded differently to bicalutamide. The estimated IC50 value demonstrated that the low-risk group in both cohorts had a better response to bicalutamide (p < 0.0001) (Figure 9B). Similarly, the level of Ki67 expression was higher in the high-risk group than in the low-risk group, suggesting that the high-risk group had more active tumor cell proliferation, faster tumor growth, and worse tissue differentiation than the low-risk group. This also raises the question on whether the response to chemotherapy is different between the high- and low-risk groups. Docetaxel, the most classic chemotherapy drug for PCa, was also found to have a better response in the low-risk group of the two cohorts, but the difference was not significant in the MSKCC cohort (p = 0.15) (Figure 9C).

To improve the accuracy of the risk score model in predicting bRFS, we combined the risk scores of the two cohorts and constructed a nomogram with conventional clinicopathological features (Figure 10A). The C-index (0.751) showed good predictive accuracy for the nomogram. The calibration curves had good linearity for the 1-, 3-, and 5-year bRFS (Figure 10B). The nomogram scores achieved a higher AUC–ROC than all other clinicopathological parameters in predicting 1-, 3-, and 5-year bRFS (AUC at 1 year, 0.774; AUC at 3 years, 0.768; AUC at 5 years, 0.789) (Figure 10C).

To further explore the mechanisms involved in the risk signature, we focused on the nine ferroptosis-related genes. The coefficients of AIFM2 and NFS1 genes were significantly higher than those of the other seven genes involved in the construction of the signature, and a larger weight often indicates that the gene is more prognostically important (Figure 11A). Therefore, AIFM2 and NFS1 became the focus of the subsequent research. First, we evaluated the differential expression of AIFM2 and NFS1 between the tumor tissues and corresponding normal tissues in multiple datasets. The results of a meta-analysis after combining 15 independent datasets showed that the expression of AIFM2 in the tumor was lower than that in normal tissues, while opposite results were obtained in NFS1 (p < 0.001) (Figure 11B). Furthermore, we found that patients with a high expression of AIFM2 or NFS1 had earlier onset of biochemical relapse and that a high expression of AIFM2 or NFS1 was associated with higher GS and WHO ISUP grade, poorer clinicopathological stage, and higher positive rate of the surgical margin, but some of the results were not significant because of the small sample size (Figure 11C). To verify these findings, IHC of these two genes was performed in 52 pairs of tumor tissues and corresponding normal tissues after radical prostatectomy. The actual results were consistent with the expected results. The expression of AIFM2 was lower in tumor tissues than in normal tissues, but the high expression of AIFM2 tended to have higher GS (p < 0.01). NFS1 is highly expressed in tumors; likewise, highly expressed NFS1 tends to have higher GS (p < 0.001) (Figure 12A).

We assessed the sensitivity of PCa cells (DU145) to ferroptosis. Erastin (ferroptosis inducer) treatment decreased the viability of DU145 cells in a dose-dependent manner. Compared with normal human prostate cells (PWR1E), the activities of DU145 cells were significantly decreased at any concentrations of erastin (p < 0.001) (Figure 12B). By using BODIPY^TM 581/591 C11 as a lipid peroxidation probe, we confirmed that the DU145 cell lines induced by erastin produced more lipid peroxidations (green) than the PWR1E cell lines (Figure 12C). In addition, the levels of Fe²⁺ release, MDA, DCF, and GSH in DU145 cells were measured for 24 h under different treatment conditions. After induction by erastin, Fe²⁺, MDA, and DCF levels in DU145 cells were significantly increased, while the GSH level was significantly decreased. This trend can be partially rescued by ferrostatin-1, a ferroptosis inhibitor (Figure 12D).

Moreover, we detected a difference in the expressions of the AIFM2 or NFS1 genes between PCa cell lines and normal cell lines, which was consistent with the analysis results of the tissue samples (Figure 13A). To evaluate the potential role of AIFM2 and NFS1 in the regulation of PCa ferroptosis, human PCa DU145 cells were treated with control shRNA or target genes shRNA (Figure 13A). Erastin-induced lipid peroxidation measured with BODIPY^TM 581/591 C11 could be promoted by sh-AIFM2 and sh-NFS1 (Figure 13B). As expected, AIFM2 and NFS1 knockdowns by shRNA notably increased the levels of Fe²⁺, MDA, and DCF, but decreased the level of GSH (but not significantly) (Figure 13C). The above results indicated that AIFM2 or NFS1 knockdown could promote the ferroptosis of PCa cell lines. To date, the phospholipid hydroperoxide-reducing enzyme glutathione peroxidase 4 (GPX4) has been known to be the main enzyme that protects against ferroptosis (Stockwell et al., 2017), and Acyl-CoA synthetase long-chain family member 4 (ACSL4) is also a key contributor and regulator of ferroptosis, which determines the sensitivity of ferroptosis (Yuan et al., 2016; Doll et al., 2017). Therefore, we tried to explore whether the AIFM2 and NFS1 knockdowns affect ferroptosis by affecting the expressions of GPX4 and ACSL4 proteins. However, we did not observe significantly positive results (Figure 13D), which suggests the involvement of other pathways. In addition, colony formation assays indicated that AIFM2 or NFS1 knockdown significantly inhibited cell colony formation. However, it is unclear whether the reduced clonogenic capacity is due to the increase in ferroptosis sensitivity or the decrease in tumor cell proliferation. Therefore, ferrostatin-1 was added after NFS1/AIFM2 knockdown to inhibit ferroptosis. The addition of ferrostatin-1 partially rescued the reduced clonogenic capacity in both sh-AIFM2 and sh-NFS1 cells (Figure 13E). This illustrated the decreased ability to form colonies after NFS1/AIFM2 knockdown, on the one hand, by promoting ferroptosis, and on the other, by inhibiting cell proliferation ability. This was indeed confirmed by subsequent detection of proliferation markers Ki67 and proliferating cell nuclear antigen (PCNA), which were significantly downregulated in both sh-AIFM2 and sh-NFS1 cells (Figure 13F). Finally, we further explored the effects of AIFM2 and NFS1 on the development of PCa in vivo. We found that AIFM2 and NFS1 knockdowns both suppressed PCa cell tumorigenicity (Figure 13G).