- 1Department of Urology, Urology Research Institute, The First Affiliated Hospital, Fujian Medical University, Fuzhou, China
- 2Department of Radiotherapy, The First Affiliated Hospital of Fujian Medical University, Fuzhou, China
- 3Fujian Key Laboratory of Precision Medicine for Cancer, The First Affiliated Hospital, Fujian Medical University, Fuzhou, China
Objective: To identify ferroptosis-related molecular clusters, and to develop and validate a ferroptosis-based molecular signature for predicting biochemical recurrence-free survival (BCRFS) and tumor immune microenvironment of prostate cancer (PCa).
Materials and Methods: The clinical data and transcriptome data of PCa were downloaded from TCGA and GEO database. Ferroptosis-related genes (FRGs) were obtained from FerrDb database. We performed consensus clustering analysis to identify ferroptosis-related molecular subtypes for PCa. Univariate and multivariate Cox regression analysis were used to establish a ferroptosis-based signature for predicting BCRFS. Internal verification, external verification and subgroup survival analysis were then successfully performed.
Results: There was a total of 40 differentially expressed FRGs in PCa. We then identified three ferroptosis-related molecular clusters of PCa, which have significantly different immune infiltrating cells, tumor immune microenvironment and PD-L1 expression level. More importantly, a novel ferroptosis-based signature for predicting BCRFS of PCa based on four FRGs (including ASNS, GPT2, NFE2L2, RRM2) was developed. Internal and external verifications were then successfully performed. Patients with high-risk score were associated with significant poor BCRFS compared with those with low-risk score in training cohort, testing cohort and validating cohort, respectively. The area under time-dependent Receiver Operating Characteristic (ROC) curve were 0.755, 0.705 and 0.726 in training cohort, testing cohort and validating cohort, respectively, indicating the great performance of this signature. Independent prognostic analysis indicated that this signature was an independent predictor for BCRFS of PCa. Subgroup analysis revealed that this signature was particularly suitable for younger or stage T III-IV or stage N0 or cluster 1-2 PCa patients. Patients with high-risk score have significantly different tumor immune microenvironment in comparison with those with low-risk score. The results of qRT-PCR successfully verified the mRNA expression levels of ASNS, GPT2, RRM2 and NFE2L2 in DU-145 and RWPE-1 cells while the results of IHC staining exactly verified the relative protein expression levels of ASNS, GPT2, RRM2 and NFE2L2 between PCa and BPH tissues.
Conclusions: This study successfully identified three ferroptosis-related molecular clusters. Besides, we developed and validated a novel ferroptosis-based molecular signature, which performed well in predicting BCRFS and tumor immune microenvironment of PCa.
Introduction
As the most prevalent non-cutaneous carcinoma of males in the world, prostate cancer (PCa) is considered as the fifth leading cause of death (Xu et al., 2020a; Huang et al., 2020). Radical prostatectomy (RP) remains the standard curable strategy for localized PCa; however, there is still a certain proportion of patients developing biochemical recurrence (BCR) (Venclovas et al., 2019), which indicated a possibility of underlying clinical metastases and poor prognosis (Xu et al., 2016; Xu et al., 2020b). The early recognition of BCR is of great significance for the subsequent treatment strategies of PCa patients. Although several studies have investigated diverse biomarkers or clinical factors to predict BCR, there are no recognized molecular subtypes and prognostic signature related to BCR for PCa (Arora and Barbieri, 2018).
Recently, immunotherapy has been regarded as an indispensable treatment modality for malignancies (Riley et al., 2019). With the Food and Drug Administration (FDA) approvals of sipuleucel-T into PCa treatment, immunotherapy has improved outcomes of carefully selected patients (Bilusic et al., 2017; Milonas et al., 2019). Furthermore, several studies have revealed that immune responses and clinical responses could be tremendously promoted by combining various immunotherapeutic agents with hormonal therapy, chemotherapy or DNA-damaging agents (Bilusic et al., 2017; Gamat and McNeel, 2017). Immunotherapy showed a promising future in patients with PCa (Bilusic et al., 2017). Therefore, in addition to discovering the optimal treatment combination, there is an urgent need to develop potential biomarkers for immunotherapy response and tumor immune microenvironment.
Ferroptosis is considered a programmed cell death pattern, which is characterized by iron-dependent peroxidation (Zuo et al., 2020). Increased cellular activity of carcinoma cells leads to high susceptibility to ferroptosis (Kwon et al., 2020). In recent years, studies related to the role of ferroptosis in carcinogenesis and progression have gained momentum (Stockwell et al., 2017). PCa has been recognized as one of the ferroptosis-related diseases (Ghoochani et al., 2021). Previous studies have demonstrated that the abnormal lipid homeostasis is associated with tumorigenesis of PCa and progression to castration-resistant prostate cancer (Blomme et al., 2020). However, to our knowledge, studies focusing on the correlation of ferroptosis with biochemical recurrence and anti-tumor immunology of PCa were rather limited.
In this study, we identified three ferroptosis-related molecular clusters for PCa. Moreover, interestingly, we successfully developed a ferroptosis-based signature for predicting biochemical recurrence-free survival (BCRFS) of PCa. Internal and external verifications were then resoundingly performed. Besides, the correlations of the ferroptosis-related molecular clusters and signature with tumor immune microenvironment, cancer stemness and drug sensitivity were particularly explored. We then conducted qRT-PCR in DU-145 and RWPE-1 cells and immunohistochemical (IHC) staining in PCa and BPH tissues to validate the mRNA and protein expression levels of four risk FRGs.
Materials and Methods
Data Collection
This study was approved by the Institutional Ethics Committee of First Affiliated Hospital of Fujian Medical University and was conducted in accordance with 1964 Helsinki declaration and its later amendments. Written informed consent was provided by all included patients before sample collection. Fresh postoperative PCa tissues and benign prostatic hyperplasia (BPH) tissues were obtained from three PCa patients and three BPH patients who have been pathologically diagnosed by biopsy and were treated in the Department of Urology, the First Affiliated Hospital of Fujian Medical University.
Transcriptome profiles of 499 PCa cases and 52 normal cases were downloaded from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). In this study, considering that there were only 70 samples that have the data type of “days_to_first_biochemical_recurrence”, we combined “biochemical_recurrence” and “A8_New_Event_Type” as the BCR status, and combined “days_to_first_biochemical_recurrence” and “A8_New_Event_Time” as the time to BCR. (1) In the 70 samples that have data type of “days_to_first_biochemical_recurrence”, we used “days_to_first_biochemical_recurrence” as the time to BCR and used “biochemical_recurrence” as the BCR status. (2) In those that did not have data type of “days_to_first_biochemical_recurrence”, we utilized “A8_New_Event_Time” as the time to BCR. (3) If we utilized “A8_New_Event_Time” as the time to BCR, we preferentially used “A8_New_Event_Type” as the BCR status; in those without data of “A8_New_Event_Type”, we used “biochemical_recurrence” as the BCR status. BCR-free survival (BCRFS) was defined as the interval from radical treatment to first BCR or death.
A total of 405 PCa cases in TCGA database had complete data of transcriptome, BCR status and BCR-free time. Dataset of GSE70770 (including 203 PCa cases with complete data of transcriptome, BCR status and BCR-free time) was obtained from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/).
FerrDb database (http://www.zhounan.org/ferrdb/) is the world’s first database related to ferroptosis regulators and markers and ferroptosis-disease associations. We downloaded a total of 382 ferroptosis-related genes (FRGs) from FerrDb database, including 150 ferroptosis drivers (genes promoting ferroptosis.), 109 ferroptosis suppressors (genes preventing ferroptosis) and 123 ferroptosis markers (genes indicating the occurrence of ferroptosis). (Supplementary Table S1).
Identification of Differentially Expressed FRGs and Functional Enrichment
The mRNA expression matrix of 382 FRGs in TCGA cohort was extracted using “limma” R package. We then utilized Wilcoxon test to filter differentially expressed FRGs (DEFRGs) between normal samples and PCa samples in TCGA cohort. The cut-off value was set as FDR (false discovery rate) < 0.05 and log2 |fold change (FC)| > 0.5. (Mock et al., 2016). All DEFRGs in PCa were shown in heatmap using “pheatmap” R package. R packages “clusterProfiler” and “org.Hs.eg.db” were utilized for perform Gene Oncology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment.
Consensus Clustering Analysis Identifying Ferroptosis-Related Molecular Clusters
We next conducted univariable Cox regression analysis to screen prognostic DEFRGs associated with BCRFS. The cut-off p value was set as 0.05. Then consensus clustering analysis was performed for identifying ferroptosis-related molecular clusters using R package “ConsensusClusterPlus”.
R packages “survival”, “survminer” and “pheatmap” were used to explore the relationship of ferroptosis-related molecular subtypes with clinicopathologic characteristics (including BCRFS, T stage, N stage and age). The relationship of ferroptosis-related molecular clusters with PD-L1 gene expression level was explored. Furthermore, the ESTIMATE algorithm was utilized to evaluate tumor microenvironment (TME) scores (Ke et al., 2021), and the CIBERSORT method was used to calculate the score of 22 types of immune infiltrating cells (Chen et al., 2019). The relationships of ferroptosis-related molecular clusters with PCa immune microenvironment and immune infiltration cells were especially investigated. p value < 0.05 was considered statistically significant.
Development and Verification of a Novel Ferroptosis-Based Signature for BCRFS
In this study, the TCGA cohort (405 PCa cases) was randomly split into training cohort (204 PCa cases) and testing cohort (201 PCa cases). The clinicopathological data should be comparable between the two groups. The 203 PCa cases in GEO database (GEO cohort) were used as validating cohort. We performed univariate and multivariate Cox regression analysis in training cohort to establish a novel ferroptosis-based signature for predicting BCRFS. The cutoff p value of univariate analysis was 0.05. On the basis of the median risk score derived from the novel ferroptosis-based signature, we divided PCa patients into high-risk and low-risk subgroup. We performed survival analysis, time-dependent the receiver operating characteristic (ROC) curve analysis and independent prognostic analysis to validate the performance of this novel ferroptosis-based signature.
More importantly, internal and external verifications were conducted in testing cohort and validating cohort, respectively. Survival analysis and time-dependent ROC curve analysis were also performed in testing cohort and validating cohort. The risk distribution of TCGA cohort, training cohort, testing cohort and validating cohort was presented using “pheatmap” R package. Especially, subgroup survival analysis was conducted in older and younger patients, stage T I-II and stage T III-IV patients, stage N0 and stage N1 patients, and different molecular subtypes. p value < 0.05 was considered statistically significant.
Association of the Ferroptosis-Based Signature with Clinicopathologic Features, Tumor Stemness Scores, Immune Infiltrating Cells and Functional Enrichment.
We investigated the association of the ferroptosis-based signature with clinicopathologic features, immune cells infiltration, immune-related pathways activity and tumor stemness scores. The scores of 16 immune infiltrating cells and 13 immune function activity in TCGA cohort were calculated by single-sample gene set enrichment analysis (ssGSEA) using R package “gsva”. Moreover, gene set enrichment analysis (GSEA) was performed to investigate the underlying mechanisms related to this ferroptosis-based signature, as previous described (Cheng et al., 2021).
Drug Sensitivity of Risk DEFRGs
CellMiner database (https://discover.nci.nih.gov/cellminer/home.do) is a database and query tool designed for the cancer research community to facilitate integration and study of molecular and pharmacological data. We applied CellMiner database to demonstrate whether these risk DEFRGs could predict the sensitivity of anticancer drugs. To enhance the relationship of ferroptosis risk genes and clinical application, only FDA approved drugs and drugs under clinical trials were included in the analysis. Spearman correlation analysis was performed to determine the correlations between the expression levels of ferroptosis risk genes and drug sensitivity.
Validation of Risk DEFRGs Using UALCAN Database
UALCAN database (http://ualcan.path.uab.edu/) is a portal for facilitating tumor gene expression. The mRNA expression levels of risk DEFRGs between normal and tumor tissues were validated using UALCAN database.
Cell Culture
Human PCa cell line DU145 (Procell CL-0075) and normal prostatic epithelial cell line RWPE-1(Procell CL-0200) were kindly provided by Procell Life ScienceandTechnology Co., Ltd. The DU145 cell line was cultured in MEM medium (Biological Industries, Beit HaEmek, Israel), supplemented with 10% fetal bovine serum in a standard humidified incubator at 37°C in a 5% CO2 atmosphere. The RWPE-1 cell line was cultured in RPMI-1640 medium (Biological Industries, Beit HaEmek, Israel) supplemented with 10% fetal bovine serum at 37°C in a humidified atmosphere containing 5% CO2.
Total RNA Extraction and Quantitative Reverse-Transcription Polymerase Chain Reaction (qRT-PCR).
Most importantly, to validate the different mRNA expression levels of risk DEFRGs, we conducted qRT-PCR in DU-145 (prostate cancer cells) and RWPE-1 (normal prostatic epithelial cells) cell lines. According to the manufacturer’s protocol, the total RNA of risk DEFRGs was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA). The RNA purity was detected using NanoDrop 2000 spectrometer (Thermo Fisher Scientific, Waltham, MA). Reverse transcription reactions using TransScript® Green One-Step qRT-PCR SuperMix (TransGen Biotech, Beijing, China) were conducted in only one step according to its specification. The qRT-PCR was carried out to detect the expression levels of the four genes with the Step One PlusTM PCR System (Applied Biosystems) using Taq Pro universal SYBR qPCR Master Mix (Vazyme, Nanjing, China) according to the manufacturer’s protocol. Results were normalized to the expression of GAPDH. The results of qRT-PCR were calculated by the 2−∆∆Ct method: ∆Ct = Ct (risk DEFRGs)—Ct (GAPDH) and ∆∆Ct = ∆Ct (DU-145 cells)—∆Ct (RWPE-1 cells).
The primer sequences were as follows:
GAPDH: forward 5-GGTGTGAACCATGAGAAGTATGA-3, reverse 5-GAGTCCTTCCACGATACCAAAG-3, product size 123 bp;
ASNS: forward 5-CTCCGCGCAGATCGAACTAC-3, reverse 5-TCTTTTGGTCGCCAGAGAATC-3, product size 197 bp;
GPT2: forward 5- CAAGAAGGTGCTGTACGAGATG-3, reverse 5- CTCCATGTAGCCTCCTCTGTAA-3, product size 118 bp;
NFE2L2: forward 5- ATGCCACAGTCAACACAGATT, reverse 5-GCCCATTTAGAAGTTCAGAGAGT-3, product size 126 bp;
RRM2: forward 5-TTGCCTGTGAAGCTCATTGG-3, reverse 5-CCTCTGATACTCGCCTACTCTC-3, product size 192 bp.
Immunohistochemical (IHC) Staining
Three PCa tissues and BPH tissues were made into paraffin-embedded wax blocks, and then processed into 5 μM thick sections, as reported previously (Broggi et al., 2021). Briefly, after roasted at 60°C for 30 min, the tissue sections were dewaxed and hydrated. Next, 3% hydrogen peroxide was used to block the endogenous peroxidase and the citric acid buffer was used to repair antigen in boiling water for 5 min. The sections were blocked with normal sheep serum for blocking (ZLI-9022, ZSGB-BIO, Beijing, China) for 15 min and then separately incubated with the anti-ASNS (14681-1-AP, Proteintech), anti-NFE2L2 (YT3189, Immunoway), anti-RRM2 (11661-1-AP, Proteintech) and anti-GPT2 antibody (16757-1-AP, Proteintech) at 4°C overnight. After that, the secondary antibody was added for 1 h. Through DAB staining (ZLI-9017, ZSGB-BIO, Beijing, China) for 3 min, hematoxylin immersion for several minutes, differentiation liquid differentiation, and finally neutral gum sealing piece. We then used the Image J v1.51k software (https://imagej.nih.gov/ij/, Wayne Rasband, United States) to convert gray values to optical density (OD) values. After adjusting the appropriate threshold, we could get the relevant data about integrated optical density (IOD) and positive area. Finally, the average optical density (AOD) was used to assess the relative expression (AOD = IOD/area).
Statistical Methods
Statistical analysis was performed utilizing R programming language. Univariate and multivariate Cox regression analyses were performed to establish a ferroptosis-related prognostic signature for predicting BCRFS of PCa. Univariate and multivariate independent prognostic analyses were used to demonstrate whether this ferroptosis-related prognostic signature was an independent predictor of BCRFS. Survival analysis and time-dependent the ROC curve were performed to explore the performance of this ferroptosis-related prognostic signature. p value < 0.05 was considered statistically significant.
Results
Differentially Expressed FRGs and Functional Enrichment
The study flow chart was presented in Figure 1. The clinicopathologic data of training and testing cohort were showed in Table 1. Finally, a total of 40 DEFRGs were finally identified, including 26 downregulated genes and 14 upregulated genes. (Table 2). The heatmap and correlation network of these 40 FRGs were showed in Figure 2A and Figure 2B, respectively. The results of functional enrichment analysis of 40 DEFRGs were presented in Figures 2C,D. According to results of the GO functional analysis, response to oxidative stress and organic acid biosynthetic process were the main GO terms of biology process (BP); organelle outer membrane and lipid droplet were the main GO terms of cellular component (CC); peroxidase activity and antioxidant activity were the main GO terms of molecular function (MF). Moreover, KEGG analysis revealed that these 40 DEFRGs were mainly enriched in response to oxidative stress, organic acid biosynthetic process, fatty acid metabolic process, carboxylic acid biosynthetic process and response to toxic substance.
FIGURE 2. The heatmap (A), correlation network (B) and GO (C) and KEGG enrichment (D) of 40 ferroptosis-related differentially expressed genes. The expression correlation of PD-L1 with 14 ferroptosis-related prognostic differentially expressed genes (E). The expression level of PD-L1 between PCa tissue and normal tissue (F).
Univariable Cox regression analysis using the whole TCGA cohort identified a total of 14 DEFRGs related to BCRFS, including PTGS2, DUSP1, ALB, ASNS, ATF3, GDF15, FTL, NFE2L2, RRM2, TP63, ACSL3, CAV1, CDO1 and ANO6. (Table 3). The expression association of PD-L1 and these 14 FRGs was demonstrated in Figure 2E. The results showed that the expression level of PD-L1 was significantly associated with the expression level of PTGS2, DUSP1, ATF3, FTL, NFE2L2, TP63, ACSL3, CAV1, and ANO6. The expression of PD-L1 in PCa tissue was significantly decreased compared with normal tissue (Figure 2F).
TABLE 3. Univariable Cox regression analysis to identify DEFRGs related to biochemical recurrence-free survival in whole TCGA cohort.
Identification of Three Ferroptosis-Related Molecular Clusters
Then, we performed consensus clustering analysis to develop a total of three ferroptosis-related molecular clusters of PCa, including 146 cases of cluster 1, 121 cases of cluster 2, 138 cases of cluster 3 (Figures 3A–C). The difference of BCRFS among these ferroptosis-related molecular clusters was not statistically significant (Figure 3D). These three molecular clusters have significantly different age (Figure 3E).
FIGURE 3. Identification of three ferroptosis-related molecular clusters using consensus clustering analysis (A,B,C). Comparison of biochemical recurrence free survival among these three clusters (D). The correlation heatmap between the ferroptosis-related molecular clusters and clinicopathologic features (E). The expression level of PD-L1 among these three clusters (F).
The expression level of PD-L1 in cluster 2 was significantly higher in comparison with that in cluster 3 and cluster 1. The expression level of PD-L1 in cluster 3 was significantly higher in comparison with that in cluster 1 (Figure 3F). The B cells naïve, plasma cells, T cell follicular helper, macrophage M2 and mast cell activated were the discrepant immune infiltrating cells among these three ferroptosis-related clusters (Figures 4A,B). The immune score in cluster 2, cluster 3 and cluster 1 was significantly decreased in sequence. The ESTIMATE score and stromal score in cluster 2 and cluster 3 were significantly increased compared with that in cluster 1; however, the ESTIMATE score and stromal score were not different between cluster 2 and cluster 3. The tumor purity in cluster 2 and cluster 3 was significantly decreased compared with that in cluster 1; however, the tumor purity was not different between cluster 2 and cluster 3. (Figures 4C–F).
FIGURE 4. Distribution of 22 immune cell among these three molecular clusters (A). Relationship between this ferroptosis-related molecular clusters and immune cell infiltration (B). Relationship between the ferroptosis-related molecular clusters and ESTIMATE score (C), immune score (D), stromal score (E), tumor purity (F) of tumor environment.
Development and Verification of a Novel Ferroptosis-Based Signature for Predicting BCRFS
First of all, we performed univariable Cox regression analysis in the training cohort and screen a total of eight DEFRGs according to the threshold p value of 0.05 (including ASNS, GPT2, NFE2L2, RRM2, TP63, CAV1, SAT1, ZEB1). Next, based on the results of univariable analysis, we performed stepwise multivariate Cox regression analysis in training cohort to develop a novel ferroptosis-related prognostic signature and finally only four DEFRGs were included in the final model (including ASNS, GPT2, NFE2L2, RRM2). The calculating formula of risk score was as follows: risk score = (0.161502387557219) ∗ ASNS + (−0.0463523973265746) ∗ GPT2 + (−0.0709927054603633) ∗ NFE2L2 + (0.144137051044523) ∗ RRM2. (Table 4). Then, we calculated the risk score of each case using above calculation formula in training cohort, testing cohort, whole TCGA cohort and validating cohort. All patients were divided into high-risk score group and low-risk score group according to the median risk score. The distribution of survival time, risk score and the expression heatmap of training cohort, testing cohort, whole TCGA cohort and validating cohort were demonstrated in Figure 5.
TABLE 4. Univariate and multivariate Cox regression analysis to developing ferroptosis based prognostic signature for prostate cancer using training cohort.
FIGURE 5. Development and validation of a novel ferroptosis-related prognostic signature for PCa. The expression heatmap, the distribution of risk score and survival time of training cohort (A), testing cohort (B), the whole TCGA cohort (C) and validating cohort (D).
The difference of BCRFS between low-risk and high-risk subgroup was statistically significant in training cohort (p < 0.001), testing cohort (p = 0.022), whole TCGA cohort (p < 0.001) and validating cohort (p < 0.001), respectively. Patients with high-risk score were associated with significant poor BCRFS and a higher possibility of BCR risk in comparison with those with low-risk score in training cohort, testing cohort, whole TCGA cohort and validating cohort, respectively. The area under time-dependent ROC curve were 0.755, 0.705, 0.724, and 0.726 in training cohort, testing cohort, whole TCGA cohort and validating cohort, respectively, indicating the great performance of this novel ferroptosis-based signature for predicting BCRFS of PCa. (Figure 6). Univariate and multivariate independent prognostic analysis indicated that this novel ferroptosis-based signature was an independent predictor for BCRFS of PCa (Table 5, p < 0.05).
FIGURE 6. The survival analysis between high and low risk group, and corresponding area under ROC curve in training cohort (A,B), testing cohort (C,D), whole TCGA cohort (E,F), and validating cohort (G,H).
Subgroup Survival Analysis
Subgroup survival analysis revealed that this ferroptosis-based signature was particularly suitable for younger or T stage III-IV or stage N0 PCa patients. High-risk score was associated with significantly poor BCRFS compared with low-risk score in PCa cases with age≤ 65 or T stage III-IV or N0 stage. However, the difference of BCRFS between low-risk and high-risk group was not statistically significant in patients with age > 65 or T stage I-II or N1 stage (Figures 7A–F).
FIGURE 7. The subgroup survival analysis between high and low risk group in patients of age ≤65 (A), age >65 (B), T stage I-II (C), T stage III-IV (D), N0 stage (E), N1 stage (F). The survival analysis between high and low risk group, and corresponding area under ROC curve in cluster 1 cohort (G,H), cluster 2 cohort (I,J), cluster 3 cohort (K,L).
We calculated the risk score for each patient in three ferroptosis-related molecular clusters, respectively. The results revealed that this ferroptosis-based signature was especially applicable in cluster 1 and cluster2 PCa patients. High-risk score was associated with significantly poor BCRFS compared with low-risk score in cluster 1 and cluster 2 PCa patients. The AUC for BCRFS prediction was 0.756, 0.694 in cluster 1 and cluster 2 patients, respectively. However, the difference of BCRFS between low-risk and high-risk group was not statistically significant in cluster 3. (Figures 7G–L).
Association of the Ferroptosis-Based Signature with Immune Cells Infiltration
We quantified the immune cells infiltration and immune functions using ssGSEA to further explore the association of this ferroptosis-based signature with immune status. Interestingly, the results revealed that patients with high-risk score had a lower infiltrating proportion of CD8+ T cells, mast cells, NK cells, and Treg in comparison with low-risk score. Moreover, the results also demonstrated that patients with high-risk score had a lower score of type II IFN response than those with low risk (Figures 8A,B).
FIGURE 8. Association of the ferroptosis-based signature with immune infiltrating cells (A), immune function activity (B).
Clinicopathologic Feature, Functional Enrichment and Cancer Stemness
The correlation heatmap among this prognostic ferroptosis-related signature, clinicopathologic features, ferroptosis-related molecular clusters and immune score was presented in Figure 9A. The results demonstrated that T stage, N stage, and subtypes were significantly different between high-risk patients and low-risk patients (p < 0.05). Besides, the correlation analysis revealed that risk score was positively associated with RNAss cancer stemness score (Figures 9B,C). KEGG functional enrichment for high-risk group and low-risk group was performed using GSEA method. RNA polymerase, base excision repair, DNA replication, pyrimidine metabolism and mismatch repair were top five KEGG enriched pathway in high-risk patients. Calcium signaling pathway, beta alanine metabolism, adipocytokine signaling pathway, nicotinate and nicotinamide metabolism and aldosterone regulated sodium reabsorption were top five KEGG enriched pathway in low-risk patients (Figure 9D).
FIGURE 9. Association of the ferroptosis-based signature with clinicopathologic features (A), cancer stemness (B,C), functional enrichment (D), anti-cancer sensitivity (E).
Drug Sensitivity of Risk DEFRGs
As revealed by Spearman correlation analysis, ASNS could predict the sensitivity of Megestrol acetate (Cor = 0.408, p = 0.001), Copanlisib (Cor = 0.361, p = 0.005), Carfilzomib (Cor = 0.308, p = 0.017), Irofulven (Cor = 0.302, p = 0.019), Oxaliplatin (Cor = 0.295, p = 0.022); RRM2 could predict the sensitivity of Nelarabine (Cor = 0.390, p = 0.002); NFE2L2 could predict the sensitivity of Cobimetinib (isomer 1) (Cor = 0.355, p = 0.005), Alectinib (Cor = 0.341, p = 0.008), Trametinib (Cor = 0.325, p = 0.011), BMN-673 (Cor = 0.324, p = 0.012), ARRY-162 (Cor = 0.309, p = 0.016), Encorafenib (Cor = 0.300, p = 0.020), Docetaxel (Cor = 0.299, p = 0.020), Nelfinavir (Cor = 0.296, p = 0.021); GPT2 could predict the sensitivity of 6-Thioguanine (Cor = 0.339, p = 0.008), Depsipeptide (Cor = 0.301, p = 0.019). (Figure 9F).
Validation of mRNA Expression Levels of Risk DEFRGs Using UALCAN Database
As indicated by UALCAN database, the mRNA expression levels of ASNS, GPT2, RRM2 in PCa tissue were significantly increased compared with that in normal tissue; however, the mRNA expression levels of NFE2L2 in PCa tissue was significantly decreased compared with that in normal tissue (Figures 10A–D).
FIGURE 10. The mRNA expression levels between PCa tissue and BPH tissue of ASNS (A), NFE2L2 (B), RRM2 (C) and GPT2 (D). The mRNA expression levels between PCa cells and normal prostatic epithelial cells of ASNS (E), NFE2L2 (F), RRM2 (G) and GPT2 (H). The protein expression levels between BPH tissue and PCa tissue of ASNS (I), NFE2L2 (J), GPT2 (K), RRM2 (L). *p < 0.05; **p < 0.01; ***p < 0.001.
Validation of mRNA Expression Levels of Risk DEFRGs in DU-145 PCa Cells Using qRT-PCR.
The results of qRT-PCR demonstrated that the mRNA expression levels of ASNS, GPT2, RRM2 in DU-145 PCa cell lines were significantly increased compared with that in RWPE-1 normal prostatic epithelial cells; however, the mRNA expression levels of NFE2L2 in DU-145 PCa cell lines was significantly decreased compared with that in RWPE-1 normal prostatic epithelial cells (Figures 10E–H). These results were consistent with UALCAN database.
Validation of Protein Expression Levels of Risk DEFRGs Using IHC
Moreover, the results of IHC suggested that the protein expression levels of ASNS, GPT2, RRM2 in PCa tissue were significantly increased compared with that in BPH tissue while the protein expression levels of NFE2L2 in PCa tissue was significantly decreased compared with that in BPH tissue (Figures 10I–L). These results were consistent with the results of qRT-PCR and UALCAN database.
Discussion
Previous studies have demonstrated an important role of ferroptosis in PCa (Ghoochani et al., 2021). The induction of ferroptosis could be a novel therapeutic strategy for advanced PCa (Ghoochani et al., 2021). Therapy-induced lipid uptake and remodeling significantly promoted the hypersensitivity of ferroptosis in PCa (Tousignant et al., 2020). Nassar et al. (2020). showed that DECR1 played a vital role in protecting PCa cells from ferroptosis by regulating polyunsaturated fatty acids oxidation. BCR is considered as a key event after radical treatment of PCa (Xu et al., 2016). Despite curative treatment, there were 20–40% of PCa patients experiencing BCR with 10 years; and BCR signified recurrent PCa. (Artibani et al., 2018). Although BCR alone might have no impact on quality of life or overall survival (Wei et al., 2019), BCR indicated the high probability of the onset of metastatic disease (Malik et al., 2019). However, there were few studies exploring whether there was a correlation between ferroptosis and BCR in PCa.
In this study, we first identified a total of three molecular clusters from the perspective of ferroptosis, which were different from previous studies. These three molecular clusters exhibited diverse clinical characteristics, PD-L1 expression levels and tumor immune microenvironment features. Moreover, interestingly, we performed multivariate Cox regression analysis to develop a novel ferroptosis-based signature for predicting BCRFS of PCa based on four DEFRGs (including ASNS, GPT2, NFE2L2, RRM2). Internal and external verifications were then resoundingly performed. Patients with high-risk score were associated with significant poor BCRFS in comparison with those with low-risk score in training cohort, testing cohort, whole TCGA cohort and validating cohort, respectively. The area under time-dependent ROC curve was 0.755, 0.705, 0.724 and 0.726 in training cohort, testing cohort, whole TCGA cohort and validating cohort, respectively, indicating the great performance of this novel ferroptosis-based signature for predicting BCRFS of PCa. Univariate and multivariate independent prognostic analyses indicated that this novel ferroptosis-based signature was an independent predictor for BCRFS of PCa.
This ferroptosis-based prognostic signature was composed of four ferroptosis-related genes: asparagine synthetase (ASNS), glutamic pyruvic transaminase 2 (GPT2), nuclear factor, erythroid 2 like 2 (NFE2L2) and ribonucleotide reductase regulatory subunit M2 (RRM2). Kanishka Sircar et al. (Sircar et al., 2012). revealed that ASNS was up-regulated castration-resistant stage of prostate cancer (CRPC), and that ASNS inhibitors might be a novel method via targeting CRPC cells. Harri M Itkonen et al. (Itkonen et al., 2016). demonstrated that coinstantaneous inhibition of GPT2 and OGT would suppress the growth and viability and additionally promote death response of PCa cells. Qiang Ju et al. (Ju et al., 2020). reported that NFE2L2 might be a potential prognostic indicator and associated with immune infiltration in brain lower grade glioma. Zhang et al. (2019) indicated that RRM2 might be involved in the ferroptosis of liver cancer. Mazzu et al. (2020) showed that RRM2 could be serve as a driver of aggressive PCa and contribute to immune escape. However, whether there is a relationship between ASNS, GPT2, NFE2L2, RRM2 and PCa ferroptosis was unknown. In this study, it was found that the mRNA expression levels of ASNS, GPT2, RRM2 in PCa tissue or DU145 PCa cells were significantly increased compared with normal tissue or RWPE-1 normal prostatic epithelial cells. Besides, the mRNA expression levels of NFE2L2 in PCa tissue or DU145 PCa cells was significantly decreased compared with that in normal tissue or RWPE-1 normal prostatic epithelial cells. Moreover, the protein levels of IHC staining showed that the protein expression levels of ASNS, GPT2 and RRM2 in PCa tissue were higher than that in BPH tissue, and that the protein expression levels of NFE2L2 in PCa tissue was lower than that in BPH tissue. As revealed by Spearman correlation analysis, ASNS could predict the sensitivity of Megestrol acetate and Copanlisib; RRM2 could predict the sensitivity of Nelarabine; NFE2L2 could predict the sensitivity of Cobimetinib (isomer 1) and Alectinib; GPT2 could predict the sensitivity of 6-Thioguanine and Depsipeptide. However, the role of ferroptosis-related genes, ASNS, GPT2, NFE2L2, RRM2, in PCa ferroptosis has not been elucidated. Additional research is required.
It has been reported that tumor-induced immunosuppressive status is crucial in PCa pathogenesis (Wong and Yu, 2021). Abnormal activating inhibitory immune checkpoint pathways have been regarded as the vital method of immune escape (Modena et al., 2016). However, it is unlikely that there were long-lasting responses from immunotherapy (including various immune checkpoint inhibitor) in the whole PCa patients (Claps et al., 2020). How to identify those suitable to receive immunotherapy is hot topic of discussion for urologist (Wong and Yu, 2021). Previous studies have recognized patient subsets who might benefit from immune checkpoint monotherapies, including those with high PD-L1 expression, hypermutation or high microsatellite instability, abundant germline/somatic DNA-repair gene mutations, intraductal/ductal histology, AR-V7-positive tumors however, it was still insufficient and these indicators could not adequately reflect the state of tumor immune microenvironment (Isaacsson Velho and Antonarakis, 2018; Vitkin et al., 2019). Interestingly, in the current study, we identified a total of three ferroptosis-related molecular clusters of PCa. Further investigation revealed that the expression levels of PD-L1 in cluster 2, cluster 3 and cluster 1 were significantly decreased in sequence. The results indicated that the sensitivity to immune check point inhibitors of cluster 2, cluster 3 and cluster 1 might decrease in sequence. Furthermore, we revealed that the B cells naïve, plasma cells, T cell follicular helper, macrophage M2 and mast cell activated were the discrepant immune infiltrating cells among these three ferroptosis-related clusters. The immune scores in cluster 2, cluster 3 and cluster 1 were also significantly decreased in sequence. The ESTIMATE score and stromal score in cluster 2 and cluster 3 were significantly increased compared with that in cluster 1. The tumor purity in cluster 2 and cluster 3 was significantly decreased compared with that in cluster 1. Patients with high-risk score have a lower infiltrating proportion of CD8+ T cells, mast cells, NK cells, and Treg and a lower score of type II IFN response in comparison with those with low-risk score.
Several limitations should be noted in this study. First of all, retrospective design with limited sample size would result in selection bias. Prospective real-world data are required to verify this ferroptosis-based molecular signature due to the retrospective methods of this studies. Besides, although qRT-PCR and IHC has been performed to validate the different mRNA and protein expression levels of risk DEFRGs, further wet experiment in vitro and in vivo into the underlying mechanism related to the correlation of ferroptosis with PCa immunotherapy and BCR are warranted.
Conclusion
This study identified a total of 40 differentially expressed FRGs and three ferroptosis-related molecular clusters in PCa. More importantly, a novel ferroptosis-based prognostic molecular signature of PCa based on four FRGs (including ASNS, GPT2, NFE2L2, RRM2) was developed, which has the great performance for predicting BCRFS. Our findings would contribute tremendously to the potential mechanism in PCa occurrence and development.
Data Availability Statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.
Ethics Statement
The studies involving human participants were reviewed and approved by the Institutional Ethics Committee of First Affiliated Hospital of Fujian Medical University. The patients/participants provided their written informed consent to participate in this study.
Author Contributions
Writing–original draft, Z-BK, QY, J-BS, and J-MZ. Writing–review and editing, X-YX and NX. Methodology, Z-BK, QY, and Q-SZ. Formal analysis, Z-BK, YW, and LS. Data curation, D-NC. Conceptualization, Z-BK. Visualization, X-DL. Project administration, X-YX and NX.
Funding
This study was supported by Natural Science Foundation of Fujian Province (Grant number: 2019J01445) and Fujian Provincial Health Technology Project (Grant number: 2019-ZQN-54).
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.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.774625/full#supplementary-material
References
Arora, K., and Barbieri, C. E. (2018). Molecular Subtypes of Prostate Cancer. Curr. Oncol. Rep. 20 (8), 58. doi:10.1007/s11912-018-0707-9
Artibani, W., Porcaro, A. B., De Marco, V., Cerruto, M. A., and Siracusano, S. (2018). Management of Biochemical Recurrence after Primary Curative Treatment for Prostate Cancer: A Review. Urol. Int. 100 (3), 251–262. doi:10.1159/000481438
Bilusic, M., Madan, R. A., and Gulley, J. L. (2017). Immunotherapy of Prostate Cancer: Facts and Hopes. Clin. Cancer Res. 23 (22), 6764–6770. doi:10.1158/1078-0432.Ccr-17-0019
Blomme, A., Ford, C. A., Mui, E., Patel, R., Ntala, C., Jamieson, L. E., et al. (2020). 2,4-dienoyl-CoA Reductase Regulates Lipid Homeostasis in Treatment-Resistant Prostate Cancer. Nat. Commun. 11 (1), 2508. doi:10.1038/s41467-020-16126-7
Broggi, G., Lo Giudice, A., Di Mauro, M., Pricoco, E., Piombino, E., Ferro, M., et al. (2021). Insulin Signaling, Androgen Receptor and PSMA Immunohistochemical Analysis by Semi-automated Tissue Microarray in Prostate Cancer with Diabetes (DIAMOND Study). Translational Res. 238, 25–35. doi:10.1016/j.trsl.2021.07.002
Chen, Y.-H., Chen, S.-H., Hou, J., Ke, Z.-B., Wu, Y.-P., Lin, T.-T., et al. (2019). Identifying Hub Genes of clear Cell Renal Cell Carcinoma Associated with the Proportion of Regulatory T Cells by Weighted Gene Co-expression Network Analysis. Aging 11 (21), 9478–9491. doi:10.18632/aging.102397
Cheng, L., Yuan, M., Li, S., Lian, Z., Chen, J., Lin, W., et al. (2021). Identification of an IFN-β-Associated Gene Signature for the Prediction of Overall Survival Among Glioblastoma Patients. Ann. Transl Med. 9 (11), 925. doi:10.21037/atm-21-1986
Claps, M., Mennitto, A., Guadalupi, V., Sepe, P., Stellato, M., Zattarin, E., et al. (2020). Immune-checkpoint Inhibitors and Metastatic Prostate Cancer Therapy: Learning by Making Mistakes. Cancer Treat. Rev. 88, 102057. doi:10.1016/j.ctrv.2020.102057
Gamat, M., and McNeel, D. G. (2017). Androgen Deprivation and Immunotherapy for the Treatment of Prostate Cancer. Endocr. Relat. Cancer 24 (12), T297–t310. doi:10.1530/erc-17-0145
Ghoochani, A., Hsu, E.-C., Aslan, M., Rice, M. A., Nguyen, H. M., Brooks, J. D., et al. (2021). Ferroptosis Inducers Are a Novel Therapeutic Approach for Advanced Prostate Cancer. Cancer Res. 81 (6), 1583–1594. doi:10.1158/0008-5472.Can-20-3477
Huang, J. B., Wu, Y. P., Lin, Y. Z., Cai, H., Chen, S. H., Sun, X. L., et al. (2020). Up‐regulation of LIMK1 Expression in Prostate Cancer Is Correlated with Poor Pathological Features, Lymph Node Metastases and Biochemical Recurrence. J. Cell Mol Med 24 (8), 4698–4706. doi:10.1111/jcmm.15138
Isaacsson Velho, P., and Antonarakis, E. S. (2018). PD-1/PD-L1 Pathway Inhibitors in Advanced Prostate Cancer. Expert Rev. Clin. Pharmacol. 11 (5), 475–486. doi:10.1080/17512433.2018.1464388
Itkonen, H. M., Gorad, S. S., Duveau, D. Y., Martin, S. E. S., Barkovskaya, A., Bathen, T. F., et al. (2016). Inhibition of O-GlcNAc Transferase Activity Reprograms Prostate Cancer Cell Metabolism. Oncotarget 7 (11), 12464–12476. doi:10.18632/oncotarget.7039
Ju, Q., Li, X., Zhang, H., Yan, S., Li, Y., and Zhao, Y. (2020). NFE2L2 Is a Potential Prognostic Biomarker and Is Correlated with Immune Infiltration in Brain Lower Grade Glioma: A Pan-Cancer Analysis. Oxidative Med. Cell Longevity 2020, 3580719. doi:10.1155/2020/3580719
Ke, Z. B., Wu, Y. P., Huang, P., Hou, J., Chen, Y. H., Dong, R. N., et al. (2021). Identification of Novel Genes in Testicular Cancer Microenvironment Based on ESTIMATE Algorithm‐derived Immune Scores. J. Cell Physiol 236 (1), 706–713. doi:10.1002/jcp.29898
Kwon, O.-S., Kwon, E.-J., Kong, H.-J., Choi, J.-Y., Kim, Y.-J., Lee, E.-W., et al. (2020). Systematic Identification of a Nuclear Receptor-Enriched Predictive Signature for Erastin-Induced Ferroptosis. Redox Biol. 37, 101719. doi:10.1016/j.redox.2020.101719
Malik, A., Srinivasan, S., and Batra, J. (2019). A New Era of Prostate Cancer Precision Medicine. Front. Oncol. 9, 1263. doi:10.3389/fonc.2019.01263
Mazzu, Y. Z., Armenia, J., Nandakumar, S., Chakraborty, G., Yoshikawa, Y., Jehane, L. E., et al. (2020). Ribonucleotide Reductase Small Subunit M2 Is a Master Driver of Aggressive Prostate Cancer. Mol. Oncol. 14 (8), 1881–1897. doi:10.1002/1878-0261.12706
Milonas, D., Venclovas, Ž., Gudinaviciene, I., Auskalnis, S., Zviniene, K., Jurkiene, N., et al. (2019). Impact of the 2014 International Society of Urological Pathology Grading System on Concept of High-Risk Prostate Cancer: Comparison of Long-Term Oncological Outcomes in Patients Undergoing Radical Prostatectomy. Front. Oncol. 9, 1272. doi:10.3389/fonc.2019.01272
Mock, A., Geisenberger, C., Orlik, C., Warta, R., Schwager, C., Jungk, C., et al. (2016). LOC283731 Promoter Hypermethylation Prognosticates Survival after Radiochemotherapy in IDH1 Wild-type Glioblastoma Patients. Int. J. Cancer 139 (2), 424–432. doi:10.1002/ijc.30069
Modena, A., Ciccarese, C., Iacovelli, R., Brunelli, M., Montironi, R., Fiorentino, M., et al. (2016). Immune Checkpoint Inhibitors and Prostate Cancer: A New Frontier? Oncol. Rev. 10 (1), 293. doi:10.4081/oncol.2016.293
Nassar, Z. D., Mah, C. Y., Dehairs, J., Burvenich, I. J., Irani, S., Centenera, M. M., et al. (2020). Human DECR1 Is an Androgen-Repressed Survival Factor that Regulates PUFA Oxidation to Protect Prostate Tumor Cells from Ferroptosis. Elife 9, e54166. doi:10.7554/eLife.54166
Riley, R. S., June, C. H., Langer, R., and Mitchell, M. J. (2019). Delivery Technologies for Cancer Immunotherapy. Nat. Rev. Drug Discov. 18 (3), 175–196. doi:10.1038/s41573-018-0006-z
Sircar, K., Huang, H., Hu, L., Cogdell, D., Dhillon, J., Tzelepi, V., et al. (2012). Integrative Molecular Profiling Reveals Asparagine Synthetase Is a Target in Castration-Resistant Prostate Cancer. Am. J. Pathol. 180 (3), 895–903. doi:10.1016/j.ajpath.2011.11.030
Stockwell, B. R., Friedmann Angeli, J. P., Bayir, H., Bush, A. I., Conrad, M., Dixon, S. J., et al. (2017). Ferroptosis: A Regulated Cell Death Nexus Linking Metabolism, Redox Biology, and Disease. Cell 171 (2), 273–285. doi:10.1016/j.cell.2017.09.021
Tousignant, K. D., Rockstroh, A., Poad, B. L. J., Talebi, A., Young, R. S. E., Taherian Fard, A., et al. (2020). Therapy-induced Lipid Uptake and Remodeling Underpin Ferroptosis Hypersensitivity in Prostate Cancer. Cancer Metab. 8, 11. doi:10.1186/s40170-020-00217-6
Venclovas, Z., Jievaltas, M., and Milonas, D. (2019). Significance of Time until PSA Recurrence after Radical Prostatectomy without Neo- or Adjuvant Treatment to Clinical Progression and Cancer-Related Death in High-Risk Prostate Cancer Patients. Front. Oncol. 9, 1286. doi:10.3389/fonc.2019.01286
Vitkin, N., Nersesian, S., Siemens, D. R., and Koti, M. (2019). The Tumor Immune Contexture of Prostate Cancer. Front. Immunol. 10, 603. doi:10.3389/fimmu.2019.00603
Wei, C., Zhang, Y., Malik, H., Zhang, X., Alqahtani, S., Upreti, D., et al. (2019). Prediction of Postprostatectomy Biochemical Recurrence Using Quantitative Ultrasound Shear Wave Elastography Imaging. Front. Oncol. 9, 572. doi:10.3389/fonc.2019.00572
Wong, R. L., and Yu, E. Y. (2021). Refining Immuno-Oncology Approaches in Metastatic Prostate Cancer: Transcending Current Limitations. Curr. Treat. Options. Oncol. 22 (2), 13. doi:10.1007/s11864-020-00808-x
Xu, N., Chen, H.-J., Chen, S.-H., Xue, X.-Y., Chen, H., Zheng, Q.-S., et al. (2016). Reduced Connexin 43 Expression Is Associated with Tumor Malignant Behaviors and Biochemical Recurrence-free Survival of Prostate Cancer. Oncotarget 7 (41), 67476–67484. doi:10.18632/oncotarget.11231
Xu, N., Chen, S. H., Lin, T. T., Cai, H., Ke, Z. B., Dong, R. N., et al. (2020). Development and Validation of Hub Genes for Lymph Node Metastasis in Patients with Prostate Cancer. J. Cell Mol Med 24 (8), 4402–4414. doi:10.1111/jcmm.15098
Xu, N., Ke, Z.-B., Chen, Y.-H., Wu, Y.-P., Chen, S.-H., Wei, Y., et al. (2020). Risk Factors for Pathologically Confirmed Lymph Nodes Metastasis in Patients with Clinical T2N0M0 Stage Prostate Cancer. Front. Oncol. 10, 1547. doi:10.3389/fonc.2020.01547
Zhang, X., Du, L., Qiao, Y., Zhang, X., Zheng, W., Wu, Q., et al. (2019). Ferroptosis Is Governed by Differential Regulation of Transcription in Liver Cancer. Redox Biol. 24, 101211. doi:10.1016/j.redox.2019.101211
Keywords: prostate cancer, ferroptosis, biochemical recurrence, tumor immune microenvironment, molecular signature
Citation: Ke Z-B, You Q, Sun J-B, Zhu J-M, Li X-D, Chen D-N, Su L, Zheng Q-S, Wei Y, Xue X-Y and Xu N (2022) A Novel Ferroptosis-Based Molecular Signature Associated with Biochemical Recurrence-Free Survival and Tumor Immune Microenvironment of Prostate Cancer. Front. Cell Dev. Biol. 9:774625. doi: 10.3389/fcell.2021.774625
Received: 12 September 2021; Accepted: 02 December 2021;
Published: 06 January 2022.
Edited by:
Markus Ritter, Paracelsus Medical University, AustriaReviewed by:
Hailiang Hu, Southern University of Science and Technology, ChinaDimitrios Korentzelos, University of Pittsburgh Medical Center, United States
Copyright © 2022 Ke, You, Sun, Zhu, Li, Chen, Su, Zheng, Wei, Xue 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.
*Correspondence: Ning Xu, drxun@fjmu.edu.cn; Xue-Yi Xue, xuexueyi@fjmu.edu.cn
†These authors have contributed equally to this work