Identification of a Costimulatory Molecule-Related Signature for Predicting Prognostic Risk in Prostate Cancer

Abstract

Costimulatory molecules have been proven to enhance antitumor immune responses, but their roles in prostate cancer (PCa) remain unexplored. In this study, we aimed to explore the gene expression profiles of costimulatory molecule genes in PCa and construct a prognostic signature to improve treatment decision making and clinical outcomes. Five prognosis-related costimulatory molecule genes (RELT, TNFRSF25, EDA2R, TNFSF18, and TNFSF10) were identified, and a prognostic signature was constructed based on these five genes. This signature was an independent prognostic factor according to multivariate Cox regression analysis; it could stratify PCa patients into two subgroups with different prognoses and was highly associated with clinical features. The prognostic significance of the signature was well validated in four different independent external datasets. Moreover, patients identified as high risk based on our prognostic signature exhibited a high mutation frequency, a high level of immune cell infiltration and an immunosuppressive microenvironment. Therefore, our signature could provide clinicians with prognosis predictions and help guide treatment for PCa patients.

Introduction

Over the past few years, prostate cancer (PCa) has become one of the most common lethal malignant tumors in men, posing a grave danger to human health (Nguyen-Nielsen and Borre, 2016). The latest global cancer statistics predicted approximately 191,930 new cases of PCa and 33,330 PCa-related deaths in America in 2020 (Siegel et al., 2020). Early stage PCa has a good prognosis after surgical resection, but once it recurs or develops into metastatic castration-resistant prostate cancer (mCRPC), there are fewer treatment options available, and the average survival time is only 2–3 years (Omlin et al., 2013; Wang et al., 2018). With advancements in medicine, various targeted therapies and immunotherapies have further improved the prognosis of PCa, especially metastatic PCa. However, only a small percentage of PCa patients can benefit from these therapies (Sharma et al., 2017). Epigenetic regulation of gene plays a critical role in cancer evolution (Grasso et al., 2012), and gene expression variations or mutations have been frequently observed in PCa (Gu et al., 2015). Therefore, exploring new genetic and epigenetic biomarkers will be useful for predicting specific survival of PCa patients and improving clinical outcomes.

Immunotherapy has become an important cancer treatment method that can activate the immune system to attack and kill tumor cells. In particular, immune checkpoint inhibitor (ICI) therapy has seen great success in multiple types of cancer (Loos et al., 2008; Jansen et al., 2019). However, only a minority of PCa patients have an obvious response to immunotherapy, stressing the need to explore more effective treatment methods (Sharma et al., 2017). Studies have shown that the tumor microenvironment (TME), which includes immune cells, stromal cells, mesenchymal cells, cytokines, chemokines, and blood vessels, plays an important role in the process of tumorigenesis and development (Casey et al., 2015). Thus, a deeper understanding of the immune microenvironment will help us improve PCa patient outcomes. In fact, many studies are currently exploring the therapeutic potential of costimulatory molecules in cancer (Wei et al., 2018). Costimulatory molecules play an important role in the regulation of tumor immunity by affecting the activation, proliferation, survival and secretion of T cells (Croft et al., 2013; Schildberg et al., 2016). At present, there are two main families of costimulatory molecules: the B7-CD28 family and the tumor necrosis factor (TNF) family. The most common ICIs target programmed cell death protein 1 (PD-1) and programmed cell death 1 ligand 1 (PD-L1), which are CD28 and B7 family members, respectively (Zhang et al., 2018). The latest research showed that TNF family members, including 19 TNF ligand superfamily (TNFSF) members and 29 TNF receptor superfamily (TNFRSF) members, are also potential therapeutic targets that play important roles in the regulation of cellular functions. TNFSF/TNFRSF members can control immune cells to coordinate various mechanisms that drive the costimulation or cosuppression of immune responses (Dostert et al., 2019). However, the molecular functions of these costimulatory molecules in PCa remain unexplored. Therefore, we aimed to develop a specific and effective prognostic signature based on several costimulatory molecules to guide treatment and improve the clinical outcomes of PCa.

High-throughput sequencing and bioinformatics technologies have facilitated the identification of key pivotal biomarkers for PCa (Wang et al., 2013; Gu et al., 2020). Here, we acquired 60 costimulatory molecules (comprising 12 B7-CD28 family members and 48 TNF family members) as candidates for potential molecular therapy targets (Janakiram et al., 2015; Dostert et al., 2019). Subsequently, the mRNA sequencing data and clinical information of PCa patients were downloaded from The Cancer Genome Atlas (TCGA) database to predict disease free survival, and a prognostic signature for PCa patients was constructed. We validated the efficiency of the prognostic signature in four external datasets. Previous studies indicated that the composition of infiltrating immune cells in the tumor microenvironment is fundamental to the outcomes of immunotherapy (Wu et al., 2020). Thus, understanding immune infiltration in the PCa microenvironment is an essential prerequisite for the success of immunotherapy. In our study, we systematically described the landscape of costimulatory molecules and classified PCa patients into different groups according to the prognostic signature. We further evaluated the immune microenvironment of PCa patients to identify patients who may benefit more from immunotherapy.

Materials and Methods

Data Acquisition and Preprocessing

We identified several pivotal prognosis-related costimulatory molecule genes in PCa using comprehensive bioinformatics analysis. The specific flow chart is shown in Figure 1. The mRNA expression profiles of PCa patients in the TCGA database were obtained from UCSC Xena¹. A total of 52 normal samples and 498 PCa samples were obtained. Disease-free survival information was obtained from the cBioPortal for Cancer Genomics². All expression data were normalized using the “RNA-Seq by Expectation-Maximization” package, and log₂(x + 1) transformation was performed. To select genes with prognostic values and establish a risk signature, samples from PCa patients with unknown survival status or without follow-up information were excluded, and a total of 491 PCa samples were included. Then, we obtained four normalized independent microarray datasets, namely, GSE21034, GSE54460, GSE70768, and GSE70769, from the Gene Expression Omnibus (GEO³) database. The GSE21034, GSE54460, GSE70768, and GSE70769 datasets included 140 PCa samples, 90 PCa samples, 111 PCa samples, and 92 PCa samples, respectively. For the GSE21034 dataset, the expression data were normalized according to the robust multi-array average (RMA) method, and the probes were annotated using the corresponding annotation file. For the GSE54460 dataset, the expression data were expressed as fragments per kilobase per million values. The GSE70768 and GSE70769 datasets were sequenced using the Illumina HumanHT-12 V4.0 expression BeadChip platform, and the probes were annotated using the corresponding “illuminaHumanv4.db” R package. The expression data in these two datasets were log2 transformed and quantile normalized. The detailed information for the datasets used in present study were showed in Supplementary Table 1.

FIGURE 1

Identification of Costimulatory Molecules With Prognostic Significance in PCa

These costimulatory molecule genes were mapped to the TCGA dataset, and univariate Cox regression analysis was performed to identify costimulatory molecule genes with prognostic significance in PCa with P < 0.05. We also used Kaplan–Meier curves and log-rank tests to verify the prognostic values of these survival-related costimulatory molecule genes using the “survival” R package (Tibshirani, 1997). Least absolute shrinkage and selection operator (LASSO) Cox regression analysis was used to screen the most significant prognostic costimulatory molecule genes in PCa using the “glmnet” R package. The optimal values of penalty parameters were determined by 10-fold cross-validation.

Consensus Clustering of Prognosis-Related Costimulatory Molecule Genes

To further explore the functions and prognostic values of the costimulatory molecule genes in PCa from the results of LASSO Cox regression analysis, we carried out consensus clustering to confirm the cluster numbers using the “Consensus ClusterPlus” R package. Then, Kaplan-Meier curves were plotted to verify the prognostic values of the clusters. In addition, gene set enrichment analysis (GSEA) was also performed to reveal the potential functional mechanisms using the h.all.v7.2.symbols.gmt file. False discovery rate (FDR) < 0.25 and normalized P-value < 0.05 were set as the threshold values for significance.

Construction of a Prognostic Signature Based on the Five Survival-Related Costimulatory Molecule Genes

We performed multivariate Cox proportional hazards regression analysis of these five survival-related costimulatory molecule genes to obtain the corresponding coefficients. Subsequently, a prognostic signature was constructed, and the risk score was calculated using the following formula:

β and Exp represent the coefficients from the multivariate Cox proportional hazards regression analysis and the expression levels of selected genes, respectively. In addition, we generated Kaplan–Meier curves and receiver operating characteristic (ROC) curves to evaluate and validate the efficiency of the signature. Then, we performed principal component analysis (PCA) to estimate the distribution patterns and confirm the cluster numbers in the TCGA dataset by using the “ggplot2” R package.

Functional Enrichment Analysis

To explore signature-related biological pathways, genes that were strongly correlated with the risk score (correlation coefficient R > 0.4) were obtained. A total of 525 positively correlated genes and 87 negatively correlated genes were generated. We performed gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses to investigate the possible molecular mechanisms of the risk signature genes using the Database for Annotation, Visualization and Integrated Discovery version 6.8⁴. P < 0.05 was regarded as the cutoff value.

Estimation of the Immune Microenvironment Composition

For quantification of the cellular composition of the immune infiltrates in each risk group, a set of metagenes, including non-overlapping sets of genes that are representative of twenty-eight specific immune cell subpopulations, was obtained from a previous study (Charoentong et al., 2017). Then, we employed single-sample gene set enrichment analysis (ssGSEA) to quantify the 28 types of immune cells based on the metagenes. In the tumor microenvironment, immune and stromal cells are the two main non-tumor components and have been proposed to be valuable for tumor treatment and prognostication. To assess the tumor microenvironment of different risk groups, the immune and stromal scores for the TCGA dataset, reflecting the infiltration levels of non-tumor cells, were calculated using the ESTIMATE package (Yoshihara et al., 2013). The differences in immune cell composition and immune and stromal scores were compared between high-risk and low-risk PCa patients.

Comparison of Significantly Mutated Genes

Tumor mutation burden (TMB) is defined as the total amount of coding errors of somatic genes, base substitutions, insertions or deletions detected per million bases (Chalmers et al., 2017). In the present study, we used 38 MB as the length of exons. TMB was calculated as the number of variants/the length of the exons for each PCa sample via Perl scripts based on the JAVA8 platform (Abida et al., 2019). The somatic mutation status data of PCa samples (workflow type: VarScan2 Variant Aggregation and Masking) were downloaded from the TCGA data portal⁵ in December 2020. Mutation data were filtered using the “maftools” R package and compared between high-risk and low-risk patients.

Statistical Analyses

We performed the t-test or Wilcoxon test for comparisons between the two groups, and used the one-way analysis of variance test or Kruskal–Wallis test for comparison between multiple samples. Univariate and multivariate Cox regression analyses were performed to evaluate the prognostic values of costimulatory molecule genes. Moreover, we used Kaplan–Meier curves and log-rank tests to assess survival differences. Pearson’s chi-square test was used to assess differences in the distribution of clinical variables for PCa patients. All procedures involved in the present study were conducted using Perl and R scripts. A P-value < 0.05 was considered to indicate significance in all statistical tests.

Results

Identification of Costimulatory Molecule Genes With Prognostic Significance in PCa

First, we extracted the expression data of 60 costimulatory molecule genes in PCa from the TCGA database. These costimulatory molecules consisted of 13 B7-CD28 family genes and 47 TNF family genes, and the gene expression correlations among the 60 costimulatory molecule genes are shown in Supplementary Figure 1. Then, we performed a univariate Cox regression analysis to assess the prognostic relevance of the expression of these costimulatory molecule genes, and the screening criterion was P < 0.05. The results showed that a total of 14 costimulatory molecule genes were significantly associated with the prognosis of PCa. Among these fourteen significant genes, ten genes (TNFSF18, TNFRSF6B, TNFRSF18, TNFRSF25, CD80, CD86, CD70, RELT, LTA, and CD276) were recognized as risk factors with HR > 1, and four genes (TNFSF10, EDA2R, TNFSF13, and LTBR) were recognized as protective factors with HR < 1 (Table 1). Moreover, Kaplan-Meier curves were used to confirm the prognostic value of each of the fourteen genes (Figure 2). The results showed that high expression of the risk genes (TNFRSF18, TNFRSF6B, TNFSF18, TNFRSF25, CD80, CD86, CD70, RELT, and LTA) was associated with a poor prognosis in PCa. High expression of protective genes (TNFSF10 and EDA2R) was associated with a good prognosis in PCa. However, patients with high expression of CD276, TNFSF13, and LTBR had no significant difference in prognosis compared with patients with low expression. In addition, LASSO Cox regression analysis was executed to select genes with the highest prognostic value; five genes were selected, namely, RELT, EDA2R, TNFSF10, TNFSF18, and TNFRSF25 (Supplementary Figures 2A,B).

FIGURE 2

Identification of Two Clusters Using Consensus Clustering

To explore the prognostic value of the five costimulatory molecule genes, we performed a consensus clustering analysis to stratify PCa patients. Then, we executed a PCA to validate the reliability of the cluster numbers. The results demonstrated that k = 2 was the most stable value in the TCGA dataset (Figures 3A,B). The clustering heatmap and PCA clearly demarcated two clusters when k = 2 (Supplementary Figures 3A,B). Nevertheless, we also showed the results of consensus clustering when k = 3 (Supplementary Figure 3C), k = 4 (Supplementary Figure 3E), and k = 5 (Supplementary Figure 3G). The clusters overlapped with each other when there were three (Supplementary Figure S3D), four (Supplementary Figure 3F), or five (Supplementary Figure 3H) clusters. Therefore, patients were divided into two clusters (cluster 1 and cluster 2). The Kaplan–Meier curves revealed that the two clusters had different prognoses, and cluster 1 showed a worse prognosis than cluster 2 (Figure 3C). In addition, GSEA revealed that several oncogenic pathways, including the inflammatory response (normalized enrichment score: NES = 1.621), the interferon alpha response (NES = 1.566), the interferon gamma response (NES = 1.596), TNFA signaling via NFKB (NES = 1.377), IL6/JAK/STAT3 signaling (NES = 1.626), IL2/STAT5 signaling (NES = 1.449), epithelial mesenchymal transition (NES = 1.426), and angiogenesis (NES = 1.392), were significantly enriched in cluster 1 (Figure 3D).

FIGURE 3

Validation of the Five Survival-Related Costimulatory Molecule Genes

The expression levels of these five costimulatory molecule genes in PCa were compared between normal and tumor samples. RELT and TNFSF10 had high expression levels in tumor tissues compared with normal tissues, while EDA2R and TNFSF18 had low expression levels in tumor tissues compared with normal tissues in TCGA dataset. TNFRSF25 expression was not significantly difference between tumor tissues and normal tissues in the TCGA dataset (Figure 4A). In addition, we assessed the correlations between the expression levels of these five genes in different clinical subgroups. RELT and TNFSF18 were expressed at high levels, and EDA2R and TNFSF10 were expression at low levels in patients with advanced disease in terms of T stage. The expression level of TNFRSF25 showed no significant difference among patients with different T stages (Figure 4B). The expression levels of RELT, TNFSF18, and TNFRSF25 were high, and the expression levels of EDA2R were low in patients with lymphatic metastasis (Figure 4C). The expression levels of RELT, TNFSF18, and TNFRSF25 were high, and the expression levels of EDA2R and TNFSF10 were low in patients with a high Gleason score (Figure 4D).

FIGURE 4

Construction and Validation of the Prognostic Signature Based on Five Costimulatory Molecule Genes

We performed a multivariate Cox proportional hazards regression analysis based on these five survival-related costimulatory molecule genes. Subsequently, we constructed a prognostic model to stratify PCa patients based on the above multivariate Cox regression analysis results, integrating the expression profiles of five survival-related costimulatory molecule genes and their corresponding regression coefficients. A risk score was calculated as shown below:

Next, we used the optimal cutoff point for survival to stratify PCa patients into the high-risk and low-risk groups in all datasets. Kaplan–Meier analysis revealed that patients in the high-risk group had a significantly poorer prognosis than patients in the low-risk group (Figure 5A). We further validated these results in four GEO datasets. Similarly, significant differences were found in the GSE21034 (Figure 5D), GSE70768 (Figure 5J), and GSE70769 (Figure 5M) datasets. The same trend, although less significant difference, was observed in the GSE54460 dataset (Figure 5G). Moreover, the ROC curve showed good performance of these models for survival prediction in the TCGA dataset, and the area under the curve (AUC) was 0.725 at 1 year, 0.705 at 2 years, 0.743 at 3 years, and 0.745 at 5 years (Figure 5B). The efficiency of the risk model was further validated in the GSE21034 (Figure 5E), GSE54460 (Figure 5H), GSE70768 (Figure 5K), and GSE70769 (Figure 5N) datasets and showed good performance. PCA was performed to determine the distribution characteristics of the high-risk and low-risk groups. Different distributions for high-risk and low-risk patients were confirmed in the TCGA (Figure 5C), GSE21034 (Figure 5F), GSE54460 (Figure 5I), GSE70768 (Figure 5L), and GSE70769 (Figure 5O) datasets.

FIGURE 5

Associations Between the Prognostic Signature and Clinicopathological Factors of PCa

The heat map shows the expression of the five selected survival-related costimulatory molecule genes and clinicopathological factors in the high- and low-risk groups (Figure 6A). Moreover, the detailed distribution of the clinicopathological data across patient subgroups is shown in Table 2. The results showed that high-risk patients tended to have an advanced T stage, high prostate-specific antigen (PSA) levels, high Gleason scores and lymphatic metastasis. We used univariate and multivariate Cox regression analyses to determine whether the prognostic signature was an independent predictor factor for disease free survival in PCa patients. Univariate Cox regression analyses showed that patient age, pathological T stage, pathological N stage, Gleason score, PSA level, and risk score were significantly associated with prognosis (Figure 6B). Multivariate Cox regression analyses revealed that the Gleason score and risk score were significantly associated with prognosis (Figure 6C). These results demonstrated that the prognostic signature is an independent risk factor that can predict the prognosis of PCa patients. To determine the relationship between the risk signature and clinicopathological factors (age, pathological T stage, pathological N stage, Gleason score, and PSA level), patients were separated into different subgroups according to clinicopathological variables. The bar charts show that PCa patients with advanced age (Figure 6D), high pathological T stage (Figure 6E), lymph node metastasis (Figure 6F), high Gleason score (Figure 6G), high PSA levels (Figure 6H), and recurrence (Figure 6I) tended to have a high risk score. These results demonstrated that our prognostic signature is closely associated with the clinical factors of PCa.

FIGURE 6

Prognostic Signature-Related Biological Processes and Pathways

To explore signature-related biological pathways, genes strongly correlated with the risk score were filtered. A total of 525 positively related genes and 87 negatively related genes in the TCGA dataset (Pearson | R| > 0.4) were obtained (Figure 7A). The detailed results of the functional enrichment analysis are shown in Supplementary Table 2. The top ten results of the GO analyses are shown in Figures 7B–D. The most enriched terms for biological process, cellular components, and molecular function were “inflammatory response,” “membrane,” and “tumor necrosis factor-activated receptor activity,” respectively. According to the KEGG analysis, the most significantly enriched term was osteoclast differentiation (Figure 7E).

FIGURE 7

The Prognostic Signature Was Associated With the Immune Microenvironment

To gain a better understanding of the associations between our prognostic signature and the immune microenvironment, we displayed the expression profiles of 28 immune cell types in a heat map for high- and low-risk groups. The results showed significant differences in the immune cell infiltration status between high- and low-risk patients (Figure 8A). The box plots showed that high-risk patients had a high percentage of immune cells, including activated B cells, activated CD8 T cells, activated CD4 T cells, CD56bright natural killer cells, CD56dim natural killer cells, central memory CD4 T cells, central memory CD8 T cells, immature B cells, natural killer cells, type 1 T helper cells, activated dendritic cells, gamma delta T cells, macrophages, mast cells, myeloid-derived suppressor cells (MDSCs), plasmacytoid dendritic cells, regulatory T cells, T follicular helper cells, and type 2 T helper cells (Figure 8B). The immune and stromal scores for the TCGA cohort were calculated according to the ESTIMATE algorithm. We found that high-risk PCa patients had significantly higher immune scores, stromal scores, and ESTIMATE scores than low-risk patients (Figures 8C–E).

FIGURE 8

Differences in Genomic Alterations Between High-Risk and Low-Risk Patients

To investigate whether there were differences in genomic alterations between the 2 clusters, we compared the TMB and gene mutation patterns between high-risk and low-risk patients with available somatic mutation data in the TCGA entire set. Patients in the high-risk group had a higher TMB than those in the low-risk group (Figure 8F). The top 20 mutated genes for low-risk (Figure 8G) and high-risk (Figure 8H) patients were compared. We found that 49.55% of samples in the low-risk group and 63.67% of samples in the high-risk group had alterations in these genes. The forest plot illustrated that the TP53, STAB2, MUC17, PCDHB7, CUBN, CACNA1A, and MXRA5 genes were mutated at a significantly higher rate in high-risk patients (Figure 8I).

Discussion

Immunotherapy, as a rapidly growing field, aims to identify corresponding biomarkers prompting the immune system to recognize and kill cancer cells (Bilusic et al., 2017; Cha et al., 2020). Accumulating evidence indicated that the potential of immunotherapy in PCa treatment (Bilusic et al., 2017). However, effective biomarkers for predicting the prognosis of PCa patients remain scarce, and a single biomarker is insufficient to predict the clinical outcomes and response to immunotherapy (Smits et al., 2017). Previous studies have revealed that costimulatory molecules play an important role in the progression of various tumors (Loos et al., 2008; Zou and Chen, 2008; Geng et al., 2015). To improve the clinical therapy outcomes of PCa, we identified five costimulatory molecules and constructed a new prognostic signature for PCa patients based on these five genes. To our knowledge, our study provides the first prognostic signature of costimulatory molecules in patients with PCa. Moreover, the efficiency of our signature was well validated in four different GEO datasets. The multivariate Cox regression analysis indicated that our signature was an independent predictive factor for PCa patients and was significantly correlated with prognosis in different clinical subgroups. These results demonstrated that our prognostic signature is an effective and reliable prognostic tool. Additionally, we found that our prognostic signature was associated with the immune microenvironment. These findings should enhance the development of immunotherapeutic strategies for PCa patients.

The costimulatory molecules play an important role in the regulation of tumor immunity (Pitt et al., 2016). Topalian et al. (2020) reported that the costimulatory molecules expressed on tumor cells or lymphocytes play vital roles in regulating the antitumor immune response and that these molecules are closely related to the progression of tumors. At present, there are two main families of costimulatory molecules: the B7-CD28 family and the TNF family (Tang et al., 2018; Dostert et al., 2019). To determine the expression of the costimulatory molecule genes in PCa patients, we acquired data on 13 B7-CD28 family members and 47 TNF family members for our study (Janakiram et al., 2015; Dostert et al., 2019). After using univariate and multivariate Cox proportional hazards regression analyses, we identified five costimulatory molecular genes (RELT, TNFRSF25, EDA2R, TNFSF18, and TNFSF10) with prognostic value. However, the functions and roles of these costimulatory molecules in PCa remain unclear. RELT, a new member of the TNFRSF family, accelerates tumor progression and regulates the infiltration of numerous immune cell types. RELT is significantly upregulated as glioma grade increases and is associated with a poor prognosis (Jin et al., 2020). TNFRSF25, also known as death receptor 3 (DR3), belongs to the TNFRSF family and is highly expressed in T cells. Recent studies have shown that DR3 is a potential immunotherapy molecular target for cancer treatment and plays essential roles in protective inflammation, autoimmune diseases, and tumor immunotherapy (Mavers et al., 2019; So and Ishii, 2019). EDA2R, also known as TNFRSF27 or XEDAR, mainly coordinates various cellular and organismal biological processes and exerts its roles by activating gene transcription. Studies have indicated that EDA2R is a direct p53 target that can be activated by p53 (Brosh et al., 2010). They also confirmed that treatment of cancer cells with the ligand EDA-A2, which can specifically activate EDA2R, leads to p53-dependent cell death. In addition, Tanikawa et al. (2010) found that EDA2R is frequently downregulated in colorectal cancer patients due to epigenetic alterations. Interestingly, this is consistent with our findings that EDA2R is also markedly downregulated in PCa patients in the TCGA dataset. TNFSF18 is one of the TNFRSF members expressed by myeloid cells and provides costimulatory signals to boost T cell activity. The blockade of TNFSF18-GITR signaling to target mesothelioma stem cells might be translated into a therapeutic strategy for mesothelioma treatment (Wu et al., 2019). Recent evidence has further shown that TNFSF10 is an apoptosis-inducing ligand and can promote the apoptosis of cancer cells (He et al., 2012). Sun et al. (2018) revealed that the abnormal expression of TNFSF10 influences the apoptosis of PCa cells. Furthermore, relevant studies have shown that TNFSF10 is a promising anticancer agent for cancer that exhibits good anticancer activity by inducing apoptosis (Janakiram et al., 2015; Dostert et al., 2019).

Based on these five costimulatory molecular genes, we developed a prognostic signature for PCa patients. The performance of our prognostic signature was validated in four different independent GEO datasets. The results showed good performance in three datasets (GSE21034, GSE70768, and GSE70769). The P-value in the GSE54460 dataset did not reach statistical significance, which might be attributed to the small sample size. In addition, we further investigated the association between the prognostic signature and clinical factors. We found the high-risk patients were older and had more advanced T stage, a higher rate of lymphatic metastasis, higher Gleason scores, higher PSA levels and a higher mortality rate than low-risk patients. These results revealed that our prognostic signature is closely correlated with clinical factors. Therefore, this prognostic signature could be applied as a supplement for guiding treatment and improving the clinical prognosis of PCa patients. In general, our prognostic signature was well validated in different GEO datasets, which encouraged our investigation of the underlying molecular mechanisms. Through GO and KEGG enrichment analysis of these costimulatory molecules, we found that the potential molecular mechanisms of these costimulatory molecules in PCa were closely related to the immune pathway, as indicated by the enrichment of terms such as “regulation of immune response,” “innate immune response,” “immune response,” “adaptive immune response,” “immunological synapse,” and “positive regulation of T cell proliferation.” These results indicate that immune heterogeneity may be the cause of the difference in prognosis between patient groups.

Subsequently, we explored the immune cell infiltration and tumor mutation profiles of the high-risk and low-risk groups to further reveal the difference in the immune microenvironment between these two groups. The results showed that the high-risk group exhibited significantly greater infiltration of immune cells, including activated B cells, activated CD8 T cells, activated CD4 T cells, CD56bright natural killer cells, CD56dim natural killer cells, central memory CD4 T cells, central memory CD8 T cells, immature B cells, natural killer cells and type 1 T helper cells, than the low-risk group. In addition, the infiltration of various immunosuppressive cells, including activated dendritic cells, gamma delta T cells, macrophages, mast cells, MDSCs, plasmacytoid dendritic cells, regulatory T cells, T follicular helper cells, and type 2 T helper cells, was greater in high-risk patients than in low-risk patients, suggesting the presence of an immunosuppressive microenvironment in high-risk patients. Complex interactions between immunosuppressive cells cooperate to suppress antitumor immune responses and promote disease progression. Targeting regulatory T cell function or the secretion of immunological processes can lead to tumor immune evasion (Vinay et al., 2015). Su et al. (2019) found that CCL2 induced the recruitment of M2-like tumor-associated macrophages and regulatory T cells, thus promoting metastasis with immune suppression and neoangiogenesis in PCa. Sullivan et al. (2021) studied Hi-Myc mice crossed to mast cell knockout mice and demonstrated that higher levels of mast cell infiltration led to the promotion of cancer cell invasion. MDSCs are known to play critical roles in tumor immune evasion. Lu et al. (2017) found that MDSCs promoted the initiation and progression of PCa in a mouse model. In addition, the level of MDSC infiltration correlated with the PSA levels and metastasis in PCa patients (Lu et al., 2017). Understanding the immune microenvironment of each PCa patient can help us identify patients who are more likely to benefit from immunotherapy.

Similarly, we found that the tumor mutation burden (TMB) in the high-risk group was also higher than that in the low-risk group, especially in terms of TP53 mutation frequency, which was obviously increased in the high-risk group. Some scholars have found that sequencing technology indicates that the TP53 mutation frequency is much higher than that reported in the TCGA database (Mateo et al., 2020). Martin et al. (2011) revealed that prostate epithelial Pten/TP53 loss leads to epithelial to mesenchymal transition. TP53 mutation is linked to a higher level of tumor-infiltrating T cells, which influences the immunotherapy response in prostate cancer (Kaur et al., 2019). Mu et al. (2017) reported that PTEN loss with TP53 mutation causes resistance to antiandrogen therapy of PCa. In addition, Jiang et al. (2018) demonstrated that TP53 mutation could result in an immunosuppressed state in gastric cancer. In the present study, the immunosuppressive microenvironment in high-risk patients might be partly due to the high frequency of TP53 mutations in high-risk patients. The above conclusions revealed that these costimulatory molecules have important prognostic values for patients with PCa. These findings give us sufficient confidence that our signature can be applied as a novel strategy for guiding treatment and improving clinical therapy outcomes.

In our study, we constructed a costimulatory molecule-related prognostic signature for PCa, which could be used to stratify patients to further guide treatment and improve clinical outcomes. This study was the first comprehensive study of the expression profiles and clinical significance of costimulatory molecule genes in PCa patients. Although our study provides important insights to better evaluate costimulatory molecules and the prognosis of PCa patients, it inevitably has some limitations that need to be noted. First, regardless of the fact that we used four different independent datasets for validation, the present study was a retrospective study. All data were obtained from the public databases. Moreover, our research was entirely conducted through a series of bioinformatics methods, so experimental and prospective studies are needed to further confirm the good predictive ability of our prognostic signature.

Conclusion

In this study, we performed the first comprehensive analysis of costimulatory molecules in patients with PCa through a series of bioinformatics analysis methods. We identified several key prognostic costimulatory molecule genes, built a reliable and valid prognostic signature, and explored the potential molecular mechanisms of this signature. Our prognostic signature could stratify PCa patients into two subgroups with different prognoses and showed high associations with the clinical features. Moreover, patients identified as high risk based on our prognostic signature exhibited a high mutation frequency, a high level of immune cell infiltration, and an immunosuppressive microenvironment. Thus, our signature could provide clinicians with prognosis predictions and treatment guidance for PCa patients.

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.

References

• 1

  AbidaW.ChengM. L.ArmeniaJ.MiddhaS.AutioK. A.VargasH. A.et al (2019). Analysis of the prevalence of microsatellite instability in prostate cancer and response to immune checkpoint blockade.JAMA Oncol.5471–478. 10.1001/jamaoncol.2018.5801

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  BilusicM.MadanR. A.GulleyJ. L. (2017). Immunotherapy of prostate cancer: Facts and Hopes.Clin. Cancer Res.236764–6770. 10.1158/1078-0432.CCR-17-0019

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  BroshR.SarigR.NatanE. B.MolchadskyA.MadarS.BornsteinC.et al (2010). P53-dependent transcriptional regulation of EDA2R and its involvement in chemotherapy-induced hair loss.FEBS. Lett.5842473–2477. 10.1016/j.febslet.2010.04.058

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  CaseyS. C.AmedeiA.AquilanoK.AzmiA. S.BenenciaF.BhaktaD.et al (2015). Cancer prevention and therapy through the modulation of the tumor microenvironment.Semin. Cancer Biol.35S199–S223. 10.1016/j.semcancer.2015.02.007

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  ChaH. R.LeeJ. H.PonnazhaganS. (2020). Revisiting immunotherapy: a focus on prostate cancer.Cancer. Res.801615–1623. 10.1158/0008-5472.CAN-19-2948

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  ChalmersZ. R.ConnellyC. F.FabrizioD.GayL.AliS. M.EnnisR.et al (2017). Analysis of 100,000 human cancer genomes reveals the landscape of tumor mutational burden.Genome. Med.9:34. 10.1186/s13073-017-0424-2

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  CharoentongP.FinotelloF.AngelovaM.MayerC.EfremovaM.RiederD.et al (2017). Pan-cancer immunogenomic analyses reveal Genotype-Immunophenotype relationships and predictors of response to checkpoint blockade.Cell. Rep.18248–262. 10.1016/j.celrep.2016.12.019

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  CroftM.BenedictC. A.WareC. F. (2013). Clinical targeting of the TNF and TNFR superfamilies.Nat. Rev. Drug. Discov.12147–168. 10.1038/nrd3930

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  DostertC.GrusdatM.LetellierE.BrennerD. (2019). The TNF family of ligands and receptors: Communication modules in the immune system and beyond.Physiol. Rev.99115–160. 10.1152/physrev.00045.2017

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  GengY.WangH.LuC.LiQ.XuB.JiangJ.et al (2015). Expression of costimulatory molecules B7-H1, B7-H4 and Foxp3+ Tregs in gastric cancer and its clinical significance.Int. J. Clin. Oncol.20273–281. 10.1007/s10147-014-0701-7

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  GrassoC. S.WuY. M.RobinsonD. R.CaoX.DhanasekaranS. M.KhanA. P.et al (2012). The mutational landscape of lethal castration-resistant prostate cancer.Nature487239–243. 10.1038/nature11125

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  GuL.FrommelS. C.OakesC. C.SimonR.GruppK.GerigC. Y.et al (2015). BAZ2A (TIP5) is involved in epigenetic alterations in prostate cancer and its overexpression predicts disease recurrence.Nat. Genet.4722–30. 10.1038/ng.3165

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  GuL.WangL.ChenH.HongJ.ShenZ.DhallA.et al (2020). CG14906 (mettl4) mediates m6A methylation of U2 snRNA in Drosophila.Cell. Discov.6:44. 10.1038/s41421-020-0178-7

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  HeW.WangQ.XuJ.XuX.PadillaM. T.RenG.et al (2012). Attenuation of TNFSF10/TRAIL-induced apoptosis by an autophagic survival pathway involving TRAF2- and RIPK1/RIP1-mediated MAPK8/JNK activation.Autophagy81811–1821. 10.4161/auto.22145

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  JanakiramM.ChinaiJ. M.ZhaoA.SparanoJ. A.ZangX. (2015). HHLA2 and TMIGD2: New immunotherapeutic targets of the B7 and CD28 families.Oncoimmunology4:e1026534. 10.1080/2162402X.2015.1026534

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  JansenC. S.ProkhnevskaN.KissickH. T. (2019). The requirement for immune infiltration and organization in the tumor microenvironment for successful immunotherapy in prostate cancer.Urol. Oncol.37543–555. 10.1016/j.urolonc.2018.10.011

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  JiangZ.LiuZ.LiM.ChenC.WangX. (2018). Immunogenomics analysis reveals that TP53 mutations inhibit tumor immunity in gastric cancer.Transl. Oncol.111171–1187. 10.1016/j.tranon.2018.07.012

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  JinX.XieH.LiuX.ShenQ.WangZ.HaoH.et al (2020). RELL1, a novel oncogene, accelerates tumor progression and regulates immune infiltrates in glioma.Int. Immunopharmacol.87:106707. 10.1016/j.intimp.2020.106707

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  KaurH. B.LuJ.GuedesL. B.MaldonadoL.ReitzL.BarberJ. R.et al (2019). TP53 missense mutation is associated with increased tumor-infiltrating T cells in primary prostate cancer.Hum. Pathol.8795–102. 10.1016/j.humpath.2019.02.006

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  LoosM.GieseN. A.KleeffJ.GieseT.GaidaM. M.BergmannF.et al (2008). Clinical significance and regulation of the costimulatory molecule B7-H1 in pancreatic cancer.Cancer Lett.26898–109. 10.1016/j.canlet.2008.03.056

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  LuX.HornerJ. W.PaulE.ShangX.TroncosoP.DengP.et al (2017). Effective combinatorial immunotherapy for castration-resistant prostate cancer.Nature543728–732. 10.1038/nature21676

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  MartinP.LiuY. N.PierceR.Abou-KheirW.CaseyO.SengV.et al (2011). Prostate epithelial Pten/TP53 loss leads to transformation of multipotential progenitors and epithelial to mesenchymal transition.Am. J. Pathol.179422–435. 10.1016/j.ajpath.2011.03.035

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  MateoJ.SeedG.BertanC.RescignoP.DollingD.FigueiredoI.et al (2020). Genomics of lethal prostate cancer at diagnosis and castration resistance.J. Clin. Invest.1301743–1751. 10.1172/JCI132031

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  MaversM.SimonettaF.NishikiiH.RibadoJ. V.Maas-BauerK.AlvarezM.et al (2019). Activation of the DR3-TL1A axis in donor mice leads to regulatory t cell expansion and activation with reduction in Graft-Versus-Host disease.Front. Immunol.10:1624. 10.3389/fimmu.2019.01624

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  MuP.ZhangZ.BenelliM.KarthausW. R.HooverE.ChenC. C.et al (2017). SOX2 promotes lineage plasticity and antiandrogen resistance in TP53- and RB1-deficient prostate cancer.Science35584–88. 10.1126/science.aah4307

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  Nguyen-NielsenM.BorreM. (2016). Diagnostic and therapeutic strategies for prostate cancer.Semin. Nucl. Med.46484–490. 10.1053/j.semnuclmed.2016.07.002

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  OmlinA.PezaroC.MukherjiD.CassidyA. M.SandhuS.BianchiniD.et al (2013). Improved survival in a cohort of trial participants with metastatic castration-resistant prostate cancer demonstrates the need for updated prognostic nomograms.Eur. Urol.64300–306. 10.1016/j.eururo.2012.12.029

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  PittJ. M.AndréF.AmigorenaS.SoriaJ. C.EggermontA.KroemerG.et al (2016). Dendritic cell-derived exosomes for cancer therapy.J. Clin. Invest.1261224–1232. 10.1172/JCI81137

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  SchildbergF. A.KleinS. R.FreemanG. J.SharpeA. H. (2016). Coinhibitory pathways in the B7-CD28 Ligand-Receptor family.Immunity44955–972. 10.1016/j.immuni.2016.05.002

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  SharmaP.Hu-LieskovanS.WargoJ. A.RibasA. (2017). Primary, adaptive, and acquired resistance to cancer immunotherapy.Cell168707–723. 10.1016/j.cell.2017.01.017

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  SiegelR. L.MillerK. D.JemalA. (2020). Cancer statistics, 2020.CA. Cancer. J. Clin.707–30. 10.3322/caac.21590

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 32

  SmitsM.MehraN.SedelaarM.GerritsenW.SchalkenJ. A. (2017). Molecular biomarkers to guide precision medicine in localized prostate cancer.Expert. Rev. Mol. Diagn.17791–804. 10.1080/14737159.2017.1345627

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  SoT.IshiiN. (2019). The TNF-TNFR family of co-signal molecules.Adv. Exp. Med. Biol.118953–84. 10.1007/978-981-32-9717-3_3

  • CrossRef
  • Google Scholar
  • View reference in article

• 34

  SuW.HanH. H.WangY.ZhangB.ZhouB.ChengY.et al (2019). The polycomb repressor complex 1 drives Double-Negative prostate cancer metastasis by coordinating stemness and immune suppression.Cancer Cell.36139–155. 10.1016/j.ccell.2019.06.009

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 35

  SullivanH. H.MaynardJ. P.HeaphyC. M.LuJ.MarzoA. M. D.LotanT. L.et al (2021). Differential mast cell phenotypes in benign versus cancer tissues and prostate cancer oncologic outcomes.J. Pathol.253415–426. 10.1002/path.5606

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 36

  SunM.GengD.LiS.ChenZ.ZhaoW. (2018). LncRNA PART1 modulates toll-like receptor pathways to influence cell proliferation and apoptosis in prostate cancer cells.Biol. Chem.399387–395. 10.1515/hsz-2017-0255

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 37

  TangJ.JiangW.LiuD.LuoJ.WuX.PanZ.et al (2018). The comprehensive molecular landscape of the immunologic co-stimulator B7 and TNFR ligand receptor families in colorectal cancer: Immunotherapeutic implications with microsatellite instability.Oncoimmunology7:e1488566. 10.1080/2162402X.2018.1488566

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 38

  TanikawaC.RiC.KumarV.NakamuraY.MatsudaK. (2010). Crosstalk of EDA-A2/XEDAR in the p53 signaling pathway.Mol. Cancer Res.8855–863. 10.1158/1541-7786.MCR-09-0484

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 39

  TibshiraniR. (1997). The lasso method for variable selection in the Cox model.Stat. Med.16385–395. 10.1002/(sici)1097-0258(19970228)16:4<385::aid-sim380<3.0.co;2-3

  • CrossRef
  • Google Scholar
  • View reference in article

• 40

  TopalianS. L.TaubeJ. M.PardollD. M. (2020). Neoadjuvant checkpoint blockade for cancer immunotherapy.Science367:eaax0182. 10.1126/science.aax0182

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 41

  VinayD. S.RyanE. P.PawelecG.TalibW. H.StaggJ.ElkordE.et al (2015). Immune evasion in cancer: Mechanistic basis and therapeutic strategies.Semin. Cancer Biol.35S185–S198. 10.1016/j.semcancer.2015.03.004

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 42

  WangG.ZhaoD.SpringD. J.DePinhoR. A. (2018). Genetics and biology of prostate cancer.Gene. Dev.321105–1140. 10.1101/gad.315739.118

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 43

  WangQ.GuL.AdeyA.RadlwimmerB.WangW.HovestadtV.et al (2013). Tagmentation-based whole-genome bisulfite sequencing.Nat. Protoc.82022–2032. 10.1038/nprot.2013.118

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 44

  WeiS. C.DuffyC. R.AllisonJ. P. (2018). Fundamental mechanisms of immune checkpoint blockade therapy.Cancer Discov.81069–1086. 10.1158/2159-8290.CD-18-0367

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 45

  WuL.Dell’AnnoI.LapidotM.SekidoY.ChanM. L.KohnoM.et al (2019). Progress of malignant mesothelioma research in basic science: A review of the 14th international conference of the international mesothelioma interest group (iMig2018).Lung. Cancer127138–145. 10.1016/j.lungcan.2018.11.034

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 46

  WuZ.ChenH.LuoW.ZhangH.LiG.ZengF.et al (2020). The landscape of immune cells infiltrating in prostate cancer.Front. Oncol.10:517637. 10.3389/fonc.2020.517637

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 47

  YoshiharaK.ShahmoradgoliM.MartínezE.VegesnaR.KimH.Torres-GarciaW.et al (2013). Inferring tumour purity and stromal and immune cell admixture from expression data.Nat. Commun.4:2612. 10.1038/ncomms3612

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 48

  ZhangC.ZhangZ.LiF.ShenZ.QiaoY.LiL.et al (2018). Large-scale analysis reveals the specific clinical and immune features of B7-H3 in glioma.Oncoimmunology7:e1461304. 10.1080/2162402X.2018.1461304

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 49

  ZouW.ChenL. (2008). Inhibitory B7-family molecules in the tumour microenvironment.Nat. Rev. Immunol.8467–477. 10.1038/nri2326

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article