Comprehensive analysis of integrin αvβ3/α6β1 in prognosis and immune escape of prostate cancer

Integrin αvβ3/α6β1 were closely associated with poor immunity in PCa patients

First, we systematically analyzed the correlations among integrin αvβ3/α6β1 with immune cell infiltration, immune functions and immunosuppressive genes. Interestingly, we found that integrin αvβ3/α6β1 was negatively associated with immune functions, including APC co-stimulation, cytolytic activity, T-cell co-stimulation and the type II IFN response, and integrin αvβ3/α6β1 was negatively correlated with most immune functional cells, including DCs, B cells, NK cells, CD8⁺ T cells, mast cells, Th cells, macrophages, and other tumor infiltrating lymphocytes (TILs) (Figure 1A). In contrast, integrin αvβ3/α6β1 was positively correlated with negative immune function, including APC co-inhibition and Tregs. Moreover, integrin αvβ3/α6β1 was negatively correlated with memory B cells, activated CD4+ memory T cells and activated mast cells, which could reduce immune cellular functions in the tumor microenvironment (TME) of PCa patients (Figure 1B). Low infiltration of immune cells in the TME and immune cell dysfunction may contribute to the immune escape of PCa cells. In addition, it also showed close positive correlations when intersecting integrin αvβ3/α6β1 with some classical immunosuppressive genes that were connected with immune escape, such as CTLA4, IL10RB, PDCD1LG2 (B7-DC), CD274 (B7-H1), TGFBR1, and PDCD1 (Figure 1C). Similarly, the protein-protein interaction network (PPI) showed tight correlations among these genes with integrin αvβ3/α6β1, which indicated mutually regulated relationships between integrin αvβ3/α6β1 and those immunosuppressive genes or their expressed proteins (Figure 1D). Together, these findings preliminarily suggested that integrin αvβ3/α6β1 might be strongly correlated with poor immunity and immune escape in PCa.

Integrin αvβ3/α6β1 were closely associated with poor prognosis of PCa patients

In our study, a heatmap was established to clearly observe the expression levels of integrin αvβ3/α6β1 and clinical characteristics, including PSA value, primary therapy outcome success, pathological T and N stages, Gleason score, treatment success and clinical T and N stages, of 496 PCa patients in the TCGA cohort, and the results implied that integrins ITGAV, ITGA6, ITGB1 and ITGB3 were closely associated with these clinical features (Figure 2A). Additionally, we found higher expression levels of these integrins in low-grade PCa, especially in medium-grade and high-grade PCa, than in normal prostate tissues, according to the immunohistochemical results from the Human Protein Atlas (HPA) (Figure 2B). To explore the effects of integrin αvβ3/α6β1 on the prognosis of PCa patients, we performed Kaplan-Meier analysis for overall survival (OS), relapse-free survival (RFS) and biochemical recurrence based on the expression levels of integrin ITGAV, ITGA6, ITGB1 and ITGB3 in these samples. The results showed that the OS for the samples with high expression levels of ITGB1 was significantly lower than that for the samples with low expression levels of ITGB1 (Figure 2C). Likewise, the RFS for the samples with high expression levels of ITGAV and ITGB1 was markedly lower than that for the samples with low expression levels of ITGAV and ITGB1, which might indicate that integrin αvβ3/α6β1 could be a potential prognostic element in PCa (Figure 2D, 2E). Furthermore, PCa patients who expressed high levels of integrin ITGAV, ITGA6, ITGB1 and ITGB3 had a higher proportion of biochemical recurrence than patients with low expression levels of integrin αvβ3/α6β1 (Figure 2F–2I). Moreover, we further investigated the correlations among the expression levels of the integrins ITGAV, ITGA6, ITGB1 and ITGB3 and the expression levels of AR in PCa patients. According to literature reports, androgen receptor (AR) is one of the prognostic factors in PCa, and contributes to the progression and recurrence of PCa. Surprisingly, there were strong correlations among the integrins ITGAV, ITGA6, ITGB1 and ITGB3 with AR in PCa patients, whose expression levels for cor-values were 0.759, 0.650, 0.575, and 0.448, respectively (Figure 2J–2M). In summary, these results suggested that integrin αvβ3/α6β1 could be a negative prognostic factor for PCa.

A prognostic risk model based on integrin αvβ3/α6β1 for PCa patients

First, we divided the PCa patients into two clusters (C1 and C2) according to the average expression levels of whole integrin αvβ3/α6β1 (ITGAV, ITGA6, ITGB1, ITGB3) in the TCGA cohort, and we performed gene cluster analysis for the two clusters to identify significant DEGs based on the whole genetic transcriptomes of 496 PCa patients in the TCGA dataset, and 13994 DEGs were identified in the two clusters (Figure 3A). In addition, information on the clinical characteristics of these PCa patients was also manifested in the heatmap. Next, 51 DEGs were identified as prominently affecting the biochemical recurrence of PCa by conducting univariate Cox analysis, including 33 DEGs (hazard ratio >1) that increased the risk of biochemical recurrence and 18 DEGs (hazard ratio <1) that decreased this risk in the TCGA cohort (Figure 3B). Subsequently, we intersected these 51 DEGs with the DEGs acquired from the GEO dataset and successfully identified eight biochemical recurrence-related genes (ASF1B, INSM2, POU4F2, MT1B, NCR1, KRTAP10-5, PCDHA13 and KIR3DL1) in both the TCGA and GEO cohorts. Next, we constructed a risk model using the optimum γ value based on the eight genes by applying LASSO Cox regression analysis (Figure 3C, 3D). In addition, we calculated the risk score of these 8 biochemical recurrence-related genes by running the formula: Risk score= (0.0741808537008443 * expression of ASF1B) + (0.0584797685631388 * expression of INSM2) + (-0.0665901093358763 * expression of POU4F2) + (0.152621619239078 * expression of MT1B) + (0.00167204457436465 * expression of NCR1) + (-0.0181031499555924 * expression of KRTAP10-5) + (0.151620509953615 * expression of PCDHA13) + (-0.0816400681584183 * expression of KIR3DL1).

Based on the median risk score, PCa patients in the TCGA cohort and GEO cohort, which had complete biochemical recurrence information, were objectively divided into the HR group and LR group (Supplementary Figure 1A). The genes ASF1B, INSM2, MT1B, NCR1 and PCDHA13 were highly expressed in the HR group, while the genes POU4F2, KRTAP10-5 and KIR3DL1 were highly expressed in the LR group in both the TCGA and GEO cohorts (Figure 3E, 3F). To further confirm the relationship between these genes and biochemical recurrence, we selected 5 cases of biochemically recurrent prostate cancer and 5 cases of nonbiochemically recurrent prostate cancer for RNA extraction. The results showed that the expression levels of ASF1B, INSM2, MT1B, NCR1 and PCDHA13 in biochemically recurrent PCa were higher than those in nonbiochemically recurrent PCa, while POU4F2, KRTAP10-5 and KIR3DL1 were significantly lower in biochemically recurrent PCa than in nonbiochemically recurrent PCa (Figure 3G). Additionally, we were able to find that there were more dead samples that were regarded as the HR group located in the area of higher risk scores with increasing risk scores of patients. In contrast, the other survivors who were regarded as the LR group were more likely to be located in areas with lower risk scores, and they had longer survival times than the HR group (Supplementary Figure 1B). Otherwise, we applied principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) analysis to demonstrate that samples in the HR and LR groups could be greatly separated (Supplementary Figure 1C, 1D).

Furthermore, we carried out Kaplan-Meier analysis for biochemical recurrence between the two groups, and the results revealed worse biochemical recurrence probabilities for PCa patients in the HR group than in the LR group (Figure 4A). ROC analysis was performed to prove the predictive ability of our model for biochemical recurrence, and the area under the ROC curve (AUC) was 0.716, 0.764, and 0.832 for 1, 3, and 5 years, respectively (Figure 4D). Likewise, worse progression-free survival (PFS) and OS were observed in the HR group than in the LR group (Figure 4B, 4C). The AUC values for PFS were 0.950 at 1 year, 0.906 at 3 years, and 0.869 at 5 years, while the AUC values for OS were 0.992 at 1 year, 0.906 at 3 years, and 0.912 at 5 years (Figure 4E, 4F). Moreover, univariate and multivariate analyses showed that the risk score was an independent risk factor for PCa patients (Figure 4G, 4H). These evidences illustrate that our risk model was powerful for predicting the prognosis of PCa patients.

The HR risk group predicted poor clinical features of PCa patients

On the other hand, we investigated the correlations of the two risk groups with clinical features in PCa patients. The results showed that the percent weight of clinical and pathologic T3-stage and T4-stage in the HR group was higher than that in the LR group (40% vs 15% and 95% vs 80%), and the risk score of clinical and pathologic T3-stage and T4-stage was higher than that of clinical and pathologic T1-stage and T2-stage, respectively (Figure 5A–5D), suggesting that more patients suffered advanced PCa in the HR group. In addition, the percent weight of fustat-1 (the dead patients) in the HR group was markedly higher than that in the LR group (74% vs 41%), which was the same result as the risk score in the two groups (Figure 5E, 5F). In addition, the percent weight of PSA (value >4), which could greatly identify PCa in patients, was dramatically higher in the HR group than in the LR group (15% vs 2%), and the risk score of PSA (value >4) was higher than that of PSA (value <4) in the two groups as well (Figure 5G, 5H).

Moreover, Kaplan-Meier survival analysis was performed for the two risk groups under subgroups of clinical characteristics of PCa patients. The results showed that PCa patients with clinical T1-stage, T2-stage and T3-stage disease probably suffered more mortality in the HR group than in the LR group (Figure 5I–5K). Similarly, the patients with Gleason score 7-9, PSA (value >4 and <4), pathologic stages of T2-T3 and N0-N1 both suffered lower survival probability in the HR group than in the LR group (Supplement Figure 3A–3I). This evidence strongly demonstrated that the HR group could well predict poor survival in PCa patients.

Genetic mutation status and functional enrichment analysis for the two risk groups

In the TCGA cohort, the GO enrichment analysis showed that GOBP: NIK (NF-κB induced kinase) NF-κB signaling, GOBP: ras protein signal transduction, GOBP: canonical Wnt signaling pathway and GOBP: positive regulation of interleukin-10 production was markedly enriched in the HR group compared with the LR group. Moreover, according to the GEO cohort, GOBP: transforming growth factor beta1 production and GOMF: heat shock protein binding was more enriched in the HR group than in the LR group. Interleukin-1-mediated signaling and FC receptor signaling were more significantly enriched in the HR group than in the LR group in both the TCGA and GEO cohorts (Figure 6A, 6B and Supplementary Tables 1, 2). In addition, KEGG enrichment analysis showed that KEGG: T-cell and B-cell receptor signaling pathway, KEGG: P53 signaling pathway, KEGG: wnt signaling pathway, TGF beta signaling pathway, chemokine signaling pathway, KEGG: cell cycle and KEGG: MAPK signaling pathway were highly enriched in the HR group compared with the LR group in the TCGA or GEO cohorts (Figure 6C, 6D and Supplementary Tables 3, 4). Based on these results, we can note that most of the functions or pathways enriched in the HR group are associated with poor patient outcomes. Therefore, these evidences might reveal the potential mechanisms for poor prognosis in the HR group.

To evaluate the difference in genetic mutation status in the two risk groups, we used Maftools to identify the whole gene mutation in 91 PCa patients with complete biochemical recurrence data from the TCGA dataset, of which the 20 most frequently mutated genes are shown in Figure 7A, 7B. We noticed a higher frequently mutated ratio in the HR group than in the LR group (75.56% vs 58.70%). Moreover, the TP53 and KMT2D genes were more often mutated in the HR group than in the LR group (18% vs 15% and 9% vs 7%). However, the gene SPOP mutation rates were higher in the LR group than in the HR group (20% vs 9%), which could be attributed to genetic heterogeneity. Collectively, these mutated genes might partially contribute to the poor clinical prognosis of PCa in the HR group.

The HR group predicted a higher risk of immune escape for PCa patients

Moreover, we further investigated the components of immune cells that infiltrated the TME of PCa between the two groups in TCGA samples (Figure 7C). The results showed lower infiltration of DCs, B cells, CD8⁺ T cells, macrophages, NK cells, Th cells and TILs in the HR group than in the LR group. In contrast, negatively regulated immune cells, such as Treg cells, were higher in the HR group than in the LR group. To further investigate the influences of immune cells on our patients, Kaplan-Meier analysis was performed for PFS, OS and biochemical recurrence in the TCGA cohort. As expected, the high expression levels of CD8⁺ T cells, NK cells and macrophages exhibited better PFS, OS or biochemical recurrence, while the high expression levels of Tregs exhibited poor PFS and OS (Supplementary Figure 2A–2H). These results were consistent with conventional wisdom.

Additionally, immune functions were also investigated in the HR and LR groups, and the results showed that poor immune cell functions, including aDCs, B cells, CD8+ T cells, NK cells and macrophages, and better Treg cell functions were found in the HR group than in the LR group. Immune cell functions such as APC co-stimulation, T-cell co-stimulation and the type II IFN response were also lower in the HR group than in the LR group (Figure 7D). These evidences illustrate worse immune functions in the HR group, which might contribute to tumor immune escape in PCa. Next, we compared the expression of immune checkpoints in the risk groups and found that the HR group had higher expression of PD-1 (CD274), PD-L1 (PDCD1), and CTLA-4 than the LR group, which might lead to a higher risk of immune escape in PCa (Figure 7F–7H). In addition, we scored the risk of immune exclusion in our patients; surprisingly, we found that the immune exclusion score was higher in the HR group than in the LR group (Figure 7E). These findings suggested that the HR risk group might suffer a higher risk of immune escape in PCa patients.

Data sources

The data of genetic transcriptomes, tumor mutational burden (TMB), SNPs and clinical characteristics for 495 PCa samples were acquired from the TCGA dataset (https://tcga-data.nci.nih.gov/tcga/), and the genetic transcriptomes and clinical features for 220 PCa cases were downloaded from the GEO dataset (http://www.ncbi.nlm.nih.gov/geo/). The results of immunohistochemistry (IHC) of prostate tissues and PCa were acquired from the Human Protein Atlas (HPA) dataset (https://www.proteinatlas.org/). The expression level of integrin αvβ3/α6β1 in each clinical PCa patient was determined by the mean of the expression levels of the integrin subunits ITGAV, ITGA6, ITGB1, and ITGB3.

Immune correlation analysis and protein-protein interaction network

We utilized CIBERSORT to acquire information on immune infiltrating cells and immune functions in PCa patients. Pearson correlation analysis was performed for integrin αvβ3/α6β1, immune cells, immune functions and immunosuppressive genes in PCa patients. In addition, the protein-protein interaction network (PPI) for integrin αvβ3/α6β1 with immunosuppressive genes was formed by applying STRING (https://string-db.org/) and Cytoscape software.

Identification of differentially expressed genes

First, the gene expression data in both the TCGA and GEO cohorts were standardized with the formula log2(x+1), and we intersected the TCGA genetic matrix with the GEO genetic matrix to acquire the expression data of the same genes in both the TCGA and GEO cohorts. Then, the data of differentially expressed genes (DEGs) were finally identified by applying the edge R package with FDR < 0.05 and |log2FC| ≥ 0.585 in both TCGA and GEO cohorts.

Establishment of the prognostic risk model

We intersected biochemical recurrence-related genes in both the TCGA and GEO databases, and eight biochemical recurrence-related genes were finally identified, which were used to establish our prognostic risk model by applying LASSO Cox regression analysis. In our work, the risk score was calculated by ${\sum }_{\text{i}}^8\text{Ai}\times \text{Bi}$ (A: coefficients, B: gene expression level), which was utilized to separate our samples into two groups (HR and LR group) in both the TCGA training cohort and GEO testing cohort.

Functional enrichment analysis

In this work, we first performed gene differential analysis of the HR group and LR group and identified 10407 and 1922 significantly different genes in the TCGA dataset and GEO dataset, respectively. These genes were used in the Database for Annotation, Visualization and Integrated Discovery (DAVID) (https://david.ncifcrf.gov/) to run gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses for the two risk groups. For gene set enrichment analysis (GSEA), TCGA and GEO transcriptome data and risk files were input into R software. The program performs enrichment analysis by sorting the expression profile data, calculating enrichment scores, estimating significance levels, and performing multiple hypothesis testing.

Statistical analysis

We employed one-way ANOVA and t tests to compare two risk groups, and comparisons of two or more component percent rates were made using the chi-square test. In addition, heatmaps, forest graphs, box plots and receiver operating characteristic (ROC) curves were completed by R software (version 3.5.1). Waterfall curves were completed by MAF tools and R software. Otherwise, Kaplan-Meier analysis was performed to estimate the survival and biochemical recurrence of PCa patients in both the TCGA cluster and two risk groups in the model. We processed all statistical analyses by utilizing SPSS 19.0 software (SPSS. Inc., Chicago, IL, USA) or R software. In our study, a P-value < 0.05 was statistically significant.

Data availability statement

All data in this study are available in the TCGA data portal (https://tcga-data.nci.nih.gov/tcga/) and GEO data portal (http://www.ncbi.nlm.nih.gov/geo/).