Role of Single Nucleotide Polymorphism-Related Genes in Tumour Immune Cell Infiltration and Prognosis of Cutaneous Melanoma

Background Advances in cancer research have allowed for early diagnosis and improved treatment of cutaneous melanoma (CM). However, its invasiveness and recurrent metastasis, along with rising resistance to newer therapies, have lent urgency to the search for novel biomarkers and the underlying molecular mechanisms of CM. Methods Single nucleotide polymorphism- (SNP-) related genes were obtained from the sequencing data of 428 CM samples in The Cancer Genome Atlas. Functional enrichment of these genes was analysed in clusterProfiler. Additionally, a protein-protein interaction (PPI) network was constructed with the Search Tool for the Retrieval of Interacting Gene (STRING) database. Gene Expression Profiling Interactive Analysis (GEPIA) was used to identify the expression and prognostic value of mutated genes. Finally, the Tumour Immune Estimation Resource (TIMER) analysed the relationship between gene expression and immune cell infiltration. Results We constructed a PPI network from the top 60 SNP-related genes. Mutated genes were mainly involved in calcium and oxytocin signalling pathways, as well as circadian entrainment. In addition, three SNP-related genes, BRAF, FLG, and SORL1, were significantly associated with patient prognosis. BRAF and SORL1 were positively associated with infiltration abundance of B cells, CD8+ T cells, CD4+ T cells, neutrophils, and dendritic cells, whereas FLG expression was negatively associated. Furthermore, higher immune cell infiltration was positively correlated with good prognosis. Conclusions Our study provides vital bioinformatic data and a relevant theoretical basis to further explore the molecular pathogenesis of CM and improve patient prognosis.


Introduction
Melanoma affects over 300,000 people worldwide annually [1], while survival in cutaneous melanoma (CM) has greatly increased with advancements in treatment, including the development of targeted therapy and immunotherapy [2]. However, a subset of patients became resistant. Identifying oncogenic biomarkers can help guide treatment decisions and elucidate the determinants of responses to immune and targeted therapies.
Single nucleotide polymorphisms (SNPs) can arise through transformation, transversion, deletion, or insertion. Because SNPs frequently occur throughout the genome, they provide an opportunity to identify mutations associated with a wide range of diseases, including cancers [3]. Increasing evidence suggests that susceptibility to malignant CM is associated with SNPs [4]. GLI-1 polymorphisms in the hedgehog pathway are risk factors for melanoma susceptibility and can be used as prognostic biomarkers [5]. Furthermore, the IRF4 SNP (rs12203592) and the MTAP (rs869330) variants are both associated with melanomaspecific survival [6]. Given their likelihood of influencing prognosis, an analysis of relevant SNPs could help identify new prognostic biomarkers for patients with CM.
The tumour microenvironment influences gene expression in tumour tissues and consequently has a strong effect on clinical outcomes [7][8][9][10]. Furthermore, immune cells in the microenvironment can be used for the prognostic assessment of multiple tumours, such as glioblastoma, breast cancer, and melanoma [11][12][13]. Notably, CM has a highly activated tumour microenvironment, with most immune system components involved in the cancer's initiation and progression [14]. The degree of immune cell infiltration significantly affects melanoma prognosis [15].

BioMed Research International
In this study, we aimed to identify novel immune-related prognostic biomarkers for CM and explore their underlying molecular mechanisms. To investigate the biological significance of SNPs in CM prognosis, we obtained SNP-related genes from The Cancer Genome Atlas (TCGA) and performed bioinformatic analysis. The results provide a theoretical basis for researchers to develop personalised treatment methods for patients with CM.  Figure 3: KEGG analysis showed that mutated genes were mainly enriched in the calcium and oxytocin signalling pathway and circadian entrainment. The redder the color of the dots, the more significant the KEGG term.

Materials and Methods
3 BioMed Research International were downloaded. Mutation data were analysed with the R Bioconductor package, "Maftools" [17].

Functional
Analysis. The R package "clusterProfiler" was used for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses [18]. The cut-off for significant enrichment was P < 0:05. Mutated genes were screened for enrichment in molecular functions, biological processes, cellular components, and biological pathways.

Protein-Protein Interaction (PPI) Network Construction.
The Search Tool for the Retrieval of Interacting Gene (STRING) database predicts physical and functional interactions between proteins [19]. Interactions, nodes, and subnetworks of the top 60 mutated genes were analysed by Cytoscape [20], a program that imports STRING networks but also integrates data from associated databases.

Survival Analyses and Identification of Potential
Prognostic Biomarkers. Patient survival was assessed using the Kaplan-Meier plots of results from gene arrays, RNA sequencing, and next-generation sequencing. These data were also used to screen for potential prognostic biomarkers by assessing the effect of mutant gene expression on patient prognosis with the Gene Expression Profiling Interactive Analysis (http://gepia.cancer-pku.cn) [21].

TIMER and Multivariate Cox Regression. The Tumour
Immune Estimation Resource (TIMER) (http://timer .cistrome.org/) can use RNA sequencing data to detect associations between mutated genes and immune cell infiltration [22]. In this study, TIMER was used to perform a multivariate Cox regression on cells involved in immune cell infiltration, including CD4+ T cells, CD8+ T cells, B cells, macrophages, and neutrophils. Hazard ratios and 95% confidence intervals were also calculated. CIBER-SORT in R was then used to identify immune cell subtypes and determine the relationship between those cells and risk scores [23]. In this study, we assessed the relationship between BRAF, filaggrin (FLG), and sortilinrelated receptor 1 (SORL1) expressions and immune cell infiltration in patients with CM.   Figure 1).

Functional Enrichment Analysis.
The results of GO analysis on mutated genes revealed enrichment in "multicellular organismal signalling", "muscle system process", and "regulation of ion transmembrane transport" terms under the biological process category. In the cellular component category, terms "sarcomere", "contractile fibre part", and "myofibril" were enriched. Lastly, in the molecular function category, the terms "gated channel activity", "calmodulin binding", and "motor activity" were enriched (Figures 2(a) and 2(b)). Results from the KEGG analysis showed that mutated genes were enriched in the calcium and oxytocin signalling pathways, as well as in circadian entrainment pathways (Figures 3(a) and 3(b)).

PPI Network and Correlation Analysis of Mutated Genes.
We explored the correlation between mutated genes using the STRING database and then constructed a PPI network ( Figure 4). We found correlations between mutations and expression of the top 60 genes.

Survival Analysis of Mutated Genes and Screening of
Prognostic Biomarkers. We then explored the relationship between the mutated genes, their expression, and prognosis. Using P < 0:05 as the significance level, three mutant genes that significantly correlated with prognosis were identified: BRAF, FLG, and SORL1. The expression of BRAF and SORL1 in the tumour samples was increased, while the expression of FLG was decreased (Figures 5(a)-5(c)). Based on the data, patients were divided into high-and lowexpression groups according to the median expression values. The Kaplan-Meier plot results showed that increased expression of BRAF and SORL1 was associated with a better prognosis (P = 3:655e − 02 and P = 4:351e − 04, respectively). In contrast, increased FLG expression resulted in poor prognosis (P = 9:672e − 03) (Figures 6(a)-6(c)). Moreover, both in vitro and in vivo experiments should be conducted in future studies.

Discussion
One of the most aggressive skin cancers, CM accounts for approximately 90% of global deaths attributed to this malignancy [24], while new treatments and early diagnosis have decreased mortality. However, 15%-20% of melanoma tumours do not respond to targeted therapies [25]. Moreover, treatment with programmed cell death protein 1 or cytotoxic T lymphocyte-associated antigen 4 antibodies does not benefit 40%-60% of patients [26]. Therefore, novel drug targets and their underlying mechanisms should be discovered to improve the assessment of malignancy and prognosis.
Our study identified three core mutated genes in patients with CM: BRAF, FLG, and SORL1. These genes were significantly associated with the prognosis of patients with CM. We thus consider mutant SNPs in these genes to be potentially carcinogenic markers that should benefit the early diagnosis and the design of individualised targeted therapy for CM.
A member of the RAF kinase family, BRAF mutations are common in many cancer types, including melanomas (60% occurrence), thyroid cancers (60%), colorectal cancers (15%), and non-small-cell lung cancers (5%-8%) [27]. BRAF inhibitors interfere with the mitogen-activated protein kinase (MAPK) signalling pathway that regulates melanoma proliferation and survival [28]. The MAPK pathway is also involved in T cell receptor signalling. Thus, BRAF interference exerts anticancer effects on the tumour microenvironment. Specifically, it increases intratumoural T cell infiltration and activity, decreases immunosuppressive cytokine levels, enhances melanoma-differentiation antigen expression, and introduces tumour antigens by regulating HLA-1 levels [29]. Through promoting tumour recognition by the immune system, these processes enhance antitumour T cell responses.

BioMed Research International
FLG is a pivotal structural protein of the stratum corneum; it encodes natural moisturising factors that are important to skin barrier function [30]. Thus, FLG mutation carriers may develop malignant melanoma [31]. Moreover, FLG is involved in allograft rejection and tumour necrosis factor-alpha signalling which is the primary pathway involved in interferon gamma response [32].
Lastly, SORL1 is involved in DNA repair and oxidative phosphorylation. SORL1 coprecipitates with endosomal receptor HER2 in cancer cells and recycles it back to the plasma membrane, thus regulating its subcellular distribution [33]. SORL1 is also involved in retromer-related endosomal trafficking and is a risk gene for Alzheimer's disease [34].
Functional analyses were performed to further investigate the molecular mechanisms of these mutant genes in CM progression. The results of GO analysis indicated enrichment in multicellular organismal signalling, muscle system processes, ion transmembrane transport, and gated channel activity. Additionally, KEGG analysis showed enrichment in the calcium and oxytocin signalling pathways, as well as in circadian entrainment.
We also found that BRAF and SORL1 expressions were positively associated with immune cell (B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils, and dendritic cells) infiltration in patients with CM. Additionally, their expression was positively correlated with improved prognosis, whereas mutated FLG expression was negatively correlated. Thus, higher levels of immune cell infiltration appear to be associated with better survival outcomes. To some extent, these patterns are consistent with previous findings suggesting that highly immune-infiltrated lung cancer evaded immune attacks by inhibiting new antigen expression [35]. Furthermore, a strong body of evidence indicates that lymphocyte infiltration in the tumour microenvironment is correlated with immunotherapeutic benefits [36][37][38][39]. Taken together, the data imply that the prognostic value of BRAF, FLG, and SORL1 is associated with immune cell infiltration. An in-depth analysis of any unique properties of immune cell infiltration may aid in the development of cancer immunotherapy targets.
It is reported that an imbalance in the immune cell component ratio is highly correlated with poor prognosis and low survival in patients with cancer [40,41]. In our study, we found that the differential expression of immune cell infiltration levels of three mutant genes depends on gene mutation types.
Our study had some limitations. We did not perform in vivo experiments to confirm the potential relationship between mutant genes, prognosis, and immune cell infiltration. Future studies that consider gene-gene and geneenvironment interactions are warranted to better clarify the molecular mechanisms of CM.

Conclusions
In summary, we identify novel biomarkers by elucidating the role of SNP-related genes in tumour immune infiltration and the prognosis of patients with CM. Our study provides new perspectives for the identification of prognostic indicators and offers an opportunity to optimise CM treatment.

Data Availability
The analyzed datasets generated during the study are available from the corresponding author on reasonable request.