Identification of biomarkers for immunotherapy response in prostate cancer and potential drugs to alleviate immunosuppression

Prediction of immune response status by tumor immune dysfunction and exclusion (TIDE)

The mRNA expression profile of PRCA patients was obtained from the Gene Expression Omnibus (GEO) database (GSE183019). We predicted the responses of 84 patients to immunotherapy and estimated the scores of 10 immune features with TIDE. Overall, 58 patients were predicted to be responsive to immunotherapy as their TIDE score was <0, whereas the other 26 patients were considered non-responders (Figure 2A). All 58 responder patients were CTL true (positive for five cytotoxic T lymphocyte markers, including CD8A, CD8B, GZMA, GZMB, and PRF1), indicating that the immune response to tumor T cell infiltration is highly consistent; that is, patients with high T cell infiltration have a high probability of responding to immunotherapy. Moreover, interferon gamma (IFNG) was a positive biomarker for the immune checkpoint blockade therapy response, and the score of IFNG in responder samples was significantly higher than that in non-responders. Except for the cancer-associated fibroblast (CAF) score, all immune feature scores showed significant differences between the groups. In addition, Merck18 (T-cell-inflamed signature), IFNG, CD8, and dysfunction signatures showed the most significant differences between the immune response and non-response groups (Figure 2B, 2C).

Identification of immune response-related genes through differential expression analysis (DEA)

By comparing the transcriptomes of 58 immune responder patients and 26 non-responder patients, 2758 immune response-related genes (IRRGs) were identified, and 99.42% (2742/2758) were expressed at low levels in immune-responsive samples (Figure 3A). The three most significant differentially expressed genes were HBA2, LOC100131257, and CCDC168. Notably, it has been reported that mutations in CCDC168 are associated with adenosquamous carcinoma of the prostate [9]. The top 30 differentially expressed genes included multiple non-coding RNAs, such as LINC00907, LINC01105, LINC00276, SNORA57, and SNORA76C (Figure 3B); these non-coding RNAs can regulate gene expression through post-transcriptional modifications, thereby affecting the immune microenvironment of the prostate.

Functional enrichment analysis of IRRGs

To explore the potential functional implications of IRRGs, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses. As a result, the most significantly enriched pathways in KEGG were “olfactory transduction,” “cAMP signaling pathway,” and “calcium signaling pathway.” Moreover, “nicotine addiction” and “cytokine-cytokine receptor interaction” were also enriched in the KEGG pathways (Figure 3C). Similarly, GO enrichment analysis suggested that IRRGs were significantly enriched in material transport processes, such as “potassium ion transport,” “regulation of membrane transport,” and “anion transmembrane transport” (Figure 3D). Various signaling processes were also enriched, including “glutamate receptor signaling pathway” and “coupled receptor signaling pathway.” Functional annotation of IRRGs revealed that ion transport and signaling pathways are implicated in patient response to immunotherapy, thus contributing to anti-tumor immunosuppression.

IRRG co-expression module identification with weighted gene co-expression network analysis (WGCNA)

In total, 2758 differentially expressed IRRGs were involved in the co-expression modules. First, the soft threshold was determined through the scale independence and mean connectivity analysis of modules with different power values ranging from 1 to 20 (Supplementary Figure 1). In this study, the power value (β) was set to 6 to produce a hierarchical clustering tree with different colors representing different modules. As a result, 16 total modules with different IRRGs were identified and displayed with different colors.

The hierarchical clustering dendrogram of the patients based on the WGCNA distance matrix is shown in Figure 4A, and the immune responding patients were grouped together. IRRGs were grouped into 16 co-expression modules (Figure 4B), and the correlation between modules and immune features is displayed in Figure 4C, which suggested that turquoise module was significantly correlated with Merck18, IFNG, and myeloid-derived suppressor cells (MDSCs). Seven modules that were significantly associated with TIDE score, with p-value < 0.05, were retained for further analysis.

To provide insight into the biological functions of the IRRGs in seven immune response-related modules, we performed GO (Supplementary Table 1) and KEGG (Supplementary Table 2) enrichment analysis. Most modules included enriched GO terms related to voltage-gated channels, such as “cation channel complex” (3/7), “ion transport” (5/7), and “sensory perception of smell” (4/7). Specifically, the IRRGs involved in the red module were significantly enriched in immune-related GO terms, such as “lymphocyte/T cell costimulation” “, immunological synapse ” “, positive regulation of adaptive immune response”, and “cytokine metabolic process” (Figure 5A). IRRGs of turquoise module were enriched in the “regulation of uterine smooth muscle contraction” (Figure 5B). Further, IRRGs of green (Figure 5C), magenta (Figure 5D), and cyan (Figure 5E) modules were mostly enriched in pathways related to “ion transport” and those of magenta module were especially enriched in “epidermal cell differentiation”. IRRGs of pink (Figure 5F) and salmon (Figure 5G) module were mostly enriched in “sensory perception of smell” and “olfactory receptor activity”, respectively.

Identification of biomarkers of immunotherapy responses

To investigate the role of IRRGs as potential marker genes in immunotherapy response, we first identified hub genes in seven modules of interest (see Methods). Ultimately, 133 genes were identified as hub genes, which were selected for survival analysis (Table 1), and 13 of them were found to be significantly associated with PRCA prognosis. BICDL2, a protein-coding gene predicted to enable small GTPase binding activity, showed the most significant association with PRCA prognosis (p-value = 0.004). Moreover, the rest of the prognosis-related genes were as follows: ABHD17A (p-value = 0.0083), ARHGAP33 (p-value = 0.0093), AP5Z1 (p-value = 0.041), ARHGEF16 (p-value = 0.026), ASMTL-AS1 (p-value = 0.029), ATAD3B (p-value = 0,029), ATG16L2 (p-value = 0.014), AXIN1 (p-value = 0.042), LLCFC1 (p-value = 0.016), LENG8 (p-value = 0.018), KMT2B (p-value = 0.05), and CACNA1H (p-value = 0.041) (Figure 6). The IRRGs described previously herein were considered biomarkers of immunotherapy response in PRCA, which might affect the immune microenvironment, thus playing an important role in patient responses to immunotherapy and determining clinical outcomes.

Construction of co-expression network

Using WGCNA, we constructed co-expression networks of the hub genes in the seven modules identified previously herein based on expression correlations to provide insights into the role of signature genes in PRCA pathogenesis. Excluding the red module, which had only one hub gene, we constructed networks for the other six modules. The results showed that in the turquoise module, INTS1 and KMT2B were at the center of the network. KMT2B was significantly related to survival and was related to multiple regulators of gene expression (Figure 7A), such as LOC100133091 and SNORD17, as well as the transcription factor (TF)-encoding genes CAPN15, SNAPC4, and CASZ1; these regulators themselves act as modular hub genes and might have regulatory functions. There were four lncRNAs and one snoRNA involved in the network of cyan module (Figure 7B). In green module, the co-expression relationship among 13 hub genes was presented (Figure 7C) and nine hub genes, including two lncRNAs and one TF-encoding gene were used in magenta module (Figure 7D). For hub genes in salmon and pink modules, the networks were constructed and presented in Figure 7E and 7F, respectively.

Identification of small molecules targeting core IRRGs

Having identified PRCA signature genes significantly associated with the immunotherapy response, most of which were downregulated in immune-responsive patients, we hypothesized that when signature genes are targeted with small molecules, the downregulated expression of these genes would sensitize cells to the drug. Such molecules could potentially reverse tolerance and facilitate therapy effectiveness when combined with traditional drugs. Data from the Library of Integrated Network-Based Cellular Signatures (LINCS) database were used to identify small molecules targeting PRCA signature genes, with their co-expressed genes were identified by WGCNA.

Based on the signature genes and the co-expression pairs identified by WGCNA, from hundreds of small molecules, we identified nine and two candidate drugs, respectively, from the expression profiles of the PC3 and LNCaP cell lines from the LINCS database (Figure 8A, Supplementary Table 3). These drugs regulated the expression of signature genes of PRCA significantly. Among them, four and two compounds, respectively, were considered for immunotherapy based on PC3 and LNCaP cells, and these significantly inhibited the expression of signature genes, making them potential candidates for combination therapies. For PC3 cells, resveratrol most significantly inhibited the expression of signature genes (Figure 8B), whereas in the LNCaP cell line, radicicol (Figure 8C) and PHA-793887 had the most significant effects. No compound was observed to have significant effectiveness in both cell lines.

Data source

The mRNA expression profiles and dataset comprising clinical information of PRCA patients of GSE183019 were acquired from the GEO (https://www.ncbi.nlm.nih.gov/geo/) database [35, 36] and The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov) database [37]. Excluding normal tissue samples, 84 PRCA samples with count profiles of 25498 expressed genes were retained for further analysis. The transcriptome profiles of the PC3 and LNCaP cell lines with different drug perturbations were obtained from the LINCS (https://lincsproject.org/LINCS/) database [38], and data of a total of 1121 compounds with the expression of 7163 genes were collected.

Prediction of immune response status

TIDE [39] is an integrated web tool used to predict patient responses to immunotherapy. Based on tumor expression profiles, TIDE can score multiple transcriptomic biomarkers of several immune features, including immune system dysfunction, T cell exclusion, MDSCs, IFNG signature, CAFs, M2 subtype of tumor-associated macrophages (TAMs), microsatellite instability (MSI), Merck18, CD274, and CD8. These features were demonstrated to be related to patient responses to immunotherapy, and we calculated the estimated score of the immune features for each patient based on their mRNA expression profiles. Specifically, patients with a TIDE score < 0 were considered responders, whereas patients with a TIDE score > 0 were considered non-responders.

Identification of immune response-related genes

Based on the immune response status predicted by TIDE, differentially expressed genes were calculated between responder and non-responder patient groups using DEA with the R package DESeq2 [40] in R software. The significantly differentially expressed genes were identified using the criteria of |log₂ fold-change| > 1 and adjusted p-value < 0.05. Genes meeting the criteria were considered IRRGs.

Functional annotation of immune response-related genes

To explore the potential biological functions and pathways associated with IRRGs, we performed GO (http://geneontology.org/) [41] and KEGG (https://www.genome.jp/kegg/) [42] pathway analysis using the R package “ClusterProfiler” [43]. GO terms and KEGG pathways with a p-value < 0.05 were considered significantly enriched.

Identification of co-expression modules for immune-related genes

IRRGs were used to construct co-expression modules with the WGCNA [44] package in R software. Pearson’s correlation matrices were calculated between each pair of retained IRRGs from the corresponding expression levels. Clusters were obtained from the dendrogram by applying the dynamic tree-cutting technique [45]. Module-trait relationships were calculated according to the correlation between modules and the predicted score of each immune feature based on Pearson’s correlation tests, and TIDE scores were selected as core features to identify significant co-expression modules related to the immune response. Ultimately, IRRGs in significant modules (p < 0.05) were exported for further analysis.

Identification of biomarkers of immunotherapy response

For co-expression modules constructed with WGCNA, hub genes of the modules of interest were identified first. For genes inside a given module, the within-module connectivity, called the module membership (MM), and the correlation with the TIDE score, called the gene significance trait (GS), were calculated as previously described [46]. Genes with high connectivity in the modules of significant immune features were identified as hub genes based on the cut-off criteria |MM| > 0.8 and |GS| > 0.2.

The expression profiles and related clinical information of PRCA samples were acquired from TCGA database. Survival analysis based on the Cox model was used to estimate the survival risk of patients under different conditions. For a specific hub gene, patients were divided into two groups according to median gene expression. Then, Kaplan–Meier curve analysis was performed with the R package “Survival” and the p-value between the two groups was also calculated. Hub genes (p-value < 0.05) that were significantly related to PRCA prognosis were considered biomarkers of the immunotherapy response.

Construction of mRNA interaction network

To better investigate the regulatory relationships between IRRGs, we constructed an IRRG co-expression network for hub genes identified in each module of interest. Significantly correlated pairs were used to construct the network based on Pearson’s correlation coefficients. Finally, the co-expression network graphs were visualized and analyzed using Cytoscape software [47].

Identification of candidate drugs from LINCS

To identify potential drugs to boost immunotherapy, we employed the method reported by Aissa et al. [48]. First, the mRNA expression profiles of hundreds of perturbations at various time points in PRCA cell lines were collected from LINCS. For each drug, different treatment groups, such as 6 h vs. 24 h, were compared to calculate the differentially expressed genes. According to the fold-change value associated with each gene, differentially expressed genes were ranked using Gene Set Enrichment Analysis (GSEA) [49]. Enrichment scores were calculated, and drugs with a p-value < 0.05 were considered significantly enriched. Specifically, a positive enrichment score shows that the drug has an upregulation effect on the target gene set, and a negative enrichment score indicates that the drug has a downregulation effect on the target gene set.