Identification of the immune subtypes associated with the prognosis and immunotherapy of metastatic melanoma

Introduction

Melanoma is a highly aggressive cancer, and while accounting for only about 2% of all skin cancer cases, it causes over 72% of the related deaths, largely due to its rapid progression and metastasis to regional lymph nodes and distant organs (1). Metastatic melanoma is associated with poor patient survival outcomes, with its 5-year survival rate being only 5–30% (2,3) and the survival time varying according to each patient. Therefore, developing screening methods and identifying biomarkers with prognostic prediction value in metastatic melanoma have garnered increased attention.

With the recent advancement in immune checkpoint blockade therapies that target programmed cell death ligand 1 (PD-L1), programmed cell death 1 (PD-1), and cytotoxic T lymphocyte-associated protein 4 (CTLA4), patients with metastatic melanoma have received benefit from the immunotherapy. However, only 40–50% of these patients demonstrate clinical benefit and improved prognosis (3-5). Moreover, a portion of patients exhibit drug resistance and adverse effects, such as diarrhea, due to the complexity of interactions between tumor cells and the tumor microenvironment. Consequently, there is an urgent need to establish molecular classifications for metastatic melanoma to identify patients who may respond to immune checkpoint blockade and thus receive a better prognosis.

The tumor microenvironment consists of multiple immune cells and stromal cells, which play an important role in regulating both the emergence and progression of cancer, as well as cellular response to therapies (6-8). Numerous studies have recently confirmed that the infiltration of different types of immune cells can serve as diagnostic and prognostic biomarkers in melanoma (9-11). Research indicates that a high level of immune cell infiltration is associated with favorable clinical outcomes in individuals with metastatic melanoma (9,12,13). Moreover, an increasing number of studies have found that the immune-related features of cancers such as the intensity of CD8⁺ T cell infiltrates, leukocyte fractions, and T-cell receptor (TCR) repertoire are correlated with the immunotherapy response in various cancers, including melanoma (14-20). Thus, multiple changes in the tumor microenvironment of metastatic melanoma may substantially influence and even account for the heterogeneity in clinical outcomes and responses to immunotherapy. Although several investigations have been conducted on the immune subtypes associated with clinical outcomes, their value in predicting response to immunotherapy has not been determined.

The purpose of this study was to systematically characterize the immune cell infiltration in the metastatic melanoma microenvironment and establish novel immune subtypes with prognostic value than can predict the response to immunotherapy. In this study, we classified patients with metastatic melanoma into two distinct subtypes using consensus clustering based on immune-related gene (IRG) expression. The subtypes showed remarkable differences in survival, signaling pathways, gene mutations, immune microenvironment composition, and immunotherapy response. We found that the Immune-High subtype had a higher response rate to immunotherapy. In addition, based on the subtypes of clustering method, we established a powerful model to predict immune subtypes and verified the reliability of the model. We believe that these newly discovered molecular subtypes can be used to predict prognosis and contribute to the development of individualized immunotherapy for patients with metastatic melanoma. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2024-2301/rc).

Methods

Data collection and processing

Microarray and RNA-sequencing gene expression [fragments per kilobase of transcript per million mapped reads (FKPM) and high-throughput sequencing analysis count] of patients with metastatic melanoma were downloaded from the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) database and The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov/), respectively. The clinical, copy number variant (CNV), and somatic mutation data of TCGA were extracted from cBioPortal (PanCancer; https://www.cbioportal.org/). We excluded patients who had received adjuvant or neoadjuvant therapy if information was available. In addition, patients missing necessary clinicopathological or follow-up data were excluded, and samples with a survival time >0 months were retained. Finally, 280 patients in TCGA-human skin cutaneous melanoma (SKCM) project and 188 patients in GSE65904 were included in this study after data filtration for further analysis. All probes of microarray were mapped based on their own Entrez Gene ID (National Center for Biotechnology Information). When multiple probes were matched to the same gene, the probe with maximum average expression value was used to represent the corresponding gene. The log₂-transformed FPKM expression measures were calculated for subsequent analysis. In addition, we recruited an immunotherapy cohort for metastatic melanoma from a previous publication (21), with FASTQ files being obtained from the European Nucleotide Archive (accession no. PRJEB23709).

Level 4 reverse-phase protein array (RPPA) data of TCGA-SKCM patients were obtained from The Cancer Proteome Atlas (TCPA) (22). Data on tumor mutation burden (TMB) and TCR expression were downloaded from a previous study (23). Moreover, human IRGs were obtained from the nCounter PanCancer Immune Profiling Panel (NanoString, Seattle, WA, USA). Various IRGs were included, such as chemokines, cytokine, and genes associated with immune response. The overlapping genes of the GEO dataset, TCGA dataset, and IRGs were selected for further analysis. All the data were obtained by December 2023. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

RNA-sequencing (RNA-seq) analysis

The RNA-seq data of the immunotherapy cohort was analyzed. After removal of adaptor contamination and low-quality reads, clean sequencing reads were aligned via STAR (v. 2.61) software (24) to human reference genome sequence (UCSC hg19 assembly). Gene expression values with FPKM were quantified via RSEM (v. 1.3.0) tool (25) against GENCODE reference annotation (v. 29) (https://www.gencodegenes.org) under default parameters.

Identification of metastatic melanoma subtypes based on IRGs

We applied the “ConsensusClusterPlus” R package (The R Foundation for Statistical Computing) for consensus clustering and molecular subtype screening of IRG expression (26). In brief, k-means clustering was applied, with 1000 iterations (each using 80% of samples) to guarantee the stability of the subtypes. To estimate the optimal number of clusters, the cumulative distribution function (CDF), the relative changes in the area under the CDF curve, and consensus matrix heat maps were evaluated (27). In addition, principal component analysis (PCA) was performed via the “ggbiplot” R package, and hierarchical clustering of gene expression profiles was performed via the “pheatmap” (v. 1.0.12) R package.

Analysis of tumor-infiltrating immune cells

Based on gene expression data, the immune and stromal scores of each sample were calculated via ESTIMATE (Estimation of Stromal and Immune Cells in Malignant Tumor Tissues Using Expression Data), an algorithm based on single-sample gene set enrichment analysis (GSEA) which provides scores for the level of stromal cells present and the infiltration level of immune cells in tumor tissue (28). To estimate the immune cell composition in each sample, Cell-type Identification by Estimating Relative Subsets of known RNA Transcripts (CIBERSORT) was used to quantify the relative proportions of different immune cell types, including B cells, T cells, natural killer (NK) cells, macrophages, and myeloid subsets (29). The algorithm was run with the LM22 signature and 1,000 permutations.

Binomial logistic regression analysis

In order to select a gene set to construct a model capable of dividing samples into Immune-High or Immune-Low subtypes with sensitivity and specificity, we used binomial logistic regression penalized by the least absolute shrinkage and selection operator (LASSO) method. Briefly, we stratified TCGA cohort according to clusters, with 50% of the samples placed in a training set and the remainder in a test set. We developed the model by training in TCGA training set and validating it in TCGA test set and GSE65904 cohort. The expression of IRGs was first z-score normalized. Subsequently, we performed logistic regression using the cv.glmnet function in the “glmnet” R package with 10-fold cross-validation to identify genes and develop the model (30). After inclusion of the z-score expression value of the selected genes, a score formula for each sample was established and weighted according to the respective estimated regression coefficients. The score of sample i was calculated with the following formula:

Scorei=∑g=1nβg∗Xgi[1]

where n is the number of selected genes, βg is the coefficient of selected gene g, and Xgi is to the z-score expression value of gene g in the sample i. If the score for sample i was higher than the cutoff, the sample was considered to be the Immune-Low subtype. In this study, for TCGA test set and GSE65904 cohort, we used the Youden index with the receiver operating characteristic curve (ROC) and median score in the immunotherapy cohort to determine the cutoff value.

Differentially expressed genes (DEGs) and gene set enrichment analysis

To identify DEGs between clusters, we used the “limma” and “DESeq2” R packages for GSE65904 and TCGA, respectively (31,32). The criteria for selecting DEGs were an adjusted P value <0.05 and fold change (FC) ≥1.5 or ≤0.5 for GSE65904, and for TCGA, they were an adjusted P value <0.05 and |log₂FC| >1. The DEGs common to both cohorts were selected as the final DEGs. Subsequently, Gene Ontology (GO) biological-process enrichment analysis and visualization were performed via the “clusterProfiler” R package (26). In addition, GSEA was performed for the two subtypes via the GSEA Desktop application in Java. The reference Hallmark gene sets (h.all.v7.2.entrez.gmt) were obtained from the Molecular Signatures Database (MSigDB; http://software.broadinstitute.org/gsea/msigdb/). The number of permutations was set to 1,000, and the threshold was a false-discovery rate <0.05.

Statistical analysis

All statistical analyses in this study were performed with R version 3.62. For comparison of continuous variable, such as gene expression, TMB, and TCR expression, the Wilcoxon test was applied. For categorical variables between groups (e.g., response rate), we used the Fisher exact test or Chi-squared test. Survival analyses were performed via the Kaplan-Meier method, with survival being compared via the log-rank test in the “survival” R package and the curves being plotted via the “survminer” R package. The ROC was plotted and the area under the curve (AUC) calculated via the “pROC” R package (33). Finally, the somatic gene mutations in the different subtypes were determined via the “maftools” R package (34). All statistical tests were two-sided, and a P value less than 0.05 was considered statistically significant.

Results

Immune clusters in metastatic melanoma

A total of 280 and 188 patients with metastatic melanoma with survival information from TCGA and GEO were included in this study, respectively. In order to group patients with metastatic melanoma according to the expression of the 709 IRGs (available online: https://cdn.amegroups.cn/static/public/10.21037tcr-2024-2301-1.xlsx), we performed consensus unsupervised clustering method on TCGA and GSE65904 cohorts. We considered the relative change in area under the CDF curve from a ranking of 2 to 10, and found that the consensus heatmaps showed stable partitions when the samples were clustered into two groups. Of the 280 metastatic melanoma samples in TCGA, 128 samples were classified into cluster1 group and 152 samples were classified into cluster2 group (Figure 1A,1B). To access the subclasses’ assignments, we performed PCA. The data showed that the two clusters were distributed in different parts of the two-dimensional coordinate systems (Figure 1C). Furthermore, we extracted the expression data of IRGs in the GSE65904 cohort with 188 samples from GEO to analyze the clusters. Consistent with the two subtypes of TCGA, the two clusters were also identified in the GSE65904 cohort, with 90 samples in cluster1 and 98 samples in cluster2 group (Figure 1D,1E). PCA of GSE65904 also confirmed the two subtypes identified by the consensus unsupervised clustering method (Figure 1F).

Figure 1 Identification of two immune subtypes in metastatic melanoma. (A) CDF curve of the consistency score for different subtype numbers (k=2–10) in TCGA cohort. (B) The consensus score matrix of all samples when k=2 in TCGA cohort. (C) PCA for the immune-gene expression profiles showing a significant difference between the two subtypes in TCGA cohort. (D) CDF curve in the GSE65904 cohort. (E) The consensus score matrix in the GSE65904 cohort. (F) PCA plot of the GSE65904 cohort. CDF, cumulative distribution function; PCA, principal component analysis; TCGA, The Cancer Genome Atlas.

To confirm that the clusters were associated with the tumor immune microenvironment, we used the ESTIMATE algorithm to calculate the immune and stromal scores, which reflect the level of infiltration into tumors by immune cells. We found that the two clusters significantly correlated with immune score and stromal score in the two cohorts, respectively (Figure 2A-2D). Consequently, we defined the cluster with low immune score as the Immune-Low subtype and the cluster with the high immune score as the Immune-High subtype and then further compared the clinical characteristics of the two subtypes. As shown in Table 1, the survival status was significantly different between the two subtypes in TCGA (Chi-squared test, P<0.001) but not quite significant in GSE65904. Based on these results, we concluded that an unsupervised clustering using IRGs via the nCounter panel could classify metastatic melanoma according to the level of immune infiltration. Moreover, tumors from the different subtypes exhibited considerable differences in their immune status.

Figure 2 Association between the tumor microenvironment and immune subtypes. The comparison of the immune score (A) and stromal score (B) in TCGA cohort, respectively, and the immune score (C) and stromal score (D) in the GSE65904 cohort, respectively, between two subtypes. (E,F) Kaplan Meier survival curves for overall survival between the two subtypes in TCGA and GSE65904 cohort, respectively. ****, P<0.0001. TCGA, The Cancer Genome Atlas.

Association of immune clusters with prognosis

To examine the potential correlation between prognosis and immune clusters, we constructed Kaplan-Meier survival curves of patients with metastatic melanoma from each subtype. The results suggested that there were considerable differences in overall survival (OS) between the two subtypes. The patients from the Immune-High subtype had better OS compared to those from the Immune-Low subtype for both TCGA (P<0.001) and GSE65904 (P=0.006) (Figure 2E,2F). For TCGA cohort, the 5-year OS rates were 46.58% (95% CI: 38.82–55.90%) in the Immune-Low subtype and 71.8% (95% CI: 63.7–80.8%) in the Immune-High subtype. For the GSE65904 cohort, the 5-year OS rate was 30.5% (95% CI: 20.7–45.1%) in the Immune-Low subtype and 42.5% (95% CI: 31.7–57.1%) in the Immune-High subtype. Therefore, in both cohorts, the immune-high subtype had longer 5-year OS, indicating that the immune-associated subtypes may be associated with prognosis.

The immune landscape and microenvironment of metastatic melanoma

Given the significant difference in immune score identified between the subtypes, we used the CIBERSORT algorithm to characterize the immune infiltration and immunologic landscape of these subtypes. We identified 11 types of immune cell populations that significantly differed between the two subtypes. Compared with the immune-low subtype, the immune-high subtype had greater proportions of naïve B cells, plasma cells, CD8 T cells, memory activated CD4 T cells, activated NK cells, and M1 macrophages. Meanwhile, the proportions of resting NK cells, M0 macrophages, M2 macrophages, resting mast cells, and eosinophils were greater in the immune-low subtype (Figure 3A and Figure S1A). We further investigated the association between subtypes and the expression of several potentially targetable immune checkpoints and MHC I/II class genes. We found that compared with the Immune-Low subtype, the Immune-High subtype exhibited significantly increased expression of PD-1, PD-L1, CTLA4, and IDO1, among other immune checkpoint genes (Figure 3B and Figure S1B). Moreover, the Immune-High subtype had a significantly increased PD-L1 reverse phase protein array (RPPA) level (Figure S2). The results suggest that the Immune-High subtype has a higher proportion of antitumor immune cells and a lower proportion of tumor-promoting immune cells as compared to the Immune-Low subtype.

Figure 3 Immune cellular and molecular features of immune subtypes in TCGA cohort. (A) Differential proportions of immune cells between the Immune-High and Immune-Low subtypes. (B) The expression of immune checkpoint molecules and MHC I/II class genes. *, P<0.05; ****, P<0.0001. TCGA, The Cancer Genome Atlas.

Comparison of TMB, CNV, and TCR levels between the Immune-High and Immune-Low subtypes

To determine the reason for cluster heterogeneity, we examined how the two subtypes differed in terms of DNA mutations. The waterfall plot of key melanoma-associated mutated genes across clusters in Figure 4A indicates that the subtypes tended to have different mutation proportions of the NF1 gene: 24% percent of samples in the Immune-Low subtype had the NF1 gene mutation, while only 14% of samples in the Immune-High subtype had this mutation (Fisher exact test, P=0.048). We compared the TMB between the Immune-High and Immune-Low subtypes, but no difference was found (Figure 4B).

Figure 4 Landscape of the mutation profile in TCGA samples. (A) The waterfall plot the important melanoma-associated mutated genes of DNA somatic mutations in the two subtypes. The comparison of the tumor mutation burden (B), the number of CNV (C), the number of CNV amplifications (D), the number of CNV deletions (E), and TCR repertoire diversity (F) between the two subtypes. ****, P<0.0001. CNV, copy number variant; TCGA, The Cancer Genome Atlas; TCR, T-cell receptor.

Previous studies indicated that an increased number of CNVs in melanoma is inversely correlated with immune cell infiltration (35-37). To determine whether CNV abundance may lead to increased heterogeneity between the immune-associated subtypes, we analyzed the number of CNVs based on segments. We found that the number of CNV segments was higher in the Immune-Low subtype than in the Immune-High subtype (Wilcox test, P<0.001; Figure 4C). Similarly, the number of amplifications and deletions were also higher in the Immune-Low subtype (Figure 4D,4E). These results suggest that CNVs are significantly different between the two subtypes.

Another previous study indicated that the immune-high tumors have higher TCR repertoire diversity (38). Therefore, we next examined the relationship between the immune subtypes and TCR repertoire diversity and found that the Immune-High subtype had higher TCR repertoire diversity (Figure 4F).

Therefore, it appears that the Immune-High subgroup and Immune-Low subgroup correspond to an immune-active phenotype and immune-resting phenotype, respectively.

Identification of DEGs between the two subtypes

To better understand the difference in molecular characteristics between the Immune-High and Immune-Low subtypes, it was necessary to determine specific genes that are enriched in each subtype. Differential expression analysis between the Immune-High and Immune-Low subtypes identified 2,751 and 1,014 DEGs in TCGA and GSE65904, respectively, with 533 DEGs common to the two cohorts (Figure 5A and available online: https://cdn.amegroups.cn/static/public/10.21037tcr-2024-2301-2.xlsx). Interestingly, most of these DEGs were upregulated in the Immune-High subtype. Furthermore, 184 immune genes in the nCounter panel were identified as the DEGs between the two subtypes. The functions of the intersecting 533 genes were analyzed. GO enrichment results revealed that these genes were mainly enriched in immune response, such as T-cell activation, regulation of lymphocyte activation, leukocyte cell-cell adhesion, and leukocyte proliferation (Figure 5B and available online: https://cdn.amegroups.cn/static/public/10.21037tcr-2024-2301-3.xlsx).

Figure 5 DEG and enrichment analyses of the two subtypes. (A) The Venn diagram of the DEGs and immune genes. (B) The top 25 most enriched GO biological processes from the 533 DEGs. (C) The enriched hallmark pathways based on gene set enrichment analysis in TCGA cohort. DEG, differentially expressed gene; GO, Gene Ontology; TCGA, The Cancer Genome Atlas.

GSEA was also used to identify the mechanism and functional differences between the Immune-High and Immune-Low subtypes. The Immune-High subtype was enriched in the following pathways: IL6/JAK/STAT3 signaling; inflammatory response; KRAS signaling up, complement; PI3K/AKT/MTOR signaling; allograft rejection; TNF-α signaling via NF-Kβ; IL2/STAT5 signaling; IFN-α response; and IFN-γ response. Meanwhile, the Immune-Low subtype was enriched in MYC targets v1 pathway (Figure 5C and Table S1). These findings suggest that the signaling pathways from the Immune-High subgroup are closely associated with immune infiltration.

Prediction of immune clusters with binomial logistic regression

Because the unsupervised clustering method depends on cohort studies with large samples, it cannot be used to predict the classification of a single sample to the immune clusters. Therefore, binomial logistic regression penalized via the LASSO method was applied to derive a gene signature for predicting metastatic melanoma immune clusters. We stratified TCGA cohort samples according to clusters, with 50% of samples being placed in the training set and the remaining samples placed in the test set. We developed a model by training in TCGA training set, and validated in TCGA test set and GSE65904 cohort. From the immune genes, a 36-gene signature was identified (Figure 6A, Table S2). To assess the specificity and sensitivity of binomial logistic regression model in predicting the immune clusters, we used ROC curve and AUC analysis. As shown in the Figure 6B,6C, the model predicted immune clusters, with AUCs of 0.984 and 0.953 in TCGA test set and GSE65904 cohort, respectively. According to the cutoff of the ROC (Youden index), we divided the samples into Immune-High and Immune-Low subtypes. Heatmaps of the two-way hierarchical clustering analysis of the 36 genes are shown in Figure S3. From the heatmaps, we found that the expression values of the 36 genes could separate the different subtypes of samples, indicating that the identified 36-gene signature can be used to distinguish metastatic melanoma Immune-High samples from Immune-Low samples. As expected, patient survival in the Immune-High subtype predicted by the LASSO model was also significantly longer than that in the Immune-Low subtype (Figure 6D,6E). Hence, these results demonstrate that the binomial logistic regression model can provide a single-sample predictor that can be applied to predict immune clusters and serve as an independent prognostic factor in patients with metastatic melanoma.

Figure 6 Construction and validation of the binomial logistic regression model. (A) Binomial deviance of the different numbers of variables revealed by the binomial logistic regression penalized by the LASSO method. The red dots represent the binomial deviance values, and the gray lines represent the SE. The two vertical dotted lines on the left and right, respectively, represent optimal values use by the minimum criterion and the 1−SE criterion. λ is the tuning parameter. (B,C) ROC curve analysis of the binomial logistic regression model for predicting sample immune subtypes in TCGA training set and GSE65904 cohort, respectively. The data in parentheses are 95% confidence intervals. (D,E) Kaplan-Meier survival curves for the overall survival of the two subtypes as predicted by the binomial logistic regression model in TCGA training set and the GSE65904 cohort, respectively. AUC, area under the curve; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic; SE, standard error; TCGA, The Cancer Genome Atlas.

Immune clusters and response to immunotherapy

Considering the differences in immune infiltration patterns and expression levels of immune checkpoint genes between the metastatic melanoma subtypes, we investigated whether our model could predict patient responses to immune checkpoint blockade therapy based on the 41 patients with metastatic melanoma that received anti-PD-1 treatment (21). According to the median of the established score, patients were divided into two groups (Immune-High and Immune-Low). As expected, the Immune-High subtype had a significantly higher immune score than did the Immune-Low subtype (Figure 7A). In addition, the messenger RNA expression of PD-1 and PD-L1 was higher in the Immune-High subtype (Figure 7B,7C). The proportion of response (complete response or partial response) to anti-PD-1 treatment in the Immune-Low (28.6%) was lower than that of the Immune-High subtype (65%) (Fisher exact test, P=0.029; Figure 7D). We also performed Kaplan-Meier survival analysis and found that patients with metastatic melanoma in the Immune-Low group, as compared to those in the Immune-High group, had shorter progression-free survival (PFS) (P=0.007) and OS (P=0.034) (Figure 7E,7F). These results suggest that patients of the Immune-High subtype may receive greater benefit from immunotherapy than patients of the Immune-Low subtype.

Figure 7 Immunotherapy analysis of the immune subtypes. (A-C) The comparison of the immune score, PD-1 expression, and PD-L1 expression in the immunotherapy cohort. (D) The proportion of patients with response to PD-1 blockade immunotherapy in the Immune-High and Immune-Low subtypes. (E,F) Kaplan-Meier survival curves for progression-free survival and overall survival of the two subtypes predicted by the binomial logistic regression model, respectively. ***, P<0.001; ****, P<0.0001. CR, complete response; PD, progressive disease; PD-1, programmed cell death 1; PR, partial response; SD, stable disease.

Discussion

Metastatic melanoma is one of most aggressive and lethal types of cancer. Immune heterogeneity with the tumor microenvironment has been confirmed to be associated with the prognosis and drug response of patients with various cancers, including melanoma (23,39,40). Emerging evidence indicates that high immune cell infiltration is associated with improved prognosis in different cancers (39,41). Therefore, in our study, we aimed to identify IRGs that may be associated with the OS and immunotherapy response of patients with metastatic melanoma by investigating the tumor microenvironment.

We first identified novel immune subtypes based on IRGs. The immune infiltration, clinical characteristics, prognosis, DNA mutations, CNV, TCR repertoire diversity, and immunotherapy response were then investigated in the different subtypes. Cluster 1 exhibited higher immune score and stromal score and was defined as the Immune-High subtype, while cluster 2 was delineated as the Immune-Low subtype. The patients from the Immune-High subtype had better prognosis than did those from the Immune-Low subtype. We identified two distinct subtypes of metastatic melanoma which may help inform the selection of patient populations responsive to immunotherapy.

We used CIBERSORT to estimate the immune-cell composition of TCGA and GSE65904 cohorts and compared the different types of cells in the Immune-High and Immune-Low subtypes. Surprisingly, it was found that the proportions of naïve B cells, plasma cells, CD8 T cells, memory activated CD4 T cells, activated NK cells, and M1 macrophages in the Immune-High subtype were higher than those in the Immune-Low subtype. Meanwhile, the proportion of resting NK cells, M0 macrophages, M2 macrophages, resting mast cells, and eosinophils were lower in the Immune-High subtype. Tumor-associated macrophages are macrophages that infiltrate into malignant tumor tissues and are mainly derived from circulating blood monocytes released from the bone marrow (42). M1 macrophages are activated by cytokines such as IFN-γ, produce proinflammatory cytokines, and are involved in inflammatory response, pathogen clearance, and antitumor immunity (43,44). M2 macrophages are well-known to promote tumor proliferation, and many studies have revealed that M2 macrophages are associated with a poor prognosis in different cancer types (45,46). In recent years, researchers have demonstrated that M2 macrophages contribute to poor prognosis in melanoma (47-49). As expected, in our study, we found that the Immune-High subtype generally had a higher abundance of CD8⁺ T cells than did the Immune-Low subtype. Cabrita et al. reported that CD8⁺ T-cell infiltration was associated with improved outcomes in patients with melanoma (12). Erdag et al. also found that patients with metastatic melanoma with higher CD8⁺ T-cell infiltration had a better prognosis (9). Moreover, it has been confirmed that higher CD8⁺ T cell trafficking is critical to anti-PD-1/anti-CTLA-4 efficacy in patients with brain metastases from melanoma (20,50). These studies support our findings regarding the immune composition and prognosis of clusters, confirming the reliability and feasibility of our research.

The TMB is an important index for reflecting the accumulation of somatic mutations and is considered to be a biomarker for the response to immune checkpoint inhibitors (51,52). Although no significant difference was found between the two groups, further study based on large samples to examine the association between the immune of metastatic melanoma and TMB is still needed.

We found that the NF1 mutation count was higher in the Immune-Low subtype. A previous report indicated that patients with melanoma harboring the NF1 mutation have poor overall survival (53). Other than TMB, CNV has been found to negatively correlate with immune-cell infiltration in melanoma, and patients with a higher CNV count are more likely to be resistant to immunotherapy and have poor survival (35-37). Our findings revealed that there was a significant difference in CNV between the two subtypes. In addition, the Immune-Low subtype had a higher number of amplifications and deletions. This suggests that increased CNV count is a key factor contributing to the differences between the two melanoma subtypes.

Additionally, the Immune-Low subtype was linked to the a lower level of TCR repertoire diversity, which is consistent with a previous study (38). TCRs are the antigen recognition structure physiologically expressed on the surface of all T cells. The type of TCRs on each T cell are unique and vary across individuals, resulting in different immune responses to a diverse range of neoantigens (54). The diversity of the TCR repertoire reflects the difference om cellular immunity, and several studies have shown that TCR diversity is associated with the immunotherapy and prognosis of patients with melanoma (55-58). Although there was no difference in the TMB between two groups in our study, the difference in TCR diversity is significantly. Moreover, we found that the Immune-High subtype had an elevated expression of HLA genes, which are linked to immune activation, immunotherapy response, and better prognosis in numerous malignancies (59,60).

GSEA and DEG and GO enrichment analyses were used to clarify the underlying mechanisms related to the immune differences in the subtypes of metastatic melanoma. A number of immune genes, immune-related functions, and pathways were significantly different between two groups. We found that IL6/JAK/STAT3, PI3K/AKT/MTOR signaling, etc., were enriched in the Immune-High subtype. Numerous studies have found that signaling via the IL-6/JAK/STAT3 pathway induces the expression of PD-1 and/or PD-L1 (61-63). In addition, IFN-γ pathway has an immunosuppressive role through the upregulation of immune-checkpoint molecules including CTLA4, PD-1, PD-L1, PD-L2, and IDO in melanoma (64-66). IFN-γ activates both dendritic cells and macrophages to enhance antigen presentation (67). Moreover, INF-γ signaling in tumor cells can also facilitate the immune recognition and apoptosis of tumor cells, thus modulating the resistance of tumor cells to immune response (68,69). Activation of the PI3K/AKT/MTOR signaling also increases PD-L1 expression (70,71). In our study, we discovered that the immune checkpoint molecules were enhanced in the Immune-High subtype. These results indicate that the Immune-High subtype might have a higher expression of immune checkpoint molecules and superior response to immunotherapy.

Furthermore, we applied the binomial logistic regression penalized by the LASSO to construct a model that predict immune clusters in metastatic melanoma. Through validation in test sets, our model was confirmed to have strong robustness in identifying the immune subtypes of melanoma and evaluating patient prognosis. It is well known that CD8⁺ T-cell infiltration, PD-1/PD-L1 expression, and immune hot/cold status are the key factors influencing melanoma response to immunotherapy (40,72). According to the above-mentioned results, we speculate that the Immune-High subtype might have a greater response to immunotherapy. Interestingly, we used the model to predict the subtype, and the results suggested that patients in the Immune-High subtype may benefit more from anti–PD-1 treatment.

Although there have been several previous studies on the tumor microenvironment of melanoma, they focused mainly on tumor-infiltrating lymphocytes (73) or DNA methylation (74). Therefore, our study possesses a degree of novelty and certain advantages. First, the sample sizes of the cohorts derived from TCGA and GEO databases were sufficiently large, ensuring the accuracy of findings. Second, we investigated the heterogeneity of metastatic melanoma with integrated multiomics data types, including genomics, transcriptomics, proteomics, and clinical data. Third, the analysis in our study was conducted via various well-known machine learning methods in clear and orderly fashion, resulting in the construction of an immune-subtype prediction model, which can be used for survival and immunotherapy response prediction.

However, our study also involved certain limitations. First, we only applied bioinformatics research, and the results lacked verification via experimental data. Further experiments are required to confirm the findings of this study. Second, the sample size was small, and the validation was performed only in anti-PD-1 treated patients, with no assessment in those receiving anti-CTLA4 or anti-PD-L1 patients. Melanoma patients treated with other immune checkpoint blockades are need to be evaluated.

Conclusions

We identified two immune subtypes of metastatic melanoma and generated a prognostic model, which might have important clinical implications for predicting patients’ survival outcome and improving the efficacy of immunotherapy.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2024-2301/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2024-2301/prf

Funding: None.

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2024-2301/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2020. CA Cancer J Clin 2020;70:7-30. [Crossref] [PubMed]
2. Weiss SA, Hanniford D, Hernando E, et al. Revisiting determinants of prognosis in cutaneous melanoma. Cancer 2015;121:4108-23. [Crossref] [PubMed]
3. van Not OJ, van den Eertwegh AJM, Jalving H, et al. Long-Term Survival in Patients With Advanced Melanoma. JAMA Netw Open 2024;7:e2426641. [Crossref] [PubMed]
4. Larkin J, Chiarion-Sileni V, Gonzalez R, et al. Five-Year Survival with Combined Nivolumab and Ipilimumab in Advanced Melanoma. N Engl J Med 2019;381:1535-46. [Crossref] [PubMed]
5. Noringriis IM, Donia M, Madsen K, et al. Long-term clinical outcome of patients with metastatic melanoma and initial stable disease during anti-PD-1 checkpoint inhibitor immunotherapy with pembrolizumab. Br J Cancer 2025;133:337-45. [Crossref] [PubMed]
6. Hui L, Chen Y. Tumor microenvironment: Sanctuary of the devil. Cancer Lett 2015;368:7-13. [Crossref] [PubMed]
7. de Visser KE, Joyce JA. The evolving tumor microenvironment: From cancer initiation to metastatic outgrowth. Cancer Cell 2023;41:374-403. [Crossref] [PubMed]
8. Glaviano A, Lau HS, Carter LM, et al. Harnessing the tumor microenvironment: targeted cancer therapies through modulation of epithelial-mesenchymal transition. J Hematol Oncol 2025;18:6. [Crossref] [PubMed]
9. Erdag G, Schaefer JT, Smolkin ME, et al. Immunotype and immunohistologic characteristics of tumor-infiltrating immune cells are associated with clinical outcome in metastatic melanoma. Cancer Res 2012;72:1070-80. [Crossref] [PubMed]
10. Yang S, Liu T, Nan H, et al. Comprehensive analysis of prognostic immune-related genes in the tumor microenvironment of cutaneous melanoma. J Cell Physiol 2020;235:1025-35. [Crossref] [PubMed]
11. Bagaev A, Kotlov N, Nomie K, et al. Conserved pan-cancer microenvironment subtypes predict response to immunotherapy. Cancer Cell 2021;39:845-865.e7. [Crossref] [PubMed]
12. Cabrita R, Lauss M, Sanna A, et al. Tertiary lymphoid structures improve immunotherapy and survival in melanoma. Nature 2020;577:561-5. [Crossref] [PubMed]
13. Connor C, Carr QL, Sweazy A, et al. Clinical Approaches for the Management of Skin Cancer: A Review of Current Progress in Diagnosis, Treatment, and Prognosis for Patients with Melanoma. Cancers (Basel) 2025;17:707. [Crossref] [PubMed]
14. Tang H, Wang Y, Chlewicki LK, et al. Facilitating T Cell Infiltration in Tumor Microenvironment Overcomes Resistance to PD-L1 Blockade. Cancer Cell 2016;30:500. [Crossref] [PubMed]
15. Li B, Cui Y, Diehn M, et al. Development and Validation of an Individualized Immune Prognostic Signature in Early-Stage Nonsquamous Non-Small Cell Lung Cancer. JAMA Oncol 2017;3:1529-37. [Crossref] [PubMed]
16. Zeng D, Li M, Zhou R, et al. Tumor Microenvironment Characterization in Gastric Cancer Identifies Prognostic and Immunotherapeutically Relevant Gene Signatures. Cancer Immunol Res 2019;7:737-50. [Crossref] [PubMed]
17. Riaz N, Havel JJ, Makarov V, et al. Tumor and Microenvironment Evolution during Immunotherapy with Nivolumab. Cell 2017;171:934-949.e16. [Crossref] [PubMed]
18. Zhao Y, Schaafsma E, Gorlov IP, et al. A Leukocyte Infiltration Score Defined by a Gene Signature Predicts Melanoma Patient Prognosis. Mol Cancer Res 2019;17:109-19. [Crossref] [PubMed]
19. Gobbini E, Hubert M, Doffin AC, et al. The Spatial Organization of cDC1 with CD8+ T Cells is Critical for the Response to Immune Checkpoint Inhibitors in Patients with Melanoma. Cancer Immunol Res 2025;13:517-26. [Crossref] [PubMed]
20. Ellingsen EB, Bounova G, Kerzeli I, et al. Characterization of the T cell receptor repertoire and melanoma tumor microenvironment upon combined treatment with ipilimumab and hTERT vaccination. J Transl Med 2022;20:419. [Crossref] [PubMed]
21. Gide TN, Quek C, Menzies AM, et al. Distinct Immune Cell Populations Define Response to Anti-PD-1 Monotherapy and Anti-PD-1/Anti-CTLA-4 Combined Therapy. Cancer Cell 2019;35:238-255.e6. [Crossref] [PubMed]
22. Li J, Lu Y, Akbani R, et al. TCPA: a resource for cancer functional proteomics data. Nat Methods 2013;10:1046-7. [Crossref] [PubMed]
23. Thorsson V, Gibbs DL, Brown SD, et al. The Immune Landscape of Cancer. Immunity 2018;48:812-830.e14. [Crossref] [PubMed]
24. Dobin A, Davis CA, Schlesinger F, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 2013;29:15-21. [Crossref] [PubMed]
25. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 2011;12:323. [Crossref] [PubMed]
26. Yu G, Wang LG, Han Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 2012;16:284-7. [Crossref] [PubMed]
27. Young JD, Cai C, Lu X. Unsupervised deep learning reveals prognostically relevant subtypes of glioblastoma. BMC Bioinformatics 2017;18:381. [Crossref] [PubMed]
28. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. [Crossref] [PubMed]
29. Chen B, Khodadoust MS, Liu CL, et al. Profiling Tumor Infiltrating Immune Cells with CIBERSORT. Methods Mol Biol 2018;1711:243-59. [Crossref] [PubMed]
30. Friedman J, Hastie T, Tibshirani R. Regularization Paths for Generalized Linear Models via Coordinate Descent. J Stat Softw 2010;33:1-22.
31. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47. [Crossref] [PubMed]
32. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 2014;15:550. [Crossref] [PubMed]
33. Robin X, Turck N, Hainard A, et al. pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics 2011;12:77. [Crossref] [PubMed]
34. Mayakonda A, Lin DC, Assenov Y, et al. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res 2018;28:1747-56. [Crossref] [PubMed]
35. Davoli T, Uno H, Wooten EC, et al. Tumor aneuploidy correlates with markers of immune evasion and with reduced response to immunotherapy. Science 2017;355:eaaf8399. [Crossref] [PubMed]
36. Asleh K, Ouellette RJ. Tumor Copy Number Alteration Burden as a Predictor for Resistance to Immune Checkpoint Blockade across Different Cancer Types. Cancers (Basel) 2024;16:732. [Crossref] [PubMed]
37. Yang X, Hu Y, Yang K, et al. Cell-free DNA copy number variations predict efficacy of immune checkpoint inhibitor-based therapy in hepatobiliary cancers. J Immunother Cancer 2021;9:e001942. [Crossref] [PubMed]
38. Song Y, Yan S, Fan W, et al. Identification and Validation of the Immune Subtypes of Lung Adenocarcinoma: Implications for Immunotherapy. Front Cell Dev Biol 2020;8:550. [Crossref] [PubMed]
39. Cesano A, Warren S. Bringing the next Generation of Immuno-Oncology Biomarkers to the Clinic. Biomedicines 2018;6:14. [Crossref] [PubMed]
40. Byrne EH, Fisher DE. Immune and molecular correlates in melanoma treated with immune checkpoint blockade. Cancer 2017;123:2143-53. [Crossref] [PubMed]
41. He B, Gao R, Lv D, et al. The prognostic landscape of interactive biological processes presents treatment responses in cancer. EBioMedicine 2019;41:120-33. [Crossref] [PubMed]
42. Komohara Y, Fujiwara Y, Ohnishi K, et al. Tumor-associated macrophages: Potential therapeutic targets for anti-cancer therapy. Adv Drug Deliv Rev 2016;99:180-5. [Crossref] [PubMed]
43. Yang L, Zhang Y. Tumor-associated macrophages: from basic research to clinical application. J Hematol Oncol 2017;10:58. [Crossref] [PubMed]
44. Chanmee T, Ontong P, Konno K, et al. Tumor-associated macrophages as major players in the tumor microenvironment. Cancers (Basel) 2014;6:1670-90. [Crossref] [PubMed]
45. Yamaguchi T, Fushida S, Yamamoto Y, et al. Tumor-associated macrophages of the M2 phenotype contribute to progression in gastric cancer with peritoneal dissemination. Gastric Cancer 2016;19:1052-65. [Crossref] [PubMed]
46. Dijkgraaf EM, Heusinkveld M, Tummers B, et al. Chemotherapy alters monocyte differentiation to favor generation of cancer-supporting M2 macrophages in the tumor microenvironment. Cancer Res 2013;73:2480-92. [Crossref] [PubMed]
47. Bronkhorst IH, Ly LV, Jordanova ES, et al. Detection of M2-macrophages in uveal melanoma and relation with survival. Invest Ophthalmol Vis Sci 2011;52:643-50. [Crossref] [PubMed]
48. López-Janeiro Á, Padilla-Ansala C, de Andrea CE, et al. Prognostic value of macrophage polarization markers in epithelial neoplasms and melanoma. A systematic review and meta-analysis. Mod Pathol 2020;33:1458-65. [Crossref] [PubMed]
49. Nasir I, McGuinness C, Poh AR, et al. Tumor macrophage functional heterogeneity can inform the development of novel cancer therapies. Trends Immunol 2023;44:971-85. [Crossref] [PubMed]
50. Taggart D, Andreou T, Scott KJ, et al. Anti-PD-1/anti-CTLA-4 efficacy in melanoma brain metastases depends on extracranial disease and augmentation of CD8(+) T cell trafficking. Proc Natl Acad Sci U S A 2018;115:E1540-9. [Crossref] [PubMed]
51. Alexandrov LB, Nik-Zainal S, Wedge DC, et al. Signatures of mutational processes in human cancer. Nature 2013;500:415-21. [Crossref] [PubMed]
52. Van Allen EM, Miao D, Schilling B, et al. Genomic correlates of response to CTLA-4 blockade in metastatic melanoma. Science 2015;350:207-11. [Crossref] [PubMed]
53. Cirenajwis H, Lauss M, Ekedahl H, et al. NF1-mutated melanoma tumors harbor distinct clinical and biological characteristics. Mol Oncol 2017;11:438-51. [Crossref] [PubMed]
54. Liu YY, Yang QF, Yang JS, et al. Characteristics and prognostic significance of profiling the peripheral blood T-cell receptor repertoire in patients with advanced lung cancer. Int J Cancer 2019;145:1423-31. [Crossref] [PubMed]
55. Lin KR, Deng FW, Jin YB, et al. T cell receptor repertoire profiling predicts the prognosis of HBV-associated hepatocellular carcinoma. Cancer Med 2018;7:3755-62. [Crossref] [PubMed]
56. Postow MA, Manuel M, Wong P, et al. Peripheral T cell receptor diversity is associated with clinical outcomes following ipilimumab treatment in metastatic melanoma. J Immunother Cancer 2015;3:23. [Crossref] [PubMed]
57. Cha E, Klinger M, Hou Y, et al. Improved survival with T cell clonotype stability after anti-CTLA-4 treatment in cancer patients. Sci Transl Med 2014;6:238ra70. [Crossref] [PubMed]
58. Robert L, Harview C, Emerson R, et al. Distinct immunological mechanisms of CTLA-4 and PD-1 blockade revealed by analyzing TCR usage in blood lymphocytes. Oncoimmunology 2014;3:e29244. [Crossref] [PubMed]
59. Johnson DB, Estrada MV, Salgado R, et al. Melanoma-specific MHC-II expression represents a tumour-autonomous phenotype and predicts response to anti-PD-1/PD-L1 therapy. Nat Commun 2016;7:10582. [Crossref] [PubMed]
60. Rodig SJ, Gusenleitner D, Jackson DG, et al. MHC proteins confer differential sensitivity to CTLA-4 and PD-1 blockade in untreated metastatic melanoma. Sci Transl Med 2018;10:eaar3342. [Crossref] [PubMed]
61. Bu LL, Yu GT, Wu L, et al. STAT3 Induces Immunosuppression by Upregulating PD-1/PD-L1 in HNSCC. J Dent Res 2017;96:1027-34. [Crossref] [PubMed]
62. Zhang N, Zeng Y, Du W, et al. The EGFR pathway is involved in the regulation of PD-L1 expression via the IL-6/JAK/STAT3 signaling pathway in EGFR-mutated non-small cell lung cancer. Int J Oncol 2016;49:1360-8. [Crossref] [PubMed]
63. Atsaves V, Tsesmetzis N, Chioureas D, et al. PD-L1 is commonly expressed and transcriptionally regulated by STAT3 and MYC in ALK-negative anaplastic large-cell lymphoma. Leukemia 2017;31:1633-7. [Crossref] [PubMed]
64. Mo X, Zhang H, Preston S, et al. Interferon-γ Signaling in Melanocytes and Melanoma Cells Regulates Expression of CTLA-4. Cancer Res 2018;78:436-50. [Crossref] [PubMed]
65. Spranger S, Spaapen RM, Zha Y, et al. Up-regulation of PD-L1, IDO, and T(regs) in the melanoma tumor microenvironment is driven by CD8(+) T cells. Sci Transl Med 2013;5:200ra116. [Crossref] [PubMed]
66. Bai X, Fisher DE, Flaherty KT. Cell-state dynamics and therapeutic resistance in melanoma from the perspective of MITF and IFNγ pathways. Nat Rev Clin Oncol 2019;16:549-62. [Crossref] [PubMed]
67. Minn AJ. Interferons and the Immunogenic Effects of Cancer Therapy. Trends Immunol 2015;36:725-37. [Crossref] [PubMed]
68. Gao J, Shi LZ, Zhao H, et al. Loss of IFN-γ Pathway Genes in Tumor Cells as a Mechanism of Resistance to Anti-CTLA-4 Therapy. Cell 2016;167:397-404.e9. [Crossref] [PubMed]
69. Ikeda H, Old LJ, Schreiber RD. The roles of IFN gamma in protection against tumor development and cancer immunoediting. Cytokine Growth Factor Rev 2002;13:95-109. [Crossref] [PubMed]
70. Hays E, Bonavida B. YY1 regulates cancer cell immune resistance by modulating PD-L1 expression. Drug Resist Updat 2019;43:10-28. [Crossref] [PubMed]
71. Lastwika KJ, Wilson W 3rd, Li QK, et al. Control of PD-L1 Expression by Oncogenic Activation of the AKT-mTOR Pathway in Non-Small Cell Lung Cancer. Cancer Res 2016;76:227-38. [Crossref] [PubMed]
72. Olbryt M, Rajczykowski M, Widłak W. Biological Factors behind Melanoma Response to Immune Checkpoint Inhibitors. Int J Mol Sci 2020;21:4071. [Crossref] [PubMed]
73. Zhang Z, Ni Y, Chen G, et al. Construction of immune-related risk signature for uveal melanoma. Artif Cells Nanomed Biotechnol 2020;48:912-9. [Crossref] [PubMed]
74. Mitra S, Lauss M, Cabrita R, et al. Analysis of DNA methylation patterns in the tumor immune microenvironment of metastatic melanoma. Mol Oncol 2020;14:933-50. [Crossref] [PubMed]

(English Language Editor: J. Gray)