The Landscape of the Tumor Microenvironment in Skin Cutaneous Melanoma Reveals a Prognostic and Immunotherapeutically Relevant Gene Signature

The tumorigenesis of skin cutaneous melanoma (SKCM) remains unclear. The tumor microenvironment (TME) is well known to play a vital role in the onset and progression of SKCM. However, the dynamic mechanisms of immune regulation are insufficient. We conducted a comprehensive analysis of immune cell infiltration in the TME. Based on the differentially expressed genes (DEGs) in clusters grouped by immune infiltration status, a set of hub genes related to the clinical prognosis of SKCM and tumor immune infiltration was explored. Methods: We analyzed immune cell infiltration in two independent cohorts and assessed the relationship between the internal pattern of immune cell infiltration and SKCM characteristics, including clinicopathological features, potential biological pathways, and gene mutations. Genes related to the infiltration pattern of TME immune cells were determined. Furthermore, the unsupervised clustering method (k-means) was used to divide samples into three different categories according to TME, which were defined as TME cluster-A, -B, and -C. DEGs among three groups of samples were analyzed as signature genes. We further distinguished common DEGs between three groups of samples according to whether differences were significant and divided DEGs into the Signature gene-A group with significant differences and the Signature gene-B group with insignificant differences. The Signature gene-A gene set mainly had exon skipping in SKCM, while the Signature gene-B gene set had no obvious alternative splicing form. Subsequently, we analyzed genetic variations of the two signatures and constructed a competing endogenous RNA (ceRNA) regulatory network. LASSO Cox regression was used to determine the immune infiltration signature and risk score of SKCM. Finally, we obtained 13 hub genes and calculated the risk score based on the coefficient of each gene to explore the impact of the high- and low-risk scores on biologically related functions and prognosis of SKCM patients further. The correlation between the risk score and clinicopathological characteristics of SKCM patients indicated that a low-risk score was associated with TME cluster-A classification (p < 0.001) and metastatic SKCM (p < 0.001). Thirteen hub genes also showed different prognostic effects in pan-cancer. The results of univariate and multivariate Cox analyses revealed that risk score could be used as an independent risk factor for predicting the prognosis of SKCM patients. The nomogram that integrated clinicopathological characteristics and immune characteristics to predict survival probability was based on multivariate Cox regression. Finally, 13 hub genes that showed different prognostic effects in pan-cancers were obtained. According to immunohistochemistry staining results, Ube2L6, SRPX2, and IFIT2 were expressed at higher levels, while CLEC4E, END3, and KIR2DL4 were expressed at lower levels in 25 melanoma specimens. Conclusion: We performed a comprehensive assessment of the immune-associated TME. To elucidate the potential development of immune-genomic features in SKCM, we constructed an unprecedented set of immune characteristic genes (EDN3, CLEC4E, SRPX2, KIR2DL4, UBE2L6, and IFIT2) related to the immune landscape of TME. These genes are related to different prognoses and drug responses of SKCM. The immune gene signature constructed can be used as a robust prognostic biomarker of SKCM and a predictor of an immunotherapy effect.

The tumorigenesis of skin cutaneous melanoma (SKCM) remains unclear. The tumor microenvironment (TME) is well known to play a vital role in the onset and progression of SKCM. However, the dynamic mechanisms of immune regulation are insufficient. We conducted a comprehensive analysis of immune cell infiltration in the TME. Based on the differentially expressed genes (DEGs) in clusters grouped by immune infiltration status, a set of hub genes related to the clinical prognosis of SKCM and tumor immune infiltration was explored.
Methods: We analyzed immune cell infiltration in two independent cohorts and assessed the relationship between the internal pattern of immune cell infiltration and SKCM characteristics, including clinicopathological features, potential biological pathways, and gene mutations. Genes related to the infiltration pattern of TME immune cells were determined. Furthermore, the unsupervised clustering method (k-means) was used to divide samples into three different categories according to TME, which were defined as TME cluster-A, -B, and -C. DEGs among three groups of samples were analyzed as signature genes. We further distinguished common DEGs between three groups of samples according to whether differences were significant and divided DEGs into the Signature gene-A group with significant differences and the Signature gene-B group with insignificant differences. The Signature gene-A gene set mainly had exon skipping in SKCM, while the Signature gene-B gene set had no obvious alternative splicing form. Subsequently, we analyzed genetic variations of the two signatures and constructed a competing endogenous RNA (ceRNA) regulatory network. LASSO Cox regression was used to determine the immune infiltration signature and risk score of SKCM. Finally, we obtained 13 hub genes and calculated the risk score based on the coefficient of each gene to explore the impact of the high-and lowrisk scores on biologically related functions and prognosis of SKCM patients further.

INTRODUCTION
Skin cutaneous melanoma (SKCM) is currently one of the most lethal human malignancies. As the deadliest form of skin cancer, it accounts for almost 75% of skin cancer lethality (Siegel et al., 2021). Although the 5-year survival of early stage SKCM patients exceeds 95% (Thompson et al., 2009), the reported survival time for advanced-stage melanoma barely exceeds 1 year (Fecher et al., 2007). Currently, patients with primary melanoma require surgical resection as a first-line therapy. However, advanced melanoma is highly aggressive, making it insensitive to radiotherapy and chemotherapy (Goodson and Grossman, 2009). The emergence of immune checkpoint inhibitors, such as ipilimumab (Jameson-Lee and Luke, 2021) and nivolumab (Zhao et al., 2020), for melanoma has revolutionized the treatment of SKCM and offers new hope for patients. However, approximately 50% of patients do not benefit from immune checkpoint inhibitors (Hodi et al., 2010;Topalian et al., 2012).
Recently, growing evidence has demonstrated that tumor microenvironment (TME) plays an important role in SKCM progression (Hanahan and Weinberg, 2011;Quezada et al., 2011). TME influences tumorigenesis and metastasis through various biological processes. In contrast to this, TME heterogeneity is also an important cause of alterations in prognosis and sensitivity to immunotherapy in various cancers (Koikawa et al., 2021).
Notably, it may regulate the immune response through a variety of mechanisms, thereby affecting the inner metabolism process and immunosuppressive state. Among these, tumor-infiltrating immune cells (TIICs) exhibit tumor-promoting effects according to the tumor type. In most cancers, CD8+ T cells play crucial roles in TME, inhibiting the proliferation and invasion of malignant cells. T cell-mediated immune responses to melanoma antigens have been extensively documented (Marron et al., 2021). Furthermore, immunosuppression may act as an additional tumor burden that fosters tumor growth or immune escape in TME. Hence, a comprehensive understanding of TME is urgently needed to improve the efficacy of immunotherapies.
Considering the positive effect of immunotherapy in SKCM patients, understanding the molecular composition and function of TME is important to facilitate effective diagnosis, prognosis, mitigation, and immunotherapeutic responsiveness of SKCM patients.
We integrated The Cancer Genome Atlas (TCGA)-SKCMindependent cohort and validated the predictive model in four additional independent cohorts (GSE8401, GSE35640, GSE15605, and GSE46517) from Gene Expression Omnibus (GEO) datasets, to develop and verify a new set of personalized immune signature models. We also analyzed clinical and pathological characteristics of all existing SKCM patients, including somatic cell copy number variation (CNV), tumor mutation burden (TMB), and gene variable splicing. Several bioinformatics methods were employed to estimate the abundance of immune cell infiltration in SKCM patients, and the correlation between genomic characteristics of the immune landscape and pathological characteristics and prognosis of SKCM was determined. Therefore, in the present study, we sought to develop and evaluate the risk-based model of six novel immune-related genes (EDN3, CLEC4E, SRPX2, KIR2DL4, Ube2L6, and IFIT2). In addition, immunohistochemistry (IHC) staining and qRT-PCR were used to verify the difference in the expression level of the novel 6-gene signature in 25 frozen SKCM tissue samples and normal tissues. Notably, gene signatures that affect the immune landscape may be closely related to different prognosis and treatment responses of SKCM. Accordingly, we established an unprecedented set of immune signatures that can be used as robust prognostic biomarkers and predictors of the immunotherapy effects in SKCM patients.

Acquisition and Preprocessing of the SKCM Expression Datasets
Data from two publicly available datasets were included in the study. RNA-seq data were extracted for 471 patients, in addition to clinical features, from TCGA-SKCM cohort (Tomczak et al., 2015). Clinical data were downloaded from the University of California Santa Cruz Xena browser 1 . Somatic CNV data and TMB were downloaded from the Genomic Data Commons Data Portal (Grossman et al., 2016). Data on alternative splicing events were downloaded from TCGASpliceSeq database. "RCircos" (Zhang et al., 2013) package was used to generate a map for the genome-wide CNV analysis of SKCM patients using 23 pairs of chromosomes. Somatic mutation data were downloaded in mutation annotation format, and the maftools package (Mayakonda et al., 2018) was used for visualization. We obtained SKCM microarray data from two GeneChips (GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array and GPL96-57554 [HG-U133A] Affymetrix Human Genome U133A Array) in GSE84014 (Xu et al., 2008), GSE356405 (Dreno et al., 2018), GSE156056 (Raskin et al., 2013), and GSE46517 (He et al., 2020) from GEO 1 as validation sets. The limma package (Ritchie et al., 2015) in R was used for gene expression normalization, while the sva (Leek et al., 2012) package (3.20) was used to correct for plate batch effects.

TIICs Analysis and Clustering of Samples
CIBERORT (Newman et al., 2015), an algorithm that quantifies the proportion of TIICs with 547 signature genes to infer the representation in bulk tumor transcriptomes, was used to evaluate the cell components of TIICs. The proportions of 22 types of TIICs in the cohort from TCGA were estimated. Cluster analysis was performed using an unsupervised hierarchical clustering method, an algorithmic approach that groups individuals with similar observations based on the Euclidean distance method. Three immune microenvironment subtypes were defined using the ConsensusClusterPlus package (Wilkerson and Hayes, 2010). The procedure was repeated 1,000 times to stabilize the stratification. 1 http://www.ncbi.nlm.nih.gov/geo

Differential Gene Expression Pattern Clustering Analysis
To determine genes associated with TME cell infiltration patterns, we divided patients into three different pattern types according to TME; these types were defined as TME cluster-A, -B, and -C. We further analyzed common differentially expressed genes (DEGs) between the three groups of samples and, based on the significance of the differential expression changes of the specific genes, divided the samples into two groups of signature genes: A and B. Signature gene-A was an immune-related specific gene with a significant difference, and Signature gene-B was an insignificantly different part of an immune-related specific gene. The limma package was used to analyze DEGs between these three groups of SKCM patients, and the significant DEGs were defined as genes with an absolute log value of fold-change >1 and FDR <0.05. The overlapping DEGs among the three groups were further analyzed as specific genes using the VennDiagram (Venn diagram) package (Chen and Boutros, 2011). The k-means clustering algorithm, which is an unsupervised clustering method, was used to cluster these specific genes into meaningful groups in the GEO datasets according to the expression of specific genes in TME cluster-A, -B, and -C. Meanwhile, based on changes in gene expression, specific genes were divided into two groups: Signature gene-A and -B. Among them, Signature gene-A and Signature gene-B were obtained by cluster analysis of the three subtypes, where Signature gene-A was highly expressed in cluster A and relatively less expressed in clusters B and C.

RNA Sequence Expression Analysis in Gene Expression Profiling Interactive Analysis
Gene Expression Profiling Interactive Analysis (GEPIA) 2 is an online data processing webpage containing RNA sequence expression information of 9,736 tumors and 8,587 normal tissue samples (Tang et al., 2019). The differential expression of hub genes between tumor and normal tissues, correlations, and survival prediction was determined using the GEPIA database. Student's t-test was used to analyze the correlation between the expression and clinicopathological features. Statistical significance was set at p < 0.05.

Gene Set Enrichment Analysis
Gene Ontology (GO) analysis was performed to illustrate the unique biological significance of gene expression of signature genes. Gene functions were categorized into three series: cellular components, molecular functions, and biological processes (BPs). Crucial pathways were identified using Kyoto Encyclopedia of Genes and Genomes pathway analysis. GO annotations were visualized using the R package clusterProfiler (Yu et al., 2012). Gene Set Enrichment Analysis (GSEA; Subramanian et al., 2005) is a calculation method for identifying the potential biological mechanisms between two biological states. GSEA was conducted in the Molecular Signatures Database (MSigDB) 3 , which provided hallmark gene sets to predict BPs between normal and SKCM samples.

Evaluation of Patient Biological Characteristics
We further analyzed correlations between different groups and some biologically related processes. Gene sets for storing genes related to certain BPs, including immune checkpoints, antigen processing, CD8+ T cells, and epithelial-mesenchymal transition (EMT) markers, such as EMT1, EMT2, and EMT3, angiogenesis characteristics, pan-fibroblast TGF-β response characteristics, WNT characteristics, DNA damage repair, mismatch repair, nucleotide excision repair, DNA replication, and antigen processing and presentation, were collected.

Establishment of the Immune Characteristic Model and Clinical Prediction Model
We adopted a two-step method to establish a signature-based risk score. First, univariate Cox regression was used to analyze the impact of signature genes on the prognosis of SKCM patients according to the cutoff, with p < 0.05. The most commonly used model to establish clinical prognosis and survival-related models is the Cox proportional hazards model owing to its flexibility and relative robustness. The Cox grade-scoring test was used for repeated prediction and prognostic analyses. The Cox gradescoring test was robust to outliers. Accordingly, outliers were removed for the sensitivity analysis. The potential prognostic hub genes were validated to enhance risk profiling after surgery and to define new targets in the prediction models for SKCM patients. To prove the significance of risk score combined with clinicopathological characteristics for a personalized evaluation of patient prognosis, we first tested the expression correlation of hub genes in different tumors in TCGA database and the predictive ability of the risk score for the prognosis of patients with different tumors. Multiple regression was the main analysis method for building immune prediction models. Multiple regressions were performed using standard statistical methods. In addition, we fitted the multivariate Cox regression model using the least absolute contraction and selection operator (LASSO). LASSO regression is a regularized approach commonly used for high-dimensional predictor selection. A system of risk score was established using the LASSO Cox proportional hazards model to identify gene signatures for predicting the overall survival (OS) of SKCM. A predictive score was developed using the weighted sum of genes, with the coefficients of the LASSO regularization.
Factors related to OS and clinical pathological characteristics were evaluated using univariate and multivariate analyses with 3 http://software.broadinstitute.org/gsea/index.jsp Cox and logistic regression, respectively. Variables with a value of p < 0.05 in the multivariate analysis were included in the prognostic model. The performance and discriminative ability were assessed using Harrell's concordance index. Nomograms were constructed to predict the 3-year, 5-year, and 10year survival rates of SKCM patients based on predictive models with identified prognostic factors. Calibration is defined as a prediction from the nomogram compared with the observed outcomes.

IHC Validation
Tumor tissues were collected from 25 melanoma patients in the Han Chinese group from 2015 to 2020. Informed consent was obtained from patients, and the study was approved by the First People's Hospital of the Foshan Subject Review Board. Paraffin-embedded tissues were sectioned to be 4 µm thick for IHC analysis. Antigen retrieval was performed by incubating the samples in citrate buffer (pH 6.0) for 15 min at 100 • C in a microwave oven. After blocking with a mixture of methanol and 0.75% hydrogen peroxide, sections were incubated overnight with appropriate dilutions of primary antibodies (CLEC4E, 1:500, Sigma; END3, 1:500, Sangon Biotech; IFIT2, 1:500, Proteintech; SRPX2, 1:600, Proteintech; Ube2L6, 1:500, Abcam; KIR2DL4, 1:500, Abcam), followed by incubation with a secondary antibody conjugated with HRP (goat antirabbit, 1:500, Cell Signaling Technology). The sections were washed three times with phosphate buffer saline and incubated with AEC (ZSGB-BIO). All specimens were examined by the cross-product (H score) of the percentage of tumor cell staining at each of three staining intensities. The intensity of immunopositivity was scored as follows: none, 0; weak, 1; moderate, 2; and strong, 3. For example, a particular tumor may have 50% cell staining at an intensity of 1 and 50% of cell staining at an intensity of 3, for a combined H score of 200 [(50 × 1) + (50 × 3) = 200], which yields a range from 0 to 300. The final score was graded by the H score as follows: low, H score 0-100; moderate, H score 101-200; and high, H score 201-300.

Statistical Analysis
All statistical analyses were performed using R version 3.6.2. Differences in continuous variables between the two groups were estimated by independent Student's t-test, and the differences between non-normally distributed variables were analyzed using the Mann-Whitney U rank-sum test. If the normality test failed, Kruskal-Wallis One Way ANOVA on Ranks was performed. Pearson's χ 2 test or Fisher's exact test was used to compare qualitative variables. p-Values and hazard ratios were obtained from univariate Cox proportional hazards regression models using the R package, survival (Zeng et al., 2019). The log-rank test was used to evaluate the significance of the difference in survival time between the two groups. Receiver operating characteristic curve analysis was conducted using the pROC package (Robin et al., 2011) to evaluate the prognostic capabilities of different risk models and the time-dependent AUC values. Univariate and multivariate Cox analyses were used to determine independent prognostic factors. All statistical tests were two-sided, and statistical significance was set at p < 0.05.

Immune Infiltration Analysis Related to SKCM Patients
We obtained the gene expression data of SKCM patients from TCGA database( Table 1) and analyzed 22 different immune cell infiltrations in each sample using the CIBERSORT algorithm (Newman et al., 2015; Figure 1A). The TME cell network depicted interactions between tumor immune cells (Supplementary Figure 1). The comprehensive status of the cell lineage and its impact on the overall survival of patients with SKCM were also analyzed ( Figure 1B). The TME cell network demonstrated that cell clusters B, C, and D were positively correlated with each other, including CD8+ T cells, naïve B cells, plasma cells, CD4 naïve T cells, activated NK cells, monocytes, resting dendritic cells, regulatory T cells, monocytes, and CD4 memory activated T cells. Moreover, there was a significant positive correlation among immune cells in cluster A, such as between gamma delta T cells, resting NK cells, neutrophils, activated mast cells, M0 macrophages, eosinophils, and activated dendritic cells. Meanwhile, M0 macrophages, activated CD4 memory T cells, and monocytes showed a significant negative correlation. NK cell resting and NK cell activation showed the same significant negative correlation. An exploratory analysis was also performed to measure survival benefits and potential risks. Results showed that CD4+ T cells, CD8+ T cells, activated NK cells, regulatory T cells, dendritic cells, and M1 macrophages were associated with shortened OS, while resting NK cells, neutrophils, M0 macrophages, and dendritic cells were activated and associated with prolonged OS. To construct the best clusters and classification, we used the ConsensusClusterPlus package to assess the stability of the clustering structure and divided SKCM patients into TME cluster-A, -B, and -C. Unsupervised hierarchical clustering was used to analyze the normalized immune cell fractions. The heatmap showed correlation between the infiltration abundance of 22 types of immune cells and immune scores in the three groups ( Figure 1C). Results showed that patients with TME cluster-A had higher immune scores, whereas patients with TME cluster-B and -C mainly had tumor purity and stromal scores. In addition, the principal component analysis results showed that based on the expression data of SKCM patients, the three TME cluster groups could be significantly distinguished ( Figure 1D). Survival analysis related to the TME phenotype showed that TME cluster-C (N = 92) was associated with a better prognosis (log-rank test, p < 0.001) ( Figure 1E). We selected and analyzed the gene expression profile data and immune cell infiltration abundance of SKCM samples from the GEO database. The ConsensusClusterPlus package was used to evaluate the stability of the clustering (Figure 2A). Subsequently, the heatmap showed that patients with TME cluster-A had higher immune scores (Figure 2B), which was consistent with the results in TCGA database. This finding also aligned with that of previous studies that have reported melanoma being a highly immunedependent malignant tumor.
Construction of the TME-Related Signature Gene in SKCM Patients TME has a crucial impact on tumor epigenetics, metastasis, and immune escape. We used the limma package to analyze the 486 DEGs among the three groups of TME cluster in TCGA database to determine the potential biological characteristics of different TME phenotypes. Similarly, 330 DEGs were analyzed in the GEO datasets ( Figure 3A). Unsupervised clustering using these DEGs provided three distinct clusters: TME cluster-A, -B, and -C ( Figure 3B). Simultaneously, according to the differential expression of DEGs in different clusters, genes were further divided into Signature gene-A and -B. GO enrichment analysis showed that Signature gene-A and -B had different unique BPs. Signature gene-A was associated with the overexpression of immune-activated genes (Figure 3C), while Signature gene-B showed upregulation of genes related to stromal and transmembrane receptors (Figure 3D). There were significant differences in the infiltration of TME cells ( Figure 3E) and the enrichment of related biological pathways ( Figure 3F) between the three GeneClusters, which was consistent with our functional enrichment results.

CeRNA Regulatory Network Construction and Signature Gene Expression in SKAM Patients
To further analyze differences between Signature gene-A and -B gene sets and determine whether these differences have a profound impact on cancer genetics, we analyzed genetic levels of single nucleotide polymorphisms. CNV and SNP analyses showed that there were obvious mutations in both gene sets. Notably, somatic SPTA1 mutations were more frequent in TME cluster-A ( Figure 4A). TNN was the most frequently observed somatic mutation in TME cluster-B ( Figure 4B). Moreover, the examination of this frequency change in CNV revealed that CNV changes were common in both groups and most were concentrated on the amplification of copy numbers. We then elucidated positions of these CNV changes on the chromosome (Figures 4C,D). In the AS analysis, Signature gene-A mainly had exon skipping (ES) in SKCM patients (Figure 4E), while Signature gene-B had no obvious AS forms ( Figure 4F). In addition, Signature gene-A and -B gene sets may have a certain significance for the prognosis of SKCM. Further research is necessary for the interaction between immune-related molecules. Based on the gene expression of Signature gene-A and -B, we screened out lncRNAs that may be related to immunity using correlation coefficients >0.4 and p < 0.05. The starBase database 4 was used to obtain targeted differentially expressed miRNAs to construct the regulatory ceRNA network of the mRNA-miRNA-LncRNA interaction (Figure 4G).

Construction of the Immune-Related Prognostic Gene Signature
To predict the impact of immune characteristics on the prognosis of patients better, we constructed a new prognostic-related risk scoring system. The Signature gene-A and-B were incorporated into univariate Cox analyses, and 233 genes related to prognosis were obtained (p < 0.05). The LASSO Cox analyses were further used for dimensionality reduction and model construction, and finally, a total of 16 hub genes were included in the risk scoring model (Figure 5A). A risk score (RS) formula was established by including individual normalized gene expression values weighted by their LASSO Cox coefficients. The risk score for each SKCM patient was calculated according to the locking of the coefficients in each gene signature. The high-and lowrisk groups were divided according to the median risk score. Kaplan-Meier analysis showed that patients with high-risk scores had a relatively poor prognosis ( Figure 5B). The time-dependent receiver operating characteristic curve analysis also showed that the risk score had a good predictive ability for OS ( Figure 5C). The area under the curve (AUC) of 1-year, 3-year, and 5-year OS was 0.713, 0.694, and 0.734, respectively. The distribution of each patient's risk score, survival status, and gene expression map is shown in Figure 5D.

GSEA
We analyzed the impact of the high-and low-risk groups on the biologically relevant functions of SKCM patients. GSEA revealed that pathways related to metabolism and oxidative phosphorylation were mainly enriched in the high-risk group (Figures 6A,C), while pathways related to immune response, including cytokine signaling pathway, JAK-STAT signaling pathway, and natural killer cell-mediated cytotoxicity, were significantly enriched in low-risk patients (Figures 6B,D). At the same time, expression levels of TME cells (Figure 6E), as well as some other pathways, such as angiogenesis, mismatchrelated features, and stromal-related features (Figure 6F), were significantly different between the high-and low-risk groups of SKCM patients (p < 0.05).

Correlation Analysis of the Risk Score and Clinicopathological Characteristics
We assessed the correlation between the risk score and the clinicopathological characteristics of SKCM. The analysis results showed that the low-risk score was correlated with TME cluster-A classification (p < 0.001; Figure 7A) and metastatic SKCM (p < 0.001; Figure 7B). In addition, the low-risk score was The interaction between immune cells in TME in TCGA-SKCM. Cell cluster-A, red; Cell cluster-B, blue; Cell cluster-C, brown; Cell cluster-D, purple. The size of the circle represents the impact of each TME cell type on survival, and the log-rank test was used for analysis. The green part in the center of the circle indicates that the cell protects against overall survival, and the black indicates the risk to overall survival. The lines connecting TME cells represent cell interactions. The thickness of the line represents the correlation strength estimated by Spearman correlation analysis, in which positive correlation is represented by red and negative correlation is represented by blue. (C) Heatmap showing the infiltration of TIICs in 471 SKCM patients in TCGA database combined with the immune score. Unsupervised clustering grouping the samples into three major clusters. (D) The PCA of the gene expression profile distinguishes patients in the TME cluster-A, -B, and -C groups in TCGA-SKCM. The results can distinguish three different immune infiltration pattern samples (TME cluster-A: blue, TME cluster-B: yellow, TME cluster-C: green). (E) The Kaplan-Meier curve of the patient's overall survival (OS) shows that TME infiltration is significantly related to the overall survival (Log-rank test, p = 0.014). correlated with gender (p = 0.031; Figure 7E), and significantly correlated with pathological stage (p = 0.0027; Figure 7F) and T stage (p < 0.001; Figure 7G). There was no significant correlation between the risk score and TMB, age, M stage, and N stage (p > 0.05; Figures 7C,D,H,I and Table 2).

Evaluation and Validation of the Prognostic Signature
Correlation analysis revealed that the expression of 13 hub genes in different tumors was significantly correlated (Figure 8A). The heterogeneity between tumors caused the risk score to have different prognostic effects on different cancers (Figure 8B). Univariate and multivariate Cox analyses showed that the risk score was an independent risk factor for predicting the prognosis of patients with SKCM ( Figure 8C). We included the index of p < 0.05 in the multivariate Cox model to construct a nomogram to predict the 1-, 2-, and 3-year survival probability of SKCM patients. The C-indexes [0.732 (95% CI: 0.697-0.767)] of the nomogram were used to calculate the discriminative ability of the nomogram, showing a high degree of discrimination ( Figure 8D). The calibration also showed a great agreement between the 1-year, 2-year, and 3-year OS estimates and the actual observed values of SKCM patients through a comparison to the nomogram ( Figure 8E).

Prognostic Value of Hub Genes in SKCM Patients
We focused on the potential prognostic value of EDN3, CLEC4E, SRPX2, KIR2DL4, UBE2L6, and IFIT2 as immune scores for melanoma patients. In the group of novel hub genes, low expression levels of UBE2L6, KIR2DL4, IFIT2, and CLEC4E were observed in tumor cells, while high expression levels of SRPX2 and EDN3 were found in SKCM. However, only UBE2L6 and IFIT2 showed significant expression differences in TCGA-SKCM cohort. A larger sample size may be required to verify the differential expression of these molecules in SKCM. The open online tool, GEPIA, was used to analyze the prognostic value of these novel hub genes. We found that the low expressions of UBE2L6, KIR2DL4, IFIT2, and CLEC4E were significantly associated with poor prognosis in SKCM.

Preliminary IHC Validation in Melanoma Specimens
We evaluated the expression of six novel genes in melanoma tissues. The IHC staining results showed that Ube2L6, SRPX2, and IFIT2 were expressed at higher levels, while CLEC4E, END3, and KIR2DL4 were expressed at lower levels in 25 melanoma specimens (Figure 9).

DISCUSSION
The incidence of SKCM continues to increase annually. The recurrence and metastasis of SKCM caused by functional variations in TME markedly contribute to the extremely poor prognosis of SKCM. TME plays a key role in the response to immunotherapy. Several treatments have improved the survival rate of patients with advanced diseases, such as radiotherapy, chemotherapy, and immunotherapy (Wagner et al., 2019). In particular, the application of immunotherapy in SKCM has enabled remarkable breakthroughs (Mitchell et al., 2018). However, there are still many patients who exhibit resistance to cancer immunotherapy and adoptive cell therapy due to the immunosuppressive barriers that exist in TME . Therefore, there is an urgent need to determine more therapeutic targets and prognostic biomarkers based on TME. In this study, we conducted a comprehensive assessment of immune cell infiltration in TME, identified novel tumor immune subtypes, and assessed the prognostic value of immune cells in SKCM patients. The prognostic characteristics proposed herein are reliable for predicting the survival of patients with SKCM.
To verify the key role of immune cell-infiltrated TME in SKCM, we performed independent microarray data analysis on four datasets in the GEO database. Notably, cross-platform research increased the reliability of our novel prognosis model. We used the CIBERSORT algorithm to evaluate the immune cell infiltration status of the SKCM cohort transcriptome data in TGCA database comprehensively and obtained 22 different immune cell infiltration abundances in each patient. We depicted the tumor immune cell interaction network, including the overall situation of the cell lineage and its impact on the overall survival of SKCM patients, which was reviewed and verified in the SKCM patient cohort from the GEO database. To evaluate the tumor heterogeneity and interaction with TME that can guide better and earlier targeted treatments better, we hierarchically clustered tumor samples based on essential differences in immune cell infiltration patterns. Principal component analysis (PCA) based on the expression profile data of SKCM patients could distinguish the three TME cluster groups well. Thus, we proceeded to analyze the immune-related scores of three groups, tumor purity, and matrix scores. Patients with TME cluster-A had higher immune scores, whereas patients with TME cluster-B and -C mainly had tumor purity and matrix scores. In the GEO database, the GSE8401, GSE35640, GSE15605, and GSE46517 datasets verified the gene expression profile data and immune cell infiltration of SKCM patients and revealed consistent clustering results. To gain a new understanding of the relationship between the above grouping and TME phenotype, we conducted a series of survival analyses and found that TME cluster-C (n = 92) was associated with a better prognosis among the three clusters. We explored the potential biological characteristics of the different TME phenotypes and performed a different analysis on the three TME clusters and displayed the intersection of 486 DEGs using a Venn diagram. Based on the above differential genes, we finally obtained 330 repeated DEGs from the four datasets in the GEO database. We used unsupervised clustering to divide SKCM patients into three different subtypes based on the expression of DEGs, namely GeneCluster-A, -B, and -C. According to the expression of DEGs in the different groups, they were divided into Signature gene-A and -B. We confirmed that the expression of these genes was quite different between the two clusters due to various tumor heterogeneities. These clusters were also clearly associated with different mutational patterns. GO enrichment analysis showed that Signature genes-A and -B had different unique BPs. Immune-activating genes are associated with the activation of immune surveillance and immune activation during tumor immunization; Signature gene-A involves the overexpression  (F) The enrichment of different pathway characteristics (immune-related characteristics, mismatch-related characteristics, and matrix-related characteristics) in highand low-risk patients ( * p < 0.05, * * p < 0.01, * * * p < 0.001).
of immune-activated genes, and the gene set overexpressing Signature gene-B is mainly demonstrated by upregulation of genes related to the matrix and transmembrane receptors. Such finding suggests that the Signature gene-A gene set may have an impact on the immunophenotypic landscape of TME. In addition, the analysis results showed that there were significant differences in the expression of TME cells and the enrichment of some related biological pathways among the three GeneCluster groups.
To explore causes of tumor immune microenvironment heterogeneity more deeply, we analyzed the tumor mutationrelated characteristics of patients. The CNV information of SKCM was downloaded in addition to the somatic mutation data and TMB of each patient in Signature gene-A and -B. We also explored SNPs and CNVs in SKCM patients. The gene variable splicing information of SKCM patients in TCGA was also obtained from TCGA Spliceseq website. The results of SNP analysis revealed that there were obvious mutations in both gene FIGURE 7 | Correlation analysis of risk score and clinicopathological characteristics of SKCM patients in TCGA. (A-C) Correlation analysis between risk score and different patient types showed that risk score was related to TME cluster grouping (p < 0.001) (A) and tumor type (p < 0.001) (B). There was no significant correlation with the patient's tumor mutation burden (TMB) (p = 0.6). (D-I) The correlation analysis between the risk score and different clinicopathological characteristics of SKCM patients showed that the risk score and gender (p = 0.01) (D), stage (p = 0.0027) (E), and T stage (p < 0.001) (F) are related, but there is no significant correlation between the patient's age (G), M stage (H), and N stage (I).  To create a better clinical prediction model for different SKCM patients, a novel prognostic-related risk scoring system was developed. First, we incorporated Signature gene-A and -B sets into the univariate Cox model and obtained 233 prognostic-related genes. A subsequent LASSO Cox analysis was performed for dimensionality reduction and model construction, and finally, 13 hub genes were obtained. This novel scoring system had good predictive ability for OS. The AUC of 1year, 3-year, and 5-year OS was 0.713, 0.694, and 0.734, respectively. Furthermore, GSEA was employed to evaluate the correlation between the high-and low-risk groups and biological characteristics. Pathways, such as metabolism-related pathways and oxidative phosphorylation, were mainly enriched in the high-risk group, while immune response-related pathways, including cytokine signaling pathway, JAK-STAT signaling pathway, and natural killer cell-mediated cytotoxicity, were significantly enriched in low-risk patients. This finding is consistent with a previous report . Moreover, the expression of TME cells, as well as some other pathways, such as angiogenesis, mismatch-related characteristics, and matrix-related characteristics, was significantly different between the high-and low-risk groups. Finally, we assessed the correlation between the risk score and the clinical classification and pathological characteristics of SKCM patients and found that low-risk scores were often associated with patients' TME cluster-A classification and metastatic SKCM. In different tumors, 13 hub genes showed different prognostic effects. Furthermore, the results of univariate and multivariate Cox analyses showed that risk scores could be identified as independent risk factors for predicting the prognosis of patients with SKCM. We subsequently selected statistically significant clinical indicators in the multivariate Cox model to construct a nomogram to predict the OS of patients with SKCM. The C-index indicates that the model had a high degree of discrimination. Validation of the calibration curve also revealed good concordance between the estimated values and the actual probability. To our knowledge, this is the first study to explore factors that directly alter the TME in SKCM. Thus, our immune model is more predictable for facilitating treatment than TMB and other predictive factors.
Among the six prognostic signatures, Ubiquitin/ISG15conjugating enzyme E2 L6 (Ube2L6) has been reported to promote insulin resistance and hepatic steatosis and is related to cisplatin resistance (Murakami et al., 2020;Wei et al., 2021). SRPX2 is invasive by upregulating the FAK/SRC/ERK pathway and can lead to pancreatic cancer drug resistance in PI3K/AKT in lung cancer. Additionally, the NFATc3/SRPX2 axis participates in human embryonic stem cell differentiation, indicating that SRPX2 plays a critical role in pan-cancers Li et al., 2020;Chen et al., 2021). SRPX2 is a component of the extracellular matrix, which is important for regulating tumor formation, as demonstrated in diverse tumors, such as colorectal cancer (Øster et al., 2013), gastrointestinal cancer (Tanaka et al., 2012), prostate cancer (Zhang et al., 2018), and pancreatic cancer (Gao et al., 2015). Interferon-induced protein with tetratricopeptide repeats protein may be a new therapeutic target for cancer therapy and can be used as a prognostic marker for cancers, such as glioblastoma and pancreatic cancer (Pidugu et al., 2019). In previous studies, CLEC4E was found to be related to tuberculosis and could be used to constrain tuberculosis through autophagy against drug-resistant strains (Kabuye et al., 2019;Pahari et al., 2020). Killer cell immunoglobulinlike receptor 2DL4 (KIR2DL4) is expressed by NK cells. According to some reports, KIR2DL4 may be an intervention for cancer immunotherapy (Attia et al., 2020). In our IHC analysis, some of these relationships were verified in SKCM patients. Only few studies have examined the relationship between SKCM and these hub genes. Thus, our study revealed their underlying relationships and found a novel important signature.
Therefore, the novel six-gene signature, which was developed to evaluate the comprehensiveness of TME in SKCM, is a reliably developed gene classifier that can be used to predict prognosis and guide more accurate molecular therapy.

CONCLUSION
In conclusion, we have developed and verified an unprecedented set of effective prognostic markers based on immune infiltrating cells, which have certain potential application value in predicting the clinical prognosis of patients with SKCM and the benefit of immunotherapy. This study provides a systematic view of the immune-related characteristics in SKCM and suggests their good prognostic performance.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by The First People's Hospital of Foshan Subject Review Board. The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this manuscript.