Data-driven analysis to identify prognostic immune-related biomarkers in BRAF mutated cutaneous melanoma microenvironment

Skin cutaneous melanoma is one of the deadly diseases, and more than 50% of the patients have BRAF gene mutations. Evidence suggests that oncogenic BRAF modulates the immune system’s ability to recognize SKCM cells. Due to the complexity of the tumor microenvironment (TME) and a lack of a rational mechanistic basis, it is urgent to investigate the immune infiltration and identify prognostic biomarkers in BRAF mutated SKCM patients. Multiple methods including ESTIMATE algorithm, differential gene analysis, prognostic analysis and immune infiltration analysis were performed to investigate the tumor microenvironment. Based on the patient’s immune score and stromal score, immune-related genes DEGs were identified. Functional analysis revealed that these genes were mainly enriched in biological processes such as immune response, defense response and positive regulation of immune system. Furthermore, we analyzed the immune infiltrating cell components of BRAF mutated patients and revealed 4 hub genes associated with overall survival time. Several cells (Monocyte, Macrophage and Gamma delta cells) have been found to be significantly decreased in immune-high BRAF mutated SKCM group. While CD4⁺T, CD8⁺T, CD4 naïve, Tr1, Th2 and many T cell subsets were significantly increased in immune-high group. These immune cells and genes were closely related to each other. This study revealed that the dysregulation of immune function and immune cells may contribute to the poor outcomes of BRAF mutated patients. It is of great significance to our further understanding of the TME and immune dysfunction in BRAF mutated SKCM.

Introduction

Skin cutaneous melanoma (SKCM) is one of the most aggressive malignancies, causing about 80% of deaths in skin cancer (Schadendorf et al., 2018). Nearly 50% of cutaneous melanoma harbor activating V600E mutations in BRAF, which is considered a prognostic indicator of tumor proliferation, metastasis, recurrence as well as an effective target for SKCM treatment. The major factor limiting the clinical benefit of BRAF inhibitor are short response duration, off-target effect and drug resistance (Ribas et al., 2019a). There is also evidence that oncogenic BRAF can modulate the ability of the immune system to recognize SKCM cells. Activating mutations in the BRAF gene activate the mitogen-activated protein kinase (MAPK) pathway, which contributes to immune escape by recruiting regulatory T cells, reducing antigen presentation, and inhibiting the release of IFN-γ and TNF-α (Ascierto and Dummer, 2018). A series of immunotherapy strategies such as anti-PD-1, anti-CTLA4 and MAGE-A3 have been applied in SKCM and result in improvement in patient survival (Bajor et al., 2018). In addition, the combination of BRAF inhibitors and anti-PD-1 has shown significant improvement in SKCM treatment response (Ascierto et al., 2019a). These results suggest that we may be able to improve the survival outcome of BRAF mutated patients by regulating their immune response and tumor microenvironment. The tumor microenvironment (TME) consists of a variety of immune cells and stromal cells, including fibroblasts, endothelial cells, extracellular matrix, cytokines, chemokines and receptors (Vigneron, 2015). Due to the complexity of TME and a lack of a rational mechanistic basis, it is urgent to investigate the tumor microenvironment and identify prognostic biomarkers in BRAF mutated SKCM patients (Gnanendran et al., 2020).

ESTIMATE algorithms have been developed to calculate tumor purity in various cancers based on the specific gene expression signature of immune and stromal cells (Yoshihara et al., 2013). In this current work, we applied the expression data of BRAF mutated SKCM cohorts and ESTIMATE algorithm to extract a list of tumor microenvironment associated genes. Most of the genes were found to be related to better survival outcomes in BRAF mutated SKCM patients. Importantly, we estimated the proportion of immune cells based on gene expression profiling in BRAF mutated samples. Finally, we identified 4 hub genes associated with prognosis and immune cell infiltration in BRAF patients.

Materials and methods

Database of BRAF mutated SKCM patients

Transcriptional data of BRAF mutated SKCM patients (n = 240)was downloaded from the TCGA (https://tcga-data.nci.nih.gov/tcga/). In addition, their age, sex, tumor stage and survival information were obtained from the clinical documents in TCGA database (Tomczak et al., 2015). Statistical information of BRAF mutated SKCM patients was downloaded from Tumor Immune Estimation Resource dataset (Li et al., 2017). As a validation dataset, transcriptional data of SKCM patients (n = 131) was download from Gene Expression Omnibus (GEO) (Barrett et al., 2013). Screening criteria include: 1) the clinical diagnosis was skin cutaneous melanoma and 2) detection of BRAF mutation character. The discharge criteria include: 1) clinical data without survival time and outcome, and 2) datasets with small sample sizes (n < 50). Finally, the datasets were eligible: accession number GSE22153 (n = 131).

Calculation of immune and stromal scores

We used the ESTIMATE method to calculate the immune score and stromal score for each patient (Yoshihara et al., 2013). It is widely used to characterize the composition of infiltrating stromal cells and immune cells in tumor tissues.

Analysis of DEGs

BRAF mutated SKCM patients were ranked and divided into top and bottom halves (high vs. low score groups) based on their immune scores. Similarly, based on the stromal scores, the SKCM samples were grouped into high-stromal group and low-stromal group. Differentially expressed genes (DEGs) between the high-immunity/high-stromal group and low-immunity/low-stromal group were identified using the “limma” package in R software. |log(Fold change)| > 2, p < 0.05 and FDR<0.05 were set as the cutoffs.

Survival analysis

Overall survival data collected from each BRAF mutated SKCM patient were used to perform Kaplan-Meier analysis to explore the prognostic genes among the above DEGs. Patients with a given gene expression above 50% were designated as the high-expression group, while those with gene expression below 50% were designated as the low-expression group. Using log-rank method to test significance. The p value <0.01 was set as the cut-off value. Then, based on the survival data from GSE22153, we verified the prognostic value of prognostic genes in TCGA. The validated prognostic genes were used for subsequent protein-protein interaction analysis.

Function annotation

In order to reveal the function of DEGs and module genes, function annotation and Genome (KEGG) pathway enrichment analysis were performed using DAVID (Huang et al., 2007). FDR< 0.05 and p < 0.01 were set as the cut-off.

Protein-protein interaction network and model analysis

Evaluation of the protein-protein interaction network coded by validated prognostic genes was constructed by STRING (Szklarczyk et al., 2015), and their co-expression network was displayed by Cytoscape (Shannon et al., 2003). Then, the plugin Molecular Complex Detection (MCODE) was applied to identify the module genes that interact most closely.

Hub genes selection, validation and their co-expression network

Hub genes were obtained in this study by using Cytohubba plugin. The top ten genes in our PPI network were calculated based on six algorithms (MCC, MNC, EPC, Closeness, Radiality, Degree) at the same time. In addition, the intersection genes contained in the results of the six algorithms are screened out by upset calculation. Then, we used these genes as hub genes for further analysis. By using the Genemania database, we constructed the co-expression network of these hub genes and investigated their function (Warde-Farley et al., 2010). We used GSE22153 data to verify the mRNA expression of hub genes. To further validate our findings, we searched the Human Protein Atlas (https://www.proteinatlas.org/) website for the immunohistochemical (IHC) staining results of nine hub genes in normal skin and tumor tissue.

Immune cell components of BRAF mutated SKCM patients

To quantify the immune cell components of BRAF mutated SKCM patients, the expression data of patients were applied to calculate the composition of infiltrating immune cells by using ImmuCellAI algorithm (Racle et al., 2017). Based on the transcribed data of tumor tissue, the deconvolution algorithm can well reflect the infiltration and composition of immune cells. In this article, 24 kinds of immune cells such as neutrophils and NKT were calculated using ImmuCellAI algorithm. In addition, we compared the difference of immune cells between immune-high group and immune-low group using t-test. Moreover, the spearman correlation coefficient was calculated between immune cells and hub genes.

Statistical analysis

Analysis of DEGs, function annotation, survival analysis, ROC curves were all performed and visualized in R software. t-test was used to calculate the significant difference of immune and stromal scores among different AJCC stages and Breslow depth. p-values<0.05 were considered a statistically significant cut-off in all tests.

Results

Clinical information of BRAF mutated cutaneous melanoma patients

According to the inclusion criteria, 240 BRAF mutated SKCM patients from TCGA and 131 SKCM patients from GSE22153 (n = 131) were collected finally. In our study, the clinicopathological characteristics of BRAF mutated SKCM patients were shown in Table 1.

TABLE 1

TABLE 1. Clinicopathological characteristics of SKCM patients.

Immune and stromal score are closely related to the prognosis of SKCM

According to the ESTIMATE results, the immune score of 240 BRAF mutated SKCM patients (TCGA) ranged from -1133.65 to 3441.88. In addition, stromal score of BRAF mutated SKCM patients ranged from -1597.24 to 1817.91. We evaluated the correlation between immune, stromal score and clinicopathological characteristics of SKCM patients. In Figures 1C,D, we found that when the Breslow depth >3 mm, the immune score was significantly lower than that in 0–1.5 mm group and 1.5–3 mm group (p < 0.05). Similar results could be found in the stromal score. We also found that there was a significant correlation between immune score, stromal score and AJCC stage (American Joint Committee on Cancer) of SKCM (p < 0.05, Figures 1A,B). When comparing the immune and stromal scores of different AJCC stages, significant differences could be observed between several groups (I vs. II, II vs. III and II vs. IV).

FIGURE 1

FIGURE 1. Immune scores and stromal scores are closely associated with BRAF mutated melanoma prognosis. (A,B) the correlation between immune/stromal score and AJCC stage. (C,D) the correlation between immune/stromal score and Breslow depth. (E,F) high immune scores and stromal scores were associated with longer survival (p < 0.05). *p < 0.05; **p < 0.01; ***p < 0.001.

We further analyzed the relationship between immune and stromal scores and the prognosis of SKCM. A total of 240 BRAF mutated SKCM patients were ranked according to their immune scores and stromal scores. Then, we divided the 240 SKCM cases into top (n = 120) and bottom halves (n = 120) based on their scores. Among them, high level of immune score and stromal score were found significantly associated with longer overall survival (Figures 1E,F, p < 0.05).

Differentially expressed genes between high vs. low group and their function annotation

In view of the fact that immune and stromal scores were closely related to SKCM prognosis. Differentially expressed genes between high vs. low group were identified. The heatmap of gene expression showed a significant difference between immune high and immune low group. Similar results could be found between stromal high and stromal low group (Figures 2A,B). As a result, there were 1310 genes upregulated and 47 genes downregulated between high immune group and the low group (|logFoldChange| >2; p < 0.05). Additionally, there were 1478 genes upregulated and 39 genes downregulated between high stromal group and low group (|logFoldChange| >2; p < 0.05) (Figures 2E,F).

FIGURE 2

FIGURE 2. (A,B) heatmaps of gene expression profiles of samples between high immune/stromal and low immune/stromal groups. (C,D) the up-regulated and down-regulated overlapped DEGs. (E,F) volcano plot of immune and stromal DEGs.

Through the intersection of the Venn diagram, there were 990 overlap genes which both upregulated in the immune and stromal groups (Figure 2C). There were only 5 genes which both downregulated in the immune and stromal groups (Figure 2D). Therefore, the overlapped 990 genes were selected for further analysis. Function annotation has been carried out among the 990 overlap genes (p < 0.01; FDR<0.01). BP category suggested that immune response, defense response, inflammatory response, positive regulation of immune system and leukocyte activation were important process of the overlap genes (Figure 3A). As expect, MF results indicated that these upregulated overlap genes were mostly involved in sugar binding, cytokine activity, chemokine receptor binding and chemokine activity (Figure 3B). The plasma membrane, intrinsic to the plasma membrane, plasma membrane part items in CC category, indicating 990 overlapped genes play their roles in the plasma membrane (Figure 3C). In addition, chemokine signaling pathway, cytokine-cytokine receptor interaction, and cell adhesion molecules were important pathways of the overlapped gene network (Figure 3D).

FIGURE 3

FIGURE 3. Top 10 GO terms (BP, MF, CC) and KEGG analysis of overlap DEGs (p < 0.01). (A) BP results. (B) MF results. (C) CC results. (D) KEGG pathways. GO, Gene Ontology; BP, biological process; MF, molecular function; CC, cellular component; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Correlation of expression of individual DEGs in overall survival

The overlap 990 upregulated genes were used to identify prognostic genes through survival analysis. Subsequently, 755 genes (76%) were found correlated with longer overall survival time (p < 0.01, Figure 4, Supplementary Table S1). These genes were considered immune-related prognostic genes for further study.

FIGURE 4

FIGURE 4. Survival curves for immune-related genes in TCGA cohort.

Survival verification in GEO cohort

We collected BRAF mutated SKCM patients from GSE22153 from GEO database. Based on their overall survival data, 755 prognostic genes were selected to further verify their survival value. As a result, a total of 107 genes out of 755 identified genes were validated (Figure 5) to be significantly linked to longer overall survival time (Supplementary Table S2). We insisted that these 107 genes were potential prognostic immune-related biomarkers for BRAF mutated SKCM patients.

FIGURE 5

FIGURE 5. Survival curves for immune-related genes in GEO cohort.

Protein-protein network among genes of prognostic value

The protein-protein interaction (PPI) networks were constructed among 107 immune-related prognostic genes to explore their potential interactions and find the co-expression network (Figure 6A). 95 nodes and 813 edges were obtained in 107 gene interactions. Moreover, we used MCODE Plug-in to select the gene modules that interact most closely in the PPI network (module nodes> 6). In module 1 (Figure 6B), 241 edges involving 25 nodes were formed in the network. TYROBP, CD86, CSF1R, ITGB2 were found most closely related to other genes. In module 2 (12 nodes and 23 edges), LAPTM5 and VSIG4 had the higher connection values, indicating their core role in the module (Figure 6C). In module 3 (7 nodes and 13 edges), several HLA-related genes such as HLA-DQA1 and HLA-DPB1 had the higher connection values (Figure 6D).

FIGURE 6

FIGURE 6. (A) PPI networks of 107 prognostic genes. (B–D) gene model 1, gene model 2 and gene model 3.

Selection of hub genes and their co-expression network

Based on the above PPI network, we evaluated the top 10 genes of BRAF patients using six algorithms (Table 2). Eventually, TYROBP, CD86, CSF1R and ITGB2 were present in six algorithms at the same time (Figure 7B). By constructing the co-expression network of hub genes, the genetic interactions and pathways were analyzed (Figure 7A). Function annotation revealed that these hub genes and their co-expression genes were mainly related to leukocyte activation and lymphocyte proliferation (Figures 7C,D).

TABLE 2

TABLE 2. Screening of hub genes using six algorithms in cytoHubba. The bold value represents Hub genes

FIGURE 7

FIGURE 7. Selection and co-expression network of hub genes. (A) co-expression network of hub genes. (B) screening hub genes based on six algorithms. (C,D) function analysis and networks of hub genes and their co-expression genes.

Immune infiltration results between high vs. low group and their association with hub genes

Based on the expression data and ImmuCellAI algorithm, we quantified the immune cell components of BRAF mutated SKCM patients (Figures 8A,C). As shown in the figure, several cells (Monocyte, Macrophage and Gamma delta cells) have been found to be significantly decreased in immune-high group. While CD4⁺T, CD8⁺T, CD4 naïve, Tr1, Th2 and many T cell subsets were significantly increased in immune-high group (Figure 8B). Similarly, when compared the stromal-high and stromal-low group, Macrophage and Gamma delta cells were found to be significantly decreased in high group. While CD4⁺T, CD8⁺T, CD4 naïve, Tr1, Th2 and many T cell subsets were significantly increased in stromal-high patients. (Figure 8D).

FIGURE 8

FIGURE 8. The immune landscape of BRAF mutated samples microenvironment. (A) The landscape of immune cells among immune high and low group. (B) Immune cell differences between immune high and low group. (C) The landscape of immune cells among stromal high and low group. (D) Immune cell differences between stromal high and low group. *p < 0.05; **p < 0.01; ***p < 0.001.

Additionally, our results revealed that the expression of these hub genes may be related to the imbalance of immune cells (Figure 9B). For example, 4 hub genes (TYROBP, CD86, CSF1R and ITGB2) were mainly positively related to CD4⁺T, CD8⁺T, CD4 naïve, Tr1, iTreg, Tfh and many T cell in SKCM. While these hub genes could be found to be significantly negatively related to Macrophage and Gamma delta cells.

FIGURE 9

FIGURE 9. Validation of hub genes expression. (A) expression levels of hub genes in GSE22153. (B) association between hub genes and immune cells. (C) expression levels of hub genes in GSE22153. (D,E), ROC curves of hub genes. (F) IHC results of hub genes.

Validation of hub genes expression and their ROC curves

In order to verify our results, transcriptional data of GSE22153 was used to analysis the expression of these hub genes. Our results showed that all of the hub gene expression results were consistent with the previous description (Figures 9A,C). In addition, these hub genes have good efficacy in the diagnosis of immune-high and low group (AUC > 0.82, Figure 9D). Similar results could be found between stromal high and stromal low group (Figure 9E). The IHC results indicated that these four hub genes were significantly differentially expressed between normal and tumor tissues (Figure 9F).

Discussion

Skin cutaneous melanoma (SKCM) is one of the dead cancers with high malignant metastasis and mortality rates (Siegel et al., 2020). Identification of oncogenes provides novel insights into the progression of cancer therapy. BRAF oncogene was found in more than 50% of skin cutaneous melanoma as well as other cancers such as colorectal cancer and papillary thyroid cancer (Pollock and Meltzer, 2002; Rajagopalan et al., 2002; Cabanillas et al., 2020). The majority of researches claimed that BRAF mutation was often associated with high risk of metastasis, recurrence and poor survival outcomes (Ascierto et al., 2019b). BRAF inhibitor was considered the foundation of BRAF mutated melanoma treatment, and have demonstrated success and enhanced patient survival. However, only about 33% of patients benefit from target therapy in 5-year overall survival and the major limitations include short response duration, development of drug tolerance, and off-target effects (Robert et al., 2019). Recently, a growing group of researches showed that oncogenic BRAF can decrease the ability of the immune system to recognize melanoma cells. And the inhibition of BRAF can restore tumor immune recognition (Boni et al., 2010). Previous studies based on mice demonstrated that BRAF inhibitor response durations in vivo were significantly longer when melanoma cell lines were grown in immunocompetent mice compared to immunocompromised (Smalley, 2020). Besides, BRAF inhibition was associated with increased infiltration of CD4⁺T, CD8⁺T cells and reduced levels of myeloid-derived suppressor cells (MDSCs), tumor-associated macrophages (TAMs). While the depletion of CD4⁺ and CD8⁺ T cells significantly blunted the BRAF inhibitors response (Ribas et al., 2019b). Over the past years, the combination of immune therapy (anti-PD-1, anti-CTLA-4) and target therapy (BRAF inhibitor) has achieved an impressive improvement of the patients’ survival (Erkes et al., 2020). All results above suggested that the tumor microenvironment and immune effects play a vital role in SKCM therapy. However, the exploration of BRAF mutated immune microenvironment and the identification of immune-related prognostic targets in SKCM patients are still lacking (Boussadia et al., 2018).

Tumor immune microenvironment (TME) is described as significantly affecting the cancer treatment and prognosis (Hinshaw and Shevde, 2019). The ESTIMATE has been applied in glioma, renal cell carcinoma and gastrointestinal tumors, showing the validity of this algorithm in estimating tumor purity (Alonso et al., 2017; Jia et al., 2018). Therefore, ESTIMATE algorithm was applied to identify immune-related prognostic genes that contributed to patients’ overall survival by investigating the TME. Our research showed that both of the immune and stromal scores were inversely correlated with Breslow depth and AJCC stage, which have been considered as classical prognostic factors for SKCM. As shown in the Kaplan-Meier survival curve that patients with a higher immune score had longer overall survival time than those with a lower immune score in the BRAF mutated SKCM.

Through the DEGs analysis, we found that there were 990 overlapped genes which both upregulated in the immune and stromal groups. Function annotation indicated that immune response, defense response, positive regulation of immune system process and cytokine binding were important biological processes of the 990 overlapped genes. Pathway analysis demonstrated that the majority of the overlapped genes served a role in chemokine signaling pathways, cytokine-cytokine receptor interaction and cell adhesion molecules. As expected, dysregulation of immune function had a significant impact on the microenvironment of BRAF mutant SKCM patients. Among them, chemokine signaling is mainly involved in the recruitment of various immune cells, and their dysregulation may be an important reason for reducing the level of immune infiltration and leading to poor prognosis in patients with BRAF mutation. Based on our results and previous researches, it is conceivable to hypothesize that chemokines and immune response play a vital role in the regulation of SKCM TME (Huang et al., 2020; Li et al., 2020).

Subsequently, 755 genes (76%) were found correlated with longer overall survival time. This further demonstrated the clinical value of these immune microenvironment-related genes in patients with BRAF mutations. Subsequently, 107 prognostic genes of BRAF patients were verified in GEO data set. PPI network and model analysis had identified 4 hub genes (TYROBP, CD86, CSF1R and ITGB2) in our study. Moreover, our hub genes occupied a central position in both the PPI network and Model 1, proving their core position and clinical value. They were mainly associated with oncogenic transformation, immune response and regulating of immune cells (D'Angelo et al., 2019; Casey et al., 2016). For example, related studies have shown that the activation of T cells requires costimulatory signals produced by the interaction of CD28 and CD86, which could increase the infiltration of T cells in tumor tissue and prolong the survival time of mice (Jia et al., 2022). The combination of anti-PD1 and anti-CSF1 receptor (CSF1R) antibodies induced the regression of melanoma in-driven transplanted mice (Neubert et al., 2018). ITGB2 is associated with immune infiltration of multiple immune cell subsets, such as CD45, CD8, CD4T cells, CD20B cells and so on (Kwak et al., 2021). Although there was no direct evidence for the association between TYROBP and melanoma, given the important association between these hub genes and immune infiltration, we regard them as potential therapeutic targets for patients with BRAF mutations.

Previous studied demonstrated that the imbalance of immune cell components was closely related to progressive disease and poor prognosis (Qiao et al., 2019). Therefore, we conducted a further immune infiltration analysis. As expected, CD4⁺T, CD8⁺T, CD4 naïve, Tr1, Th2 and many T cell subsets were significantly increased in immune-high group. While several cells (Monocyte, Macrophage and Gamma delta cells) have been found to be significantly decreased in immune-high BRAF mutated SKCM group. There were significant differences in immune-infiltration between the two groups, which may help to identify groups that are more responsive to BRAF inhibitors. Monocyte-lymphocyte ratio (MLR) is considered to be an important indicator of tumor prognosis (Garcia et al., 2022). It has been reported that cancer-associated Macrophage play a key role in tumor progression, angiogenesis, invasion and recruitment of immunosuppressive cells (Samain and Sanz-Moreno, 2020). Persistent immune-related gene expression and T-cell penetration were associated with clinical benefit in SKCM patients (Shoushtari et al., 2022). The infiltrating levels of various effector T cells, such as CD4⁺ and CD8⁺ T, were significantly higher in the immune-high group than in the control group. After binding to MHC class I antigens on tumor cells via T cell receptors, CD8⁺ T cells can produce granzymes and perforin to destroy cancer cells (Tsukumo and Yasutomo, 2018). It is well known that CD8⁺ T cells have an antitumor effect, and the increase of CD8⁺ T cells can significantly improve the prognosis of SKCM patients (Chen et al., 2018). It is important to note the emerging role of CD4⁺ T cells in antitumor immunity, and in particular, their functional versatility in the context of the tumor immune microenvironment. In actual tumor therapy, the single immune function of CD8⁺ cells is not enough to destroy tumor cells, as the immune checkpoint inhibitors (anti-PD-1, anti-CTLA-4) are only 30% effective (Gellrich et al., 2020). Recent studies have found that the best initiation and maturation of MHC-I-restricted CD8⁺T cells is CTL (cytotoxic T lymphocytes), which requires the response of CD4⁺T cells (Alspach et al., 2019). By secreting interferon and promoting the proliferation and lethality of CD8⁺T cells in TME, CD4⁺T cells play a vital role (Zhu et al., 2015; Borst et al., 2018). In preclinical studies, it was found that BRAF inhibition led to increased CD40L expression and IFN-γ release from CD4⁺T cells, and decreased levels of multiple cytokines including IL1, IL6, and IL10 (Ott et al., 2013). Therefore, we speculate that increasing the proportion of CD4⁺T cells to enhance the lethality of CD8⁺T cells in TME may be a potential strategy to improve the prognosis of BRAF mutated SKCM patients.

Moreover, the expression of these hub genes was related to the imbalance of multiple immune cells. For example, 4 hub genes (TYROBP, CD86, CSF1R and ITGB2) were mainly positively related to CD4⁺T, CD8⁺T, CD4 naïve, Tr1, iTreg, Tfh and many T cell in SKCM. While these hub genes could be found to be significantly negatively related to Macrophage and Gamma delta cells. These results were consistent with their previous association with longer overall survival. There are few study on the relationship between the hub genes and BRAF mutated SKCM treatment. Therefore, we had identified several immune-related prognostic biomarkers for BRAF mutated patients. Finally, we preliminarily validated the expression of hub genes in another dataset and evaluated their diagnostic value. Our results showed that all of the hub genes significantly up-regulated in immune-high group. These data provide reference for further development of treatment for patients with BRAF mutations.

We must acknowledge the limitations in this study. First, more patients should be collected in the future to expand the sample size, which is conducive to a deeper understanding of the mechanisms of BRAF mutated SKCM and immune dysfunction. Second, we have limited experimental data and further function validation is required to investigate the interaction between the prognostic genes and immune cells.

Conclusion

For the first time in this study, we try to explore the TME to better understand the potential prognostic immune-related targets and mechanisms in BRAF mutated SKCM patients. This study revealed that the dysregulation of immune function and immune cells may contribute to the poor outcomes of BRAF mutated patients. It is of great significance to our further understanding of the TME and immune dysfunction in BRAF mutated SKCM.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

BH and WS designed and carried out the bioinformatics analyses and drafted the manuscript; BH and WS did the statistical analyses and drew the figures; DY initialed study. All authors read and approved the final manuscript.

Funding

This work was supported by China National Nuclear Corporation Medical Department “Nuclear Medicine Technology Innovation” Project (ZHYLYB2021009, ZHYLTD2021002). National Natural Science Foundation of China (32071238, 82073477). The Fundamental Research Funds for the Central Universities and Young Talent Program of China National Nuclear Corporation (CNNC2021136) and Natural Science Foundation of Sichuan Province (2020YJ0194). Foundation of China BaoYuan (CBYI202104).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

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.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.1081418/full#supplementary-material

Abbreviations

TME, tumor microenvironment; SKCM, skin cutaneous melanoma; GEO, Gene Expression Omnibus; TCGA, The Cancer Genome Atlas; DEGs, differentially expressed genes; PPI, protein−protein interaction; NS, normal skin; AJCC, American Joint Committee on Cancer; FDR, false discovery rate; FC, fold changes.

References

Alonso, M. H., Aussó, S., Lopez-Doriga, A., Cordero, D., Guinó, E., Solé, X., et al. (2017). Comprehensive analysis of copy number aberrations in microsatellite stable colon cancer in view of stromal component. Br. J. Cancer 117 (3), 421–431. doi:10.1038/bjc.2017.208

PubMed Abstract | CrossRef Full Text | Google Scholar

Alspach, E., Lussier, D. M., Miceli, A. P., Kizhvatov, I., DuPage, M., Luoma, A. M., et al. (2019). MHC-II neoantigens shape tumour immunity and response to immunotherapy. Nature 574 (7780), 696–701. doi:10.1038/s41586-019-1671-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Ascierto, P. A., Capone, M., Grimaldi, A. M., Mallardo, D., Simeone, E., Madonna, G., et al. (2019). Proteomic test for anti-PD-1 checkpoint blockade treatment of metastatic melanoma with and without BRAF mutations. J. Immunother. Cancer 7 (1), 91. doi:10.1186/s40425-019-0569-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Ascierto, P. A., and Dummer, R. (2018). Immunological effects of BRAF+MEK inhibition. Oncoimmunology 7 (9), e1468955. doi:10.1080/2162402X.2018.1468955

PubMed Abstract | CrossRef Full Text | Google Scholar

Ascierto, P. A., Ferrucci, P. F., Fisher, R., Del Vecchio, M., Atkinson, V., Schmidt, H., et al. (2019). Dabrafenib, trametinib and pembrolizumab or placebo in BRAF-mutant melanoma. Nat. Med. 25 (6), 941–946. doi:10.1038/s41591-019-0448-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Bajor, D. L., Mick, R., Riese, M. J., Huang, A. C., Sullivan, B., Richman, L. P., et al. (2018). Long-term outcomes of a phase I study of agonist CD40 antibody and CTLA-4 blockade in patients with metastatic melanoma. Oncoimmunology 7 (10), e1468956. doi:10.1080/2162402X.2018.1468956

PubMed Abstract | CrossRef Full Text | Google Scholar

Barrett, T., Wilhite, S. E., Ledoux, P., Evangelista, C., Kim, I. F., Tomashevsky, M., et al. (2013). NCBI GEO: Archive for functional genomics data sets-update. Nucleic Acids Res. 41, D991–D995. doi:10.1093/nar/gks1193

PubMed Abstract | CrossRef Full Text | Google Scholar

Boni, A., Cogdill, A. P., Dang, P., Udayakumar, D., Njauw, C. N., Sloss, C. M., et al. (2010). Selective BRAFV600E inhibition enhances T-cell recognition of melanoma without affecting lymphocyte function. Cancer Res. 70 (13), 5213–5219. doi:10.1158/0008-5472.CAN-10-0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Borst, J., Ahrends, T., Bąbała, N., Melief, C. J. M., and Kastenmüller, W. (2018). CD4(+) T cell help in cancer immunology and immunotherapy. Nat. Rev. Immunol. 18 (10), 635–647. doi:10.1038/s41577-018-0044-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Boussadia, Z., Lamberti, J., Mattei, F., Pizzi, E., Puglisi, R., Zanetti, C., et al. (2018). Acidic microenvironment plays a key role in human melanoma progression through a sustained exosome mediated transfer of clinically relevant metastatic molecules. J. Exp. Clin. Cancer Res. 37 (1), 245. doi:10.1186/s13046-018-0915-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Cabanillas, M. E., Dadu, R., Iyer, P., Wanland, K. B., Busaidy, N. L., Ying, A., et al. (2020). Acquired secondary RAS mutation in BRAF(V600e)-mutated thyroid cancer patients treated with BRAF inhibitors. Thyroid 30 (9), 1288–1296. doi:10.1089/thy.2019.0514

PubMed Abstract | CrossRef Full Text | Google Scholar

Casey, S. C., Tong, L., Li, Y., Do, R., Walz, S., Fitzgerald, K. N., et al. (2016). MYC regulates the antitumor immune response through CD47 and PD-L1. Science 352 (6282), 227–231. doi:10.1126/science.aac9935

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, H. Y., Xu, L., Li, L. F., Liu, X. X., Gao, J. X., and Bai, Y. R. (2018). Inhibiting the CD8(+) T cell infiltration in the tumor microenvironment after radiotherapy is an important mechanism of radioresistance. Sci. Rep. 8 (1), 11934. doi:10.1038/s41598-018-30417-6

PubMed Abstract | CrossRef Full Text | Google Scholar

D'Angelo, A., Sobhani, N., Roviello, G., Bagby, S., Bonazza, D., Bottin, C., et al. (2019). Tumour infiltrating lymphocytes and immune-related genes as predictors of outcome in pancreatic adenocarcinoma. PLoS One 14 (8), e0219566. doi:10.1371/journal.pone.0219566

PubMed Abstract | CrossRef Full Text | Google Scholar

Erkes, D. A., Cai, W., Sanchez, I. M., Purwin, T. J., Rogers, C., Field, C. O., et al. (2020). Mutant BRAF and MEK inhibitors regulate the tumor immune microenvironment via pyroptosis. Cancer Discov. 10 (2), 254–269. doi:10.1158/2159-8290.CD-19-0672

PubMed Abstract | CrossRef Full Text | Google Scholar

Garcia, J. S., Nowosh, V., López, R. V. M., and Massoco, C. O. (2022). Association of systemic inflammatory and immune indices with survival in canine patients with oral melanoma, treated with experimental immunotherapy alone or experimental immunotherapy plus metronomic chemotherapy. Front. Vet. Sci. 9, 888411. doi:10.3389/fvets.2022.888411

PubMed Abstract | CrossRef Full Text | Google Scholar

Gellrich, F. F., Schmitz, M., Beissert, S., and Meier, F. (2020). Anti-PD-1 and novel combinations in the treatment of melanoma-an update. J. Clin. Med. 9 (1), E223. doi:10.3390/jcm9010223

PubMed Abstract | CrossRef Full Text | Google Scholar

Gnanendran, S. S., Turner, L. M., Miller, J. A., Hwang, S. J. E., and Miller, A. C. (2020). Cutaneous adverse events of anti-PD-1 therapy and BRAF inhibitors. Curr. Treat. Options Oncol. 21 (4), 29. doi:10.1007/s11864-020-0721-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Hinshaw, D. C., and Shevde, L. A. (2019). The tumor microenvironment innately modulates cancer progression. Cancer Res. 79 (18), 4557–4566. doi:10.1158/0008-5472.CAN-18-3962

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, B., Han, W., Sheng, Z. F., and Shen, G. L. (2020). Identification of immune-related biomarkers associated with tumorigenesis and prognosis in cutaneous melanoma patients. Cancer Cell. Int. 20, 195. doi:10.1186/s12935-020-01271-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, D. W., Sherman, B. T., Tan, Q., Collins, J. R., Alvord, W. G., Roayaei, J., et al. (2007). The DAVID gene functional classification tool: A novel biological module-centric algorithm to functionally analyze large gene lists. Genome Biol. 8 (9), R183. doi:10.1186/gb-2007-8-9-r183

PubMed Abstract | CrossRef Full Text | Google Scholar

Jia, D., Li, S., Li, D., Xue, H., Yang, D., and Liu, Y. (2018). Mining TCGA database for genes of prognostic value in glioblastoma microenvironment. Aging (Albany NY) 10 (4), 592–605. doi:10.18632/aging.101415

PubMed Abstract | CrossRef Full Text | Google Scholar

Jia, H., Guo, J., Liu, Z., Chen, P., Li, Y., Li, R., et al. (2022). High expression of CD28 enhanced the anti-cancer effect of siRNA-PD-1 through prompting the immune response of melanoma-bearing mice. Int. Immunopharmacol. 105, 108572. doi:10.1016/j.intimp.2022.108572

PubMed Abstract | CrossRef Full Text | Google Scholar

Kwak, M., Erdag, G., Leick, K. M., Bekiranov, S., Engelhard, V. H., and Slingluff, C. L. (2021). Associations of immune cell homing gene signatures and infiltrates of lymphocyte subsets in human melanomas: Discordance with CD163(+) myeloid cell infiltrates. J. Transl. Med. 19 (1), 371. doi:10.1186/s12967-021-03044-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, B. H., Garstka, M. A., and Li, Z. F. (2020). Chemokines and their receptors promoting the recruitment of myeloid-derived suppressor cells into the tumor. Mol. Immunol. 117, 201–215. doi:10.1016/j.molimm.2019.11.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, T., Fan, J., Wang, B., Traugh, N., Chen, Q., Liu, J. S., et al. (2017). Timer: A web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res. 77 (21), e108–e110. doi:10.1158/0008-5472.CAN-17-0307

PubMed Abstract | CrossRef Full Text | Google Scholar

Neubert, N. J., Schmittnaegel, M., Bordry, N., Nassiri, S., Wald, N., Martignier, C., et al. (2018). T cell-induced CSF1 promotes melanoma resistance to PD1 blockade. Sci. Transl. Med. 10 (436), eaan3311. doi:10.1126/scitranslmed.aan3311

PubMed Abstract | CrossRef Full Text | Google Scholar

Ott, P. A., Henry, T., Baranda, S. J., Frleta, D., Manches, O., Bogunovic, D., et al. (2013). Inhibition of both BRAF and MEK in BRAF(V600E) mutant melanoma restores compromised dendritic cell (DC) function while having differential direct effects on DC properties. Cancer Immunol. Immunother. 62 (4), 811–822. doi:10.1007/s00262-012-1389-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Pollock, P. M., and Meltzer, P. S. (2002). A genome-based strategy uncovers frequent BRAF mutations in melanoma. Cancer Cell. 2 (1), 5–7. doi:10.1016/s1535-6108(02)00089-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Qiao, F., Pan, P., Yan, J., Sun, J., Zong, Y., Wu, Z., et al. (2019). Role of tumor-derived extracellular vesicles in cancer progression and their clinical applications (Review). Int. J. Oncol. 54 (5), 1525–1533. doi:10.3892/ijo.2019.4745

PubMed Abstract | CrossRef Full Text | Google Scholar

Racle, J., de Jonge, K., Baumgaertner, P., Speiser, D. E., and Gfeller, D. (2017). Simultaneous enumeration of cancer and immune cell types from bulk tumor gene expression data. Elife 6, e26476. doi:10.7554/eLife.26476

PubMed Abstract | CrossRef Full Text | Google Scholar

Rajagopalan, H., Bardelli, A., Lengauer, C., Kinzler, K. W., Vogelstein, B., and Velculescu, V. E. (2002). Tumorigenesis: RAF/RAS oncogenes and mismatch-repair status. Nature 418 (6901), 934. doi:10.1038/418934a

PubMed Abstract | CrossRef Full Text | Google Scholar

Ribas, A., Lawrence, D., Atkinson, V., Agarwal, S., Miller, W. H., Carlino, M. S., et al. (2019). Combined BRAF and MEK inhibition with PD-1 blockade immunotherapy in BRAF-mutant melanoma. Nat. Med. 25 (6), 936–940. doi:10.1038/s41591-019-0476-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Ribas, A., Lawrence, D., Atkinson, V., Agarwal, S., Miller, W. H., Carlino, M. S., et al. (2019). Publisher Correction: Combined BRAF and MEK inhibition with PD-1 blockade immunotherapy in BRAF-mutant melanoma. Nat. Med. 25 (8), 1319. doi:10.1038/s41591-019-0535-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Robert, C., Grob, J. J., Stroyakovskiy, D., Karaszewska, B., Hauschild, A., Levchenko, E., et al. (2019). Five-year outcomes with dabrafenib plus trametinib in metastatic melanoma. N. Engl. J. Med. 381 (7), 626–636. doi:10.1056/NEJMoa1904059

PubMed Abstract | CrossRef Full Text | Google Scholar

Samain, R., and Sanz-Moreno, V. (2020). Cancer-associated fibroblasts: activin A adds another string to their bow. EMBO Mol. Med. 12 (4), e12102. doi:10.15252/emmm.202012102

PubMed Abstract | CrossRef Full Text | Google Scholar

Schadendorf, D., van Akkooi, A. C. J., Berking, C., Griewank, K. G., Gutzmer, R., Hauschild, A., et al. (2018). Lancet 392 (10151), 971–984. doi:10.1016/S0140-6736(18)31559-9

PubMed Abstract | CrossRef Full Text

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res. 13 (11), 2498–2504. doi:10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Shoushtari, A., Olszanski, A. J., Nyakas, M., Hornyak, T. J., Wolchok, J. D., Levitsky, V., et al. (2022). Pilot study of ONCOS-102 and pembrolizumab: Remodeling of the tumor micro-environment and clinical outcomes in anti-PD1-resistant advanced melanoma. Clin. Cancer Res., OF1–OF10. doi:10.1158/1078-0432.CCR-22-2046

PubMed Abstract | CrossRef Full Text | Google Scholar

Siegel, R. L., Miller, K. D., and Jemal, A. (2020). Cancer statistics. Ca. Cancer J. Clin. 70 (1), 7–30. doi:10.3322/caac.21590

PubMed Abstract | CrossRef Full Text | Google Scholar

Smalley, K. S. M. (2020). Two worlds collide: Unraveling the role of the immune system in BRAF-MEK inhibitor responses. Cancer Discov. 10 (2), 176–178. doi:10.1158/2159-8290.CD-19-1441

PubMed Abstract | CrossRef Full Text | Google Scholar

Szklarczyk, D., Franceschini, A., Wyder, S., Forslund, K., Heller, D., Huerta-Cepas, J., et al. (2015). STRING v10: Protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 43, D447–D452. doi:10.1093/nar/gku1003

PubMed Abstract | CrossRef Full Text | Google Scholar

Tomczak, K., Czerwińska, P., and Wiznerowicz, M. (2015). The cancer genome Atlas (TCGA): An immeasurable source of knowledge. Contemp. Oncol. 19, A68–A77. doi:10.5114/wo.2014.47136

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsukumo, S. I., and Yasutomo, K. (2018). Regulation of CD8(+) T cells and antitumor immunity by notch signaling. Front. Immunol. 9, 101. doi:10.3389/fimmu.2018.00101

PubMed Abstract | CrossRef Full Text | Google Scholar

Vigneron, N. (2015). Human tumor antigens and cancer immunotherapy. Biomed. Res. Int. 2015, 948501. doi:10.1155/2015/948501

PubMed Abstract | CrossRef Full Text | Google Scholar

Warde-Farley, D., Donaldson, S. L., Comes, O., Zuberi, K., Badrawi, R., Chao, P., et al. (2010). The GeneMANIA prediction server: Biological network integration for gene prioritization and predicting gene function. Nucleic Acids Res. 38, W214–W220. doi:10.1093/nar/gkq537

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4, 2612. doi:10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, Z., Cuss, S. M., Singh, V., Gurusamy, D., Shoe, J. L., Leighty, R., et al. (2015). CD4+ T cell help selectively enhances high-avidity tumor antigen-specific CD8+ T cells. J. Immunol. 195 (7), 3482–3489. doi:10.4049/jimmunol.1401571

PubMed Abstract | CrossRef Full Text | Google Scholar