Identification of hypoxia-related differentially expressed genes (DE-HRGs) in PCa
To assess the hypoxia status, we first calculated the hypoxia score of each patient with PCa using the maximally selected rank statistics method. The optimal cut-off value of hypoxia score was 0·48, which was used to assign patients with PCa from the TCGA database to high- and low-hypoxia score groups(Figure 1A, Table S1). As shown in Figure 1B, patients with PCa and a high hypoxia score showed poor disease-free survival (DFS) compared to patients with a low score (P = 0.0021). To further investigate the association between gene expression and hypoxia scores, we performed a differential expression analysis between high- and low-hypoxia scores. A total of 69 DEGs, including 63 upregulated and 6 downregulated genes, were obtained (Figures 1C and 1D, Table S2). These DE-HRGs were involved in the HIF-1 signalling pathway, various cancers, and several immune-related pathways in the progression of PCa (Figures 1E and 1F).
Identification of immune-related DEGs (DE-IRGs) in PCa
Furthermore, we determined the immune scores of 475 PCa samples with complete survival data using the ESTIMATE algorithm. As a result, the optimal cut-off value of the immune score was -103.9, which divided the PCa cases into high- and low-immune score groups (Figure 1G, Table S3). To establish a possible correlation between the immune score and DFS, we conducted a Kaplan-Meier (K-M) analysis between the high- and low-immune score groups. The results showed that a high-immune score was significantly associated with a poorer outcome (P = 0.01) (Figure 1H). In addition, a total of 981 DEGs, including 946 upregulated genes and 35 downregulated genes, were identified between the high- and low-immune score groups (Figures 1I and 1J, Table S4). These DEGs were associated with immune-related biological processes and signalling pathways, such as T cell activation, MHC protein complex, chemokine receptor binding, chemokine signalling pathways, Th17 cell differentiation, and NF-kappa B signalling pathway (Figures 1K and 1L).
Identification of hypoxia- and immune-related DEGs for PCa
In total, 28 hypoxia- and immune-related DEGs were obtained using Venn analysis (Figure 2A). To better understand the biological roles of the 28 hypoxia- and immune-related DEGs, the clusterProfilerR package was employed to perform GO annotation and KEGG pathway enrichment analyses. Figures 2B and 2C illustrate the top 20 enriched BPs in terms of GO and KEGG pathways. GO analysis revealed that these DEGs were markedly enriched in the regulation of leukocyte migration and chemotaxis, and the regulation of vasculature development. In addition, the significantly enriched pathways for these 28 DEGs were TNF signalling pathway, colorectal cancer, AGE-RAGE signalling pathway in diabetic complications, malaria, p53 signalling pathway, Chagas disease, osteoclast differentiation, and IL-17 signalling pathway.
Development of the hypoxia- and immune-related risk signature in PCa
In the training set, five hypoxia- and immune-related genes were identified through univariate Cox regression analysis, which were significantly associated with DFS of patients with PCa (P < 0.005) (Figure 2D, Table S5). Two genes, ISG15 and ZFP36, were selected through stepwise multivariate regression analysis (Figure 2E, Table S6), and used to develop a prognostic gene signature. The K-M curves of the two hypoxia- and immune-related genes are shown in Figure 2F.
The risk score of each patient with PCa was calculated using the formula mentioned above, and the patients were categorised into high- and low-risk groups according to the median risk score. The distribution of risk scores for each patient and the expression of the two prognostic genes are presented in Figures 2G and 2H, respectively. Patients with PCa and a high-risk score had a shorter survival time than those with a low-risk score (Figure 2I) (P = 0.0053). The predictive accuracy of the two hypoxia- and immune-related signatures was assessed using time-dependent ROC curves. The area under the curve (AUC) scores of ROC curves at 1, 3, and 5-year were 0.674, 0.649, and 0.605, respectively, indicating that the prognostic signature had a good predictive performance (Figure 2J).
Validation of the hypoxia- and immune-related risk signature in MSKCC database
To further demonstrate the predictive accuracy of the hypoxia- and immune-related risk signature, an independent external dataset, the Memorial Sloan Kettering Cancer Centre (MSKCC database), was applied to validate the results obtained from the training set. The K-M curves of the two hypoxia- and immune-related genes in MSKCC are shown in Figure 3A. In the validation dataset, we divided patients with PCa into high- and low-risk groups based on the median risk score calculated with the same equation. The risk score distribution and gene expression patterns for patients in the validation set are shown in Figures 3B and 3C. The patients in the high-risk group had markedly shorter DFS than those in the low-risk group (Figure 3D, P = 0.00018), which was consistent with the results of the training dataset. The AUCs of the hypoxia- and immune-related risk signature were 0.720, 0.720, and 0.747 at 1-, 3-, and 5-years, respectively. This further suggested a significant predictive value of the prognostic signature for clinical outcomes (Figure 3E).
Identification of independent prognostic indicators for PCa
To further explore the relationship between the risk score and clinical outcomes, we analysed the differences in clinical features of the high- and low-risk groups in the training and validation sets, respectively. As illustrated in Tables 1 and 2, the risk score was associated with the Gleason score, T stage, and N stage of patients with PCa (P < 0.05). Based on these clinical features, univariate and multivariate Cox regression analyses were conducted to identify independent predictors of PCa prognosis. Gleason score and risk score were independent prognostic factors for PCa (Figures 4A and 4B). Gleason score and risk score were integrated to generate a nomogram in the training set (Figure 4C). In addition, the calibration curves and the nomogram showed good agreement, which provided a quantitative method for clinicians to predict the survival probability of 1-, 3-, and 5-year survival rates (Figure 4D).
Functional enrichment analysis of prognostic hypoxia- and immune-related genes
As shown in Figure 5A, the expression of the two prognostic hypoxia- and immune-related genes was compared between tumour tissues and normal tissues. Compared with normal tissues, the expression of ISG15 was significantly increased in tumour tissues, whereas the expression of ZFP36 was significantly decreased. To further explore the functions of the two prognostic genes in PCa progression, GSEA was employed to search for pathways enriched in TCGA samples. A total of 72 KEGG pathways were enriched, from which 46 were enriched in genes positively correlated with ISG15, and 26 were enriched in genes negatively correlated with ISG15 (Table S7). The top 10 pathways are visually displayed in Figure 5B. Amyotrophic lateral sclerosis, Epstein-Barr virus infection, hepatitis C, Huntington’s disease, influenza A, oxidative phosphorylation, Parkinson’s disease, prion disease, proteasome, and systemic lupus erythematosus were differentially enriched in ISG15 high expression samples. In addition, 197 KEGG pathways were significantly enriched, from which 151 were enriched in genes that were positively related to ZFP36, and 46 were enriched in genes that were negatively related to ZFP36 (Table S8). Figure 5C illustrates that upregulation of ZFP36 was associated with viral protein interaction with cytokine and cytokine receptor, NF-kappa B signalling pathway, cytokine-cytokine receptor interaction, MAPK signalling pathway, chemokine signalling pathway, and calcium signalling pathway. However, downregulation of ZFP36 was correlated with RNA transport, aminoacyl-tRNA biosynthesis, oxidative phosphorylation, and ribosomes.
The risk score was associated with immune cell infiltration
The immune cell infiltration level of each PCa sample in the TCGA database was evaluated using the ssGSEA algorithm. The enrichment score of 24 immune-related gene sets for each PCa sample was quantified to estimate the abundance of immune cells in the tumour immune microenvironment. Subsequently, we analysed the differences in immune infiltration levels between the high- and low-risk groups. The infiltration levels of eosinophils, neutrophils, Tcm, Tem, TFH, Th1 cells, and Th17 cells were significantly lower in the high-risk group than in the low-risk group. Conversely, aDC, pDC, T helper cells, and Treg cells were significantly higher in the high-risk group (Figure 6A) (P < 0.05). In general, the results indicated that the risk score was associated with immune cell infiltration in PCa; thus, affecting cancer progression. The correlation between the 24 immune cells is shown in Figure 6B.
To further explore the relationship between the two prognostic hypoxia- and immune-related genes and the infiltration of immune cells, we conducted Pearson’s correlation test to calculate and analyse the correlation coefficients of the two genes in the 24 immune-related subsets (Figure 6C). We found that the two genes were closely correlated with the tumour-infiltrating immune cell subset, revealing that the prognostic genes mainly participate in the immune response during PCa progression (P < 0.05).
Expression of ISG15 and ZFP36 is increased and decreased in PCa respectively and ISG15 is positively correlated with hypoxia
We demonstrated the differential expression of ISG15 and ZFP36 in PCa and paracancerous tissue by immunohistochemical staining and in PCa cell lines and normal prostate epithelial cell lines by real-time quantitative PCR (RT-qPCR) and western blot (WB). Immunohistochemical results of PCa and paracancerous tissues showed differences in the expression of ISG15 and ZFP36, respectively(Figure 7A). As shown in Figure 7B-C, the expression level of ISG15 was significantly higher in PC-3, 22RV-1, and LNCaP than in RWPE-1. ZFP36 showed opposite results. After hypoxic treatment of 22RV-1, PC-3, and LNCaP cells for 24h, 48h, 72h, respectively, we performed RT-qPCR and WB to demonstrate that the expression level of ISG15 and ZFP36 was positively and negatively correlated with hypoxia, respectively. (Figure 7D-K).
Overexpression of ISG15 promote migration and invasion ability of PCa cells
After verifying that ISG15 expression levels were positively regulated by hypoxia, we performed wound-healing and transwell invasion assays using 22RV-1 and PC-3 cells after transfection (Figure 8A-D). The results showed that a high expression of ISG15 could promote the migration and invasion ability of PCa cells, thereby suggesting that high expression of ISG15, under hypoxic conditions, may promote poor prognosis of PCa.