Development and validation of novel inflammatory response-related gene signature to predict prostate cancer recurrence and response to immune checkpoint therapy

1.   Introduction

Prostate cancer (PCa), a highly malignant tumor of the urinary tract, is an important public health problem. 2020 saw the addition of 1,414,000 new cases and the death of 375,000 people worldwide ^[1]. In China, PCa has become the ninth most prevalent malignancy in 2022 ^[2]. It was estimated that more than 260,000 people in the United States would develop into PCa in 2022 ^[3]. The current common treatments for PCa include radical prostatectomy (RP), endocrinotherapy, radiotherapy, chemotherapy, and immunotherapy. Nevertheless, more than 30% of PCa patients with post-RP would undergo biochemical recurrence (BCR), which could facilitate the development of castration-resistant prostate cancer (CRPC), eventually progressing into distant metastases or PCa-specific mortality ^[4]. The existing clinical characteristics were not adequate for clinicians to determine the prognosis of the increasing number of PCa patients. In light of this, it was imperative to develop more accurate methods of predicting prognosis and treatment outcomes.

FDA-approved immune checkpoint therapies (ICTs) have been determined as the routine treatment of locally advanced and other metastatic cancers such as melanoma, lung cancer, bladder cancer ^[5,6,7]. Targeting immune checkpoints such as programmed cell death ligand 1 (PD-L1), programmed cell death protein 1 (PD-1), and cytotoxic T-lymphocyte-associated protein 4 (CTLA4) to boost anti-tumor immunity was currently one of the most prospective therapeutic regimens in immunotherapy ^[8]. However, multiple preclinical studies undertaken to harness its benefit for PCa have met with limited success ^[9,10]. Meanwhile, clinically useful biomarkers to guide treatment selection remained unclear, restricting patients from getting better treatments. ICTs were promising strategies in cancer treatments that were strongly influenced by immune microenvironmental factors ^[11]. Novel treatment strategies utilizing novel molecular targets, drugs, and combination of immunotherapy exhibited an enormous potential in PCa treatment, especially for CRPC ^[12,13]. It was well known that when the immune microenvironment was deficient in the immune infiltrates of CD8+ T cell, anti-PD-1 treatment results would not be encouraging, while patients would receive favorable results when such a high infiltration was present ^[14]. Identifying a reliable biomarker involved with the tumor immune microenvironment (TIME) for predicting the therapeutic sensitivity and effectiveness of ICTs appeared to be imminent.

A growing body of literatures reported that inflammation played a crucial role in the occurrence, development, and metastasis of PCa ^[15,16]. A significant host defense mechanism, inflammation occurred when body tissues were confronted with harmful stimuli such as irritants, pathogens, or damaged cells ^[17]. Nevertheless, it has also been revealed that cancer development was linked to chronic inflammation through the activation of a range of inflammatory cytokines and signals. An increased likelihood in PCa risk was paralleled with chronic inflammation, be it local or systemic. Additionally, carcinogenic alteration might induce chronic inflammation in the microenvironment, and had a variety of cancer-promoting effects in PCa, including angiogenesis, proliferation, invasion, and distant metastasis ^[18]. Hence, a comprehensive insight of inflammatory response-related genes (IRRGs) in PCa could enhance our understanding of the tumor immune status and provide novel and improved immunotherapeutic strategies for patients.

In this research, a unique genomic signature involved in IRRGs that could independently predict the progression and recurrence of PCa patients was determined. Superior to conventional clinical features, the signature was closely correlated with clinicopathological characteristics. Then a nomogram was constructed to apply in clinical practices. Moreover, comprehensive studies of TIME were conducted. The signature played an essential role in predicting response to immunotherapy. In all, this was an innovative approach to survival prediction and clinical treatments of PCa patients.

2.   Materials and methods

2.1.   Multi-omics data acquisition of PCa patients

Firstly, the relevant clinical information and RNA sequencing transcriptome with HTSeq-FPKM format, containing 499 PCa and 52 corresponding adjacent normal samples, were extracted from the Cancer Genome Atlas database (TCGA-PRAD, https://portal.gdc.cancer.gov/projects/TCGA-PRAD). The detail information of samples' IDs was attached to our supplementary file (Table S1). All ensembl gene IDs were transferred to gene symbols using an annotation GTF file for the following analysis. Clinicopathological data included age, grade, AJCC-TNM stage, BCR-free time, and BCR status. After excluding samples without BCR events or recorded BCR time, 427 PCa samples were screened out for the following investigation. Moreover, tumor mutation burden (TMB) data of 475 PCa samples in Mutation Annotation Format was retrieved from the TCGA database and was implemented through the maftools R package ^[19]. Meanwhile, we downloaded copy number variation (CNV) data, including 497 PCa patients, from the TCGA database, which was analyzed based on Perl scripts and the RCircos R package ^[20]. The TCGA-PRAD dataset was appointed as the training cohort to build the signature. 200 IRRGs were obtained from hallmark gene sets in the gene set enrichment analysis (GSEA) Molecular Signatures Database (MSigDB, http://www.gsea-msigdb.org/gsea/msigdb/index.jsp). We extracted the GSE21034 cohort (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE21034), including genomic profile of 140 PCa patients (Table S2), from the Gene Expression Omnibus (GEO) and the DKFZ cohort (https://www.cbioportal.org/study/summary?id=prostate_dkfz_2018) of 82 PCa samples (Table S3) from the cbioportal website to validate the signature. Additionally, we distinguish the mRNAs and lncRNAs from the RNA sequencing profiles, including 4 of BPH-1 (the benign prostate cell line) and 4 of DU145 (the PCa cell line). The mRNAs were applied to evaluate the differentially expressed genes (DEGs) of our signature. The total RNA was isolated and obtained from above two cell lines using TRIzol and chloroform. On the whole, we removed growth media and added 1.0 mL of TRIzol reagent directly to the culture dish to lyse a total of 3 × 10⁶ cells. Then, we pipette the lysate up and down several times to homogenize and incubate for 5 minutes. 0.2 mL of chloroform was added for lysis and incubating for 2–3 minutes. Finally, we centrifuged the sample for 15 minutes at 12,000 × g at 4 ℃ and then transferred the aqueous phase containing the RNA to a new tube for the following study. The entire process of RNA sequencing was carried out by Sinotech Genomics (China) and the RNA seq data presented in our study has been deposited in GEO (GSE210205).

2.2.   Screening of differentially expressed IRRGs

The limma package was utilized to single out differentially expressed IRRGs with a threshold of P value < 0.05 in the TCGA-PRAD dataset and RNA sequencing profile respectively ^[21]. Venn diagram was performed to determine the upregulated or downregulated IRRGs between the above-mentioned two datasets. In addition, the position changes of the IRRGs on the human chromosomes were depicted by R package RCircos ^[20]. Furthermore, the correlation network and protein-protein interaction (PPI) were constructed using the igraph and reshape2 R packages and the STRING website separately.

2.3.   Establishment and evaluation of an IRRGs signature

We determined the IRRGs with prognostic value through Cox regression analysis and developed an IRRGs signature using the least absolute shrinkage and selection operator (Lasso) regression test.

PCa patients were further stratified into the high- and low-risk groups, with the median risk score as the threshold. The receiver operating characteristic (ROC) curve was employed to assess the area under the curve (AUC) value and the effectiveness of prediction based on survminer, survival, and timeROC R packages. The independence of the signature was confirmed by univariate and multivariate analysis through survival R package.

2.4.   Distribution analysis, a predictive nomogram building, and functional enrichment analysis

Principal component analysis (PCA) and t-distributed stochastic neighbour embedding (t-SNE) were illustrated by ggplot2, scatterplot3d, and Rtsne R packages. The rms R package was adopted to set up a nomogram integrating risk score and diverse clinical variables for predicting BCR-free survival ^[22]. The calibration curves were depicted to evaluate the accuracy of the nomogram. We also assessed the superiority of our signature compared to other prognostic models constructed by the researchers using concordance index (C-index). Additionally, GO and KEGG enrichment analysis was analyzed by org.Hs.eg.db, enrichplot, ggplot2, and clusterProfiler R packages ^[23,24,25].

2.5.   Immune landscape of our signature

In order to evaluate the immune status in both groups, including estimate score, stromal score, immune score, and corresponding tumor purity, the ESTIMATE algorithm was performed by limma and ggpubr R packages. We further evaluated the immune cell infiltration levels of the core IRRGs based on TIMER database (https://cistrome.shinyapps.io/timer/) ^[26]. A differential immune cell abundance of hub IRRGs in PCa was determined using the QUANTISEQ, XCELL, CIBERSORT, CIBERSORT-ABS, TIMER, EPIC, and MCPCOUNTER algorithms ^[27]. The analysis of Wilcoxon test revealed the difference between the high-risk and low-risk groups regarding immune infiltration. The lollipop graph depicted the connection of the risk score with immune cell subtypes through Spearman correlation analysis. The discrepancy of immune landscape, comprising of immune function and immune cell subtypes, between both subgroups was assessed by the single-sample gene set enrichment approach (ssGSEA). Additionally, the tumor immune dysfunction and exclusion (TIDE) was computed by the online tool (http://tide.dfci.harvard.edu/).

2.6.   Statistical analysis

A series of above-mentioned R packages were carried out to conduct statistical analyses using R version 4.1.1. Continuous variables were analyzed by Student t-test or Wilcoxon test. Categorical variables, including grade and AJCC-TNM stage, were calculated and classified as 1 (grade 1 or stage Ⅰ), 2 (grade 2 or stage Ⅱ), 3 (grade 3 or stage Ⅲ), and 4 (grade 4 or stage Ⅳ) in Cox regression and nomogram analysis. The P-value < 0.05 indicated the statistical significance.

3.   Results

3.1.   Identification of differentially expressed IRRGs in PCa

As shown in the workflow of Figure 1, we first obtained 114 DEGs via investigating the expression patterns of 200 IRRGs in PCa and normal tissues in the TCGA database (Figure 2A), and screened out 8663 DEGs between BPH-1 and DU145 cell lines by RNA sequencing (Figure 2B). A total of 41 intersection genes between TCGA-PRAD and RNA-seq profiles, containing 10 upregulated and 31 downregulated DEGs (Tables S4 and S5), were further displayed in the Venn diagram (Figure 2C). The frequency of CNV indicated that most DEGs appeared as deletion copy numbers (Figure 2D). The CNV alteration sites of DEGs were annotated on chromosomes (Figure 2E). The correlation and PPI network of these DEGs were presented in the Figure 3A, B. As a result, a list of DEGs were identified as the hub genes, including RELA, IL1A, and IL1B, the prognostic values of which were validated in PCa ^[28,29]. Functional enrichment analyses indicated that these DEGs were mainly involved in regulation of inflammatory response, NF-κB signaling pathway, cytokine production, and PD-L1 expression and PD-1 checkpoint pathway in cancer (Figure 3C, D).

3.2.   Construction of an IRRGs signature in PCa

To obtain genes with prognostic significance from IRRGs, we determined five IRRGs associated with prognosis for the next investigation (Figure 4A). Lasso regression is widely used to eliminate multi-collinearity. After feature selection by Lasso regression and backward elimination regression, correlations within five prognostic IRRGs were dramatically reduced. In our risk signature, the coefficient of LPAR1 was 0 so that LPAR1 was excluded from the model. The formula to calculate the signature was as follow: Riskscore = -0.347 × BTG2 - 0.183 × CX3CL1 + 0.646 × KCNJ2 + 1.171 × RELA (Figure 4B, C). We stratified the patients into two groups (ie, high-risk and low-risk). It was found that there was an apparent difference in BCR-free survival between the two groups where PCa patients' mortality gradually increased with rising risk scores (Figure 4D). PCa patients in the high-risk group conferred undesirable BCR-free survival outcomes (Figure 4E). Furthermore, time-dependent ROC analysis illustrated the specificity and sensitivity of our signature in BCR-free survival prediction and the AUC values were 0.686, 0.692, and 0.674 at 1-, 3-, and 5-years respectively (Figure 4F). The results of univariate and multivariate analysis revealed our signature could act as an independent prognostic predictor for PCa patients (Figure 4G, H).

3.3.   Verification of the signature in survival prediction

In addition, to confirm the predictive ability of the signature, we verified it in two independent external datasets. The corresponding results in the GSE21034 cohort were similar to those in the TCGA-PRAD cohort (Figure 5A). Analogously, the results of ROC, survival, and Cox regression analysis were consistent with those of the TCGA-PRAD cohort, indicating that the high-risk patients represented poorer prognosis and the validity and generalizability of our signature (Figure 5B, C, E). Additionally, the DKFZ cohort further validated our signature (Figure 5D) and the random effects model of these three cohorts based on the meta-analysis further confirmed the role of the signature as an independent prognostic factor in PCa patients (Figure 5F).

3.4.   Correlation between the signature and disease progression

We have demonstrated the ability of the signature to predict BCR-free survival in PCa, but the connection of the signature to several clinicopathologic characteristics in the course of the disease remained to be elucidated. The results demonstrated that there were significant differences in grade, TNM stage, and survival state between the two subgroups (Figure 6A–D). Further exploring the relationship between the signature and clinical features, we found that as risk score increased, the progression was more like to develop into high-grade or advanced-stage PCa (Figure 6E). Collectively, the signature had a crucial impact on disease progression and was closely related to PCa prognosis.

3.5.   Distribution patterns and construction of a nomogram

t-SNE and PCA analyses showed that both subgroups were distributed in opposite directions (Figure 7A, B). Nevertheless, the IRRGs expression and the genome-wide expression profiles didn't display a clear separation (Figure 7C, D). To facilitate the assessment of prognosis in PCa patients, we constructed a nomogram integrating multiple clinical variables to quantitatively predict BCR-free survival (Figure 7E). A higher point represented a higher mortality rate for PCa patients. In addition, there was a steady consistency between nomogram-predicted probability and actual observations regarding the 1-, 2-, and 3-years BCR-free survival (Figure 7F). Furthermore, the C-index was higher than those of any other five published models ^{[30,31,32,33,34]} (Figure 7G, H). Therefore, our results indicated that the signature could serve as an effective and reliable indicator for predicting the prognosis of PCa in clinical applications.

3.6.   Functional analyses and TMB

We found that the signature was associated with many biological processes of immunity including neutrophil activation, immune response, and neutrophil mediated immunity by GO enrichment analysis (Figure 8A). Additionally, the result of KEGG pathway indicated that our signature was enriched in TNF signaling pathway, which were related to inflammatory response (Figure 8B). Next, we explored gene mutations to further obtain biological perspective into immunological characteristic of two risk groups. The waterfall plot revealed that missense variations were more frequently than other mutation types and higher mutation counts occurred in the high risk group. We then determined top 20 mutated genes in both subgroups (Figure 8C, D).

3.7.   Analysis of immune status in PCa

We used the ESTIMATE algorithm to find that the high-risk group possessed lower levels of infiltrating stromal cells and immune cells (Figure 9A–E). Meanwhile, the four genes that constructed the signature were also robustly related to a series of immune cell subtypes (Figure 10A–D). To further elucidate the impact of our signature on TIME, a thorough examination of immune cell infiltration in PCa was conducted. The heatmap further depicted the differences in the distribution of immune cell subtypes between both risk subgroups (Figure 11A). Multiple algorithms were applied to analyze the infiltration patterns of immune cell subtypes, the results of which indicated that the signature had a close relationship with a range of immune cell subtypes (Figure 11B). The low-risk group exhibited a significantly higher infiltrated abundance of neutrophils, macrophages, regulatory cells (Tregs), B cells, and CD8+ T cells compared to the high-risk group (Figure 11C). The results also revealed that significant differences were observed in promoting inflammation, parainflammation, human leukocyte antigen (HLA), cytolytic activity, and checkpoint between two subgroups (Figure 11C).

3.8.   The signature in immunotherapy

We examined the expression of ICT-related genes in two subgroups to understand the impact of our signature on immunotherapy. The expression of immune checkpoints such as PD-1, CTLA-4, and PD-L1 was higher in the low-risk group (Figure 12A). Moreover, analysis of HLAs gene expression disclosed that the low-risk group exhibited higher expression levels (Figure 12B). Next, The TIDE algorithm was adopted to determine the likelihood of immunotherapy response, and the results indicated that low-risk patients were suitable for ICTs (Figure 12C). PCa patients in the low-risk group might make a sensitive treatment response in ICTs applications and they should receive more attention from clinicians.

3.9.   Verification of the expression levels of four IRRGs

We performed RNA sequencing on BPH-1 (the benign prostate cell line) and DU145 (the PCa cell line). The results indicated higher expression of BTG2, CX3CL1 and KCNJ2 in BPH-1 and higher expression of RELA in DU145 (Figure 13A). The mRNA expression in normal samples and tumor tissues from TCGA-PRAD also showed the same results as RNA sequencing (Figure 13B). The increased expression of BTG2 predicted a favorable survival outcome in PCa patients with post-RP ^[35]. The overexpression of CX3CL1 enhanced spinal metastasis via activating the Src/FAK pathway in PCa ^[36]. Interference with KCNJ2 played a vital role in the proliferation and invasion in multiple malignancies ^[37,38]. Additionally, RELA was an independent predictor of BCR in PCa ^[28].

4.   Discussion

The currently available treatments could not meet the increasing number of PCa patients who needed better outcomes. About 34% of patients with localized disease at initial diagnosis eventually suffered BCR at ten years after RP without neoadjuvant or adjuvant therapy ^[39]. The median survival time from BCR to inoperable locally advanced or metastatic stage and from metastasis to death was approximately 8 and 5 years, separately ^[40]. However, the most clinically used TMN staging system primarily relied on anatomical information regardless of any biological information, which was out of step with the needs for clinical application. It was of crucial significance to predict the patients' survival for guiding the treatments, necessitating further insights to ascertain a biomarker that could accurately predict PCa patients' prognosis.

With the advances of high-throughput methodologies like microarrays and RNA sequencing, gene expression profiling has become an effective and practical tool for the discovery of molecular biomarkers associated with the prognosis of malignancies ^[41]. In our research, we have constructed an IRRGs signature including four IRRGs in the TCGA-PRAD cohort. As PCa patients were separated into the high- and low-risk groups by the intermediate risk score, the high-risk patients were obviously at poorer clinical prognosis. We determined the signature was an independent prognostic predictor for PCa patients. Subsequently, the results were well verified in two independent external datasets and the statistical power of our study was improved by a meta-analysis. Additionally, we found that the signature exhibited excellent prediction accuracy compared to five published models. Clinicopathological factors were clinically used to diagnose and differentiate PCa. However, the outcomes of some patients with similar clinical characteristics were totally different. Hence, we established a nomogram to explore the effect of the signature on patients with the same clinical features. What excited us was that our signature could still accurately distinguish the different prognoses of patients with the same clinical characteristics. The nomogram could help clinicians better formulate treatment strategies based on the patient's condition and help personalize treatments. Finally, we analyzed the differential expression between BPH-1 and DU145 cell lines through RNA sequencing, and the expression of four core IRRGs constituting the signature differed significantly in cancer and adjacent non-tumor tissues from TCGA, finding that BTG2, CX3CL1, and KCNJ2 were lowly expressed in cancer, while the expression of RELA was exactly the opposite. Collectively, the potential of our signature to predict survival of PCa patients independently and correctly at various clinical stages was comprehensively and robustly elucidated, making it possible for patients to receive better treatments.

The chronic inflammation promoted the development of PCa and contributed to prostate carcinogenesis, the microenvironment of which was composed of immune cells, stromal cells, and tumor cells ^[18]. The complex and diverse interactions between tumor cells and immune cells influenced tumor development and response to therapeutics ^[14,42]. The results revealed that our signature was enriched in neutrophil activation, immune response, and neutrophil mediated immunity. Additionally, immune function analysis further identified our signature was related to inflammation-promoting, parainflammation, and other immune processes or cell subtypes. Therefore, to further explore the effects of immune cells on the immunotherapy of PCa, a variety of algorithms were applied to explore the differences in immune cells between two subgroups. CD8+ T cells as major effector immune cells played a vital role in antitumor effects and were critical for tumorigenesis and progression ^[43,44].

CD8+ T cell was considered to act as an anti-tumor factor in most conditions ^[45,46], the cut down of which might promote BCR in PCa ^[47]. High abundance of CD8+ T cells in the TIME was associated with improved survival or prolonged BCR in PCa patients ^[47,48]. Additionally, the activation of CD8+ T cell could induce prostate inflammation ^[49]. Here, a close relationship between the prognosis of PCa and immune cells was revealed. This was paralleled with our result that the high-risk patients conferred a worse prognosis along with lower infiltrated levels of CD8+ T cells. Our signature not only stratified patients with different prognosis, but also effectively analyzed the status of immune infiltration cells in different PCa patients, facilitating the choice of ICTs for clinicians.

With the introduction of ICTs into the treatment of PCa, ICTs had dramatically changed the treatment landscape of metastatic CRPC patients and the scenario of immunotherapy was promising ^[50]. Detecting patients who could benefit most from ICTs and the combination of other treatment modalities were essential, given the relatively limited success and low overall response rates ^[8,51]. Some studies had found that anti-PD-1 therapy was broadly taken into clinical practice in PCa ^[52,53]. The expression of PD-L1 was related to the infiltrated level of CD8+ T cell and aggressive clinical progression or BCR whether in localized PCa or metastatic disease ^[54,55,56]. PD-L1 was not only an an independent prognostic predictor for PCa patients but also achieved higher treatment response to ICTs ^[57,58]. Increased immune checkpoint gene suppressed the anti-tumor immune response of T cells by increasing the expression of PD-1, CTLA4, and PD-L1 reporters. With the help of the TIDE, as a predictor of ICTs efficacy, the low-risk patients would be likely to escape from tumor immune surveillance. Our study showed that patients in the low-risk group had higher expression of PD-1, CTLA4, and PD-L1, suggesting that the low-risk patients might have a favorable response to anti-PD-1/CTLA4/PD-L1 therapy. Considering that HLA-related genes played an essential role in regulating the immune response, we compared the expression of HLA-related genes between two risk groups. The immune score and expression levels of HLA-related genes in the low-risk group were significantly higher than those in the high-risk group, while the tumor purity exhibited the opposite trend, likely explaining why the low-risk group patients were more suitable for immunotherapy. Collectively, our signature offered innovative insights into improving the clinical reaction of PCa patients to ICTs, which helped patients get more suitable and accurate treatments.

However, there were still some deficiencies in our study. Firstly, more prostate cell line quantities and types were further needed to verify our conclusions in the future. Besides, immunohistochemical staining of clinical samples from our hospital would verify the expression levels of four hub genes in PCa on the other hand.

5.   Conclusions

In summary, our signature was an effective and promising IRRGs prognostic biomarker for PCa patients. The widespread regulatory mechanism of the signature in the immune landscape and immunotherapy of PCa was comprehensively elucidated. We provided a novel insight into survival prediction and immunotherapy in PCa.

Acknowledgments

This study was supported by Guangdong Medical Science and Technology Research Fund Project (B2022332, B2020127) granted to Qingbiao Chen, Foshan Science and Technology Innovation Project (2020001005794) granted to Jingbo Lin, National Natural Science Foundation of China (82072813, 81571427) and Guangzhou Municipal Science and Technology Project (201803040001) granted to Weide Zhong.

Conflict of interest

The authors declare no competing interests.