Immune-related biomarker risk score predicts prognosis in prostate cancer

In this study, we constructed a model using a Cox proportional hazards model based on the expression of eight immune-related genes that were associated with prognosis in prostate cancer: EDNRB, ANGPTL2, TNFSF15, TNFRSF10D, EDN2, BMP2, NLRP14, and PLK1. We then identified associations between risk scores calculated with the model, tumor microenvironment characteristics, and immune cell infiltration. Prostate cancer patients in the high score group had poorer prognoses, and validation with the external GSE54460 dataset confirmed that the scoring model predicted biochemical recurrence with AUC values of 0.749 at 1 year, 0.804 at 3 years, and 0.774 at 5 years. Proportions of infiltrated M2 macrophages and regulatory T cells were increased in the high risk group, while CD8+ T cells were increased in the low risk group. Network analysis revealed that PLK1 may be a key regulator of the immune-suppressive microenvironment in prostate cancer. Double immunofluorescence labeling of a prostate cancer tissue microarray indicated that PLK1 expression correlated positively with numbers of infiltrating macrophages. These results indicate that an immune- related, gene-based risk score effectively reflects immune microenvironment characteristics and predicts prognosis in prostate cancer.


INTRODUCTION
Prostate cancer poses a serious threat to the health of men all over the world [1,2]. Although radical surgical excision and radiation therapy are effective in treating prostate cancer, only lung cancer has a higher mortality rate. In China, where prostate cancer screening is relatively uncommon, most patients have locally advanced or metastatic prostate cancer upon diagnosis. In these cases, endocrine therapy is typically used as the primary treatment. After about 18 months of endocrine therapy, most prostate cancers become resistant to hormone treatments, and no other effective clinical treatments are currently available. Novel immunotherapies are an important development in tumor treatment, and immune checkpoint inhibitors may be particularly effective. However, recent phase II clinical trials indicate that immune check point inhibitors are only effective against specific types of prostate cancer, and the disease control rate does not exceed 20% [3,4].
The development of resistance to immunotherapies is the main reason for their poor efficacy in treating prostate cancer. PD-L1 + or VISTA + M2 macrophages are the major drivers of prostate immunotherapy resistance [5]. M1 macrophages activated by the classical pathway can induce tissue inflammation, and T cells activated by inflammation can effectively kill or suppress prostate cancer cells. In contrast, M2 macrophages with anti-inflammatory characteristics can AGING promote tissue repair and immunosuppression. M2 macrophages exert protective effects in tumors by inhibiting tumor antigen presentation. This prevents T cells from recognizing tumor antigens and immune checkpoint inhibitors, which regulate T cell function, from effectively treating prostate cancer [6][7][8]. Additional causes of insensitivity to immunotherapy in prostate cancer remain poorly understood. Further explorations of the immune characteristics of prostate cancer are therefore needed.
Interactions between tumor cells and immune cells are complex, and tumor status can alter the tumor immune microenvironment. In this study, we examined relationships between gene expression and immune cell infiltration in prostate cancer. Immune genes associated with biochemical recurrence of prostate cancer were screened to construct a prognostic model. The effectiveness of the model was validated using an external dataset. Relationships between immune characteristic genes and immune cell infiltration in prostate cancer were analyzed. Finally, key genes that might induce crucial changes in the prostate cancer immune microenvironment were identified ( Figure 1). The relationship between expression of these potential key genes and M2 macrophage infiltration was verified in tissue samples.

RESULTS
Weighted correlation network analysis (WGCNA) identified prostate cancer-related genes associated with biochemical recurrence (BCR) The TCGA database includes 5132 genes that are differentially expressed in prostate cancer compared to normal prostate tissue. Among them, 602 were immuneregulated genes (IRGs). Weighted correlation network analysis (WGCNA) was performed using these 602 candidate genes based on TCGA PRAD cancer sample expression profiles ( Figure 2). Eight samples with abnormal clustering were removed during the screening process. The soft-thresholding power in the WGCNA (β=5) was determined based on scale-free R2 > 0.85 ( Figure 3A). Eight modules were identified by average AGING linkage hierarchical clustering based the soft-thresholding power ( Figure 3B). According to the hierarchical clustering of modules, there were significant correlations between most modules ( Figure 3C). A network heatmap of all genes is shown in Figure 3D.

Construction of an immunogenetic risk score associated with BCR in prostate cancer
We analyzed associations between all modules and clinical characteristics of prostate cancer patients. The brown, turquoise, pink, and yellow modules were highly associated with time to BCR, and the pink and red modules were also associated with BCR ( Figure 4). Further screening of the genes included in these modules revealed a correlation coefficients of greater than 0.1 between their expression and BCR in prostate cancer and of greater than 0.5 for gene expression correlations within the modules. In total, 221 immunerelated genes were associated with BCR in prostate cancer (Supplementary Table 3). We then performed Cox univariate regression analysis to identify correlations between the above genes and BCR based on clinical information for prostate cancer patients in TCGA. A log-rank test revealed that 53 genes were associated with BCR in prostate cancer. Among these 53 genes, eight survival-associated IRGs (EDNRB, ANGPTL2, TNFSF15, TNFRSF10D, EDN2, BMP2, NLRP14, PLK1) were also identified using lasso regression analysis ( Figure 5A, 5B). KM curves for PLK1, NLRP14, TNFRSF10D, and FGFR2 are shown in Figure 5C-5F; lasso regression coefficients for all eight IRGs are show in Table 1.

Immunogenetic risk score is associated with BCR in prostate cancer
Risk scores for BCR in prostate cancer based on corresponding lasso coefficients were calculated for each sample in the training cohort (GSE54460). Risk score was significantly associated with BCR in prostate cancer patients; BCR occurred sooner in the high risk group ( Figure 6A). An ROC curve was used to assess the effect of the Risk score. The AUC (Area Under Curve) was 0.749 at one year, 0.804 at 3 years, and 0.774 at 5 years in the testing cohort ( Figure 6B). The ROC curve in the TCGA cohort was consistent with the curve in the training cohort ( Figure 6C). The AUC was  Figure 6D). Risk score was therefore a useful predictor of BCR in prostate cancer. We also performed a multivariate correlation analysis to determine the impact of other clinical factors on the prognostic power of the risk score. Immune risk score was still associated with BCR in prostate cancer after multivariate adjustment (p=0.042) ( Table 2).

Differences in immune characteristics between high and low risk groups
Next, we examined differences in the tumor microenvironment between low and high risk patients based on prostate cancer mRNA expression data from the TCGA database. Genes that were differentially expressed between the high and low risk score groups are shown in the volcano plot in Supplementary Figure  1. No significant differences in immune and stroma genes were observed between the high and low risk score groups (Supplementary Figure 2). However, differences were observed in infiltration levels for 22 types immune cells in CIBERSORT data. Memory B cell, regulatory T cell, M2 macrophage, and dendritic cell infiltration were higher in the high risk group, while plasma cell, CD8 + T cell, monocyte, and activated mast cell infiltration were higher in the low score group ( Figure 7A). Overall immune cell infiltration data for the high and low risk score groups is shown in Figure  7B. Expression of immunoregulatory factors in the AGING tumor immune microenvironment also differed between the high and low risk groups. ( Figure 8A); the 28 factors that differed significantly are shown in Supplementary Figure 3.

Identification of pathways and gene ontology (GO) terms associated with high and low risk score groups
The ClusterProfiler R package was used to perform GSEA enrichment analysis on DEGs between the high and low score groups using gene sets from MsigDB as background genes. Log2(fold change) values from the differential expression analysis were used as the sorting criterion, and P value < 0.05 was used as the screening criterion (Supplementary Table 4). The results indicated that there were significant differences in "immune effector process," "immune response," "immune system process," "innate immune response," "positive regulation of immune system process," "regulation of immune system process," and "regulation of immune system process" between low and high risk groups ( Figure 8B). GO enrichment analysis was also performed using p value < 0.05 and overlap > 0.75 as screening criteria (Supplementary Figure 4).

PLK1 may be a key regulator of the immune microenvironment in prostate cancer
The LnCeVar database contains transcription factor regulation data validated in the literature. As shown in Figure 9, PLK1 is part of a rich regulation system. In a previous study, we found that infiltration of M2 macrophages promotes prostate cancer progression. Here, risk score was also associated with M2 macrophages. Among the genes included in the risk score, PLK1 had the largest coefficient, indicating that it may be the most critical factor in the model. We therefore performed double immunofluorescence labeling of CD163, a marker of M2 macrophages, and PLK1 in a prostate cancer tissue microarray ( Figure  10A). PLK1 was positively correlated with CD163 in prostate cancer (Pca) samples (r=0.69, p<0.01), but not in benign prostatic hypertrophy (BPH) samples (r=0.12, p=0.63) ( Figure 10B). Furthermore, PLK1 expression was significantly increased in prostate cancer compared to prostatic hyperplasia. (p=0.022) ( Figure 10C).

DISCUSSION
Current treatments for prostate cancer include surgery, radiation therapy, and endocrine therapy. Recent studies  and non-functional T cells (CD4 + and CD8 + ) can induce immunosuppression in the tumor microenvironment [9]. Such immunosuppression is the main obstacle to the efficacy of immunotherapy-induced anti-tumor immune responses. Multi-drug combination therapy, which can     AGING increase the impact of immune checkpoint inhibitors by regulating the tumor microenvironment, is an important method for improving the efficacy of cancer immunotherapy [10]. In this study, we examined the immune characteristics of prostate cancer to identify potential key genes that regulate the tumor immune microenvironment.
We found that 53 immune-related genes were associated with biochemical recurrence in prostate cancer. Among these, eight were identified as key genes by lasso regression. Validation with an external dataset indicated that the prediction model was highly accurate (5 year AUC=0.774). Moreover, biochemical recurrence occurred sooner and survival rates were lower in the high risk score group. Analysis of the LnCeVar database indicated that PLK1 had a particularly rich regulatory system in prostate cancer among the eight key genes.
Decreased expression of EDNRB, TNFSF15, TNFRSF10D, EDN2, and NLRP14, as well as increased expression of ANGPTL2, BMP2, and PLK1, were associated with higher risk scores; EDNRB promoter methylation status was also associated with risk score. A meta-analysis of 11 similar studies indicated that the frequency of EDNRB methylation was substantially higher in prostate cancer compared with normal prostate tissues (OR = 5.42, 95 % CI = 1.98-14.88, P = 0.001), suggesting that EDNRB promoter methylation might increase the risk of prostate cancer [11]. Here, TNFSF15 expression was inversely associated with prostate cancer risk. This is consistent with the role of TNFSF15 as a downstream effector of AMPK that AGING inhibits prostate cancer growth [12]. TNFRSF10D expression is also associated with prostate cancer and with the direct p53 effectors and ERK signaling pathways; here, TNFRSF10D was inversely related to prostate cancer biochemical recurrence risk score. Endothelins are involved in the regulation of various physiological processes, including plumage development in chickens, pigmentation, neural crest cell proliferation, differentiation, migration, cardiovascular development and functions, and pulmonary hypertension [13]. The endothelin EDN2 was inversely associated with prostate cancer immunological risk score in our model. As a member of a family of molecules that belong to a signal-induced multiprotein complex termed the inflammasome that activates proinflammatory caspase-1 and caspase-5, NLRP14 may play a regulatory role in the innate immune system. NLRP14 is considered an oncogene, and increased expression of NLRP14 is associated with increases in prostate cancer mortality [14]. This contradicts the inverse association found here between NLRP14 and risk score, and further study is warranted. ANGPTL2 is a secreted glycoprotein with homology to angiopoietins that may exert autocrine or paracrine effects on endothelial cells. ANGPTL2 also promotes M2 polarization of macrophages in non-small cell lung cancer [15]. In addition, ANGPTL2 may promote acquisition of androgen independence and tumor progression in prostate cancer by exerting autocrine and/or paracrine effects via the integrin α5β1 receptor [16]. Here, we found a positive association between ANGPTL2 expression and immune-related risk score in prostate cancer. Finally, several studies suggest that BMP2 promotes progression and induces biochemical recurrence in prostate cancer [17][18][19], which is consistent with the positive association found here between BMP2 and risk score.
Our present results demonstrate that memory B cell, regulatory T cell, M2 macrophage, and dendritic cell infiltration were significantly increased in the high risk score group, while plasma cell, CD8+T cell, monocyte, and activated mast cell infiltration were higher in the low risk score group. M2 macrophages promote prostate cancer progression and help establish an immunosuppressive state in tumors [20][21]; this might help explain the increased M2 macrophage infiltration observed in the higher risk score group. In contrast, local increases in the density of infiltrating CD8+T cells in tumors is a marker of good prognosis; this might account for the increased CD8+T cell infiltration observed here in the low risk score group with better prognoses. Infiltration of these cell types reflects the immune microenvironment and can also predict prognosis in prostate cancer [22].
The role of PLK1 in prostate cancer is not clear at present [23]. However, numerous studies indicate that PLK1 can act as an oncogene, and PLK-1 inhibitors can effectively inhibit prostate cancer progression [24][25][26]. PLK1 can also inactivate other tumor suppressors [27,28]. Another study suggests that PLK1 is a carcinogenic factor [29]. Here, we found that PLK1 expression was higher in the high risk group that experienced earlier biochemical recurrence. Moreover, downstream pathways regulated by PLK1 comprised the single largest pathway group in our biochemical recurrence prediction model for prostate cancer. PLK1 might therefore be one of the most important immune genes that contribute to biochemical recurrence in prostate cancer.

Data acquisition
The training dataset was obtained from TCGA and included 498 cancer samples and 52 normal control samples; clinical information for those samples is shown in Supplementary Table 1. The GSE54460 testing dataset containing 106 samples was downloaded from the gene expression omnibus database (GEO: https://www.ncbi. nlm.nih.gov/geo/); clinical information for those samples is shown in Supplementary Table 2. Immune genes that were included in our analyses were obtained from InnateDB (https://www.innatedb.ca/) and ImmPort (https://www.immport.org/home).

Identification of differentially expressed genes
Genes that were differentially expressed between normal and prostate cancer samples were identified based on the screening criteria of p value > 0.05 and log(fold-change) > 1.5.

Weighted gene co-expression network analysis
TCGA expression data for 602 immune-related genes was used for WGCNA (weighted gene co-expression network analysis) to identify associations between gene expression modules and clinical characteristics [30]. During sample clustering, eight samples with abnormal clustering were identified and removed from the analysis (Supplementary Figure 1). The coexpression network was then constructed and divided into modules. A scale-free network coefficient of greater than 0.85 was used to ensure that the coexpression network conformed to the scale-free network standards. Gene significance (GS) indicates the strength of linear correlations between the expression of different gene modules and clinical features. Modules with P≤0.01 and higher GS values AGING were identified as survival-related modules and included in subsequent analysis.

Univariate Cox regression and lasso regression
The Survival package for R was used to perform univariate Cox regression and the KM test. The Glmnet package was used for lasso regression. Prognosis related genes were screen and correlation regression coefficients were obtained. The resulting risk score was validated in the GEO dataset (GSE54460). "TimeROC" was used to draw receiver operating characteristic curves (ROC), and the area under the curve (AUC) was calculated. Samples were divided between the high and low risk groups using the median risk score (RS) as a cutoff value. KM curves were used to evaluate the survival of prostate cancer patients.

Comparison of immune cell infiltration between high and low risk groups
The "ESTIMATE" package was used to calculate the microenvironment score, and the "CIBERSORT" package was used to assess the proportions of 22 leukocyte subtypes based on differences in mRNA expression between the high and low risk groups [31]. The Wilcox test was used to evaluate differences between the high low risk score groups.

Statistical analysis
All statistical analyses were performed using R software (version 3.6.3, http://www.R-project.org). A two-sided P < 0.05 indicated a statistically significant difference.

Supplementary Figures
Supplementary Figure 1. Volcano map showing genes that were differential expressed between the high and low risk score groups. Red dots represent genes that are up-regulated in the high risk group, and blue dots represent genes that are downregulated in the high risk group.

Supplementary Tables
Please browse Full Text version to see the data of Supplementary Tables 3 and 4.