Overall, 937 patients from six suitable skin cutaneous melanoma (SKCM) cohorts [GSE65904, GSE59455, GSE19234, GSE78220, GSE91061, GSE54467, and The Cancer Genome Atlas (TCGA)-SKCM] were subjected to analyses. Patients without survival information and RNA-seq data were excluded. For the dataset from Gene Expression Omnibus (GEO), we downloaded preprocessed clinical and transcriptome data using the R ‘GEOquery’ package (19) and merged the related GEO dataset using the ‘ComBat’ algorithm from the R ‘sva’ package (20). As for the dataset from TCGA, we used the R ‘TCGAbiolinks’ package (21) to download all available transcriptome [fragments per kilobase of transcript per million fragments mapped (FPKM) value] and clinical data from the Genomic Data Commons (GDC, https://portal.gdc.cancer.gov/). Furthermore, clinical information pertaining to ulceration in TCGA patients was extracted from the UCSC Cancer Browser (XENA, https://xenabrowser.net/datapages/). FPKM values were then converted to transcripts per million (TPM) values, which were used in subsequent analyses.

We used the R ‘TCGAbiolinks’ package (21) to download raw read counts for TCGA-SKCM patients from the GDC. The R ‘DESeq2’ package (22) was applied to standardize the read counts and perform differential gene analysis. And for the multiclass DEseq2, we only took the genes in the top 25% of variance for further differential gene expression analysis. Ulcer-immunity-related DEGs between the “ulcer_low-immunity” and “nonulcer_high-immunity” groups were identified as those with fold change > 2 and false discovery rate (FDR) < 0.05. FDR < 0.05 was the significance criterion for cluster-specific genes.

ssGSEA and GSVA were performed using the R ‘gsva’ package (23). The gene set ‘c2.cp.kegg.v7.4.symbols.gmt’ was downloaded from MSigDB v7.4 for GSVA. For ssGSEA, 24 immune cell signatures (24) were used to describe immune cell populations in individual patients.

To identify prognostic genes with P < 0.05, univariate Cox regression analysis was performed using DEGs identified on comparing the “ulcer_low-immunity” and “nonulcer_high-immunity” groups. Based on 53 prognostic genes, we used the R ‘ConsensuClusterPlus’ package (25) to generate robust clusters of TCGA-SKCM patients. The delta area and cumulative distribution function were used to identify the optimal K value of 3. To validate the subtype in the three merged GEO cohorts (GSE65904, GSE59455, GSE19234), The Nearest Template Prediction (NTP) method provides a convenient way to make classification predictions using only a list of signature genes and a test data set for predictive confidence assessment of gene expression data for each individual patient (26). In our study, the top 200 upregulated of each cluster in TCGA (600 genes in total) were applied to predict the clusters in merged GEO dataset by applying the R ‘MOVICs’ package (27).

Somatic mutations and somatic copy number alternations (CNAs) were downloaded from the GDC using the R ‘TCGAbiolinks’ package (21). The R ‘Maftools’ package (28) was then used to visualize and analyze data pertaining to somatic mutations and CNAs (GISTIC output). GISTIC 2.0 (29) was used to identify significant copy number amplifications and deletions. Chi-square and Fisher’s exact test were used to detect differential mutated genes and differentially copy number gains and losses.

The EPIC (30), CIBERSORT (31, 32), MCPcounter (33), and Quantiseq (34) algorithms were used to estimate the fraction of immune cells based on the transcriptome data (TPM value) of TCGA-SKCM patients.

The PPIs between differentially expressed genes were searched using the Search Interacting Genes/Proteins (STRING) v.11.5 database with a confidence level of 0.7, and the interaction network was visualized using Cytoscape software. In a PPI network, nodes with higher connectivity are more important for maintaining the stability of the entire network and are therefore considered more relevant to the overall biological process. The radiality method is a topological analysis method that plays an important role in predicting key nodes (35). Using the radiality method of CytoHubba plugin (35) in Cytoscape software, the top 10 genes in the network were screened as hub genes by radiality method.

We herein used the A375 (ATCC, Cellcook Biotechnology, Guangzhou, China) and SK-MEL-28 (ATCC, Cellcook Biotechnology) cell lines. A375 and SK-MEL-28 cells were cultured in DMEM (Gibco, Carlsbad, CA, USA) supplemented with 10% fetal bovine serum (Gibco) and Penicillin Steptomycin (Invitrogen, Carlsbad, CA, USA). For siRNA transfection, A375 and SK-MEL-28 cells were transfected with EIF3B siRNA (RiboBio, Guangzhou, China) using Lipofectamine 3000 (Invitrogen), according to manufacturer instructions. The sequence of siRNAs for targeting EIF3B was as follows: CCTTAGCGTTTGTGGACACTT (EIF3B siRNA1) and CGGGAAGATTGAACTCATCAA (EIF3B siRNA2). Total RNA extraction, purification, and qRT-PCR were performed as previously reported. The sequences of primers used for qRT-PCR were as follows: (1) β-actin: forward primer, CTCGCCTTTGCCGATCC and reverse primer, TTCTCCATGTCGTCCCAGTT and (2) EIF3B: forward primer, CTGGGTGCCTGAAGACAAAGA and reverse primer, CGTTCTTCTGCCAATGGAGC.

Western blotting was performed as previously described (36) using primary antibodies against α-tubulin (A11126, Invitrogen) and EIF3B (MA5-36159, Thermo Fisher).

For the transwell migration assay, after 16h serum-free starvation, the cells were resuspended in serum-free media; subsequently, 250000 cells in 300 µL serum-free media were seeded into transwell inserts (Corning) with 8-µm pore size. Different treatment media were then added in the lower chamber. For the invasion assay, the inserts were coated with 40 µL matrigel (1 mg/ml, BD Biosciences). The cells were then seeded onto the coated inserts and incubated with different treatment media. After 24 h or 48 h of incubation at 37°C in a CO₂ incubator, the cells and media were carefully removed from the top of the insert, and the migration or invasion inserts were placed into a clean well. The cells were then fixed in 4% paraformaldehyde and stained with 0.2% crystal violet. After wiping out upper cells in the insert, the cells which grew through the porous membrane were photographed by an inverted light microscope (×100). The relative numbers of migrating and invasive cells were counted by using ImageJ software.

To demonstrate the efficacy of immunotherapy, we integrated two datasets (GSE78220 and GSE91061) of patients with melanoma who had received anti-PD-1 therapy, and only therapy-naïve patients further analyzed. The submap tool (29) was used to predict responsiveness to immunotherapy depending on EIF3B expression.

Cells were firstly seed at 1×10⁵/well in 12-well plate and the supernatant was collected at 48h. The amount of TGF-β1 in supernatant of EIF3B knockdown cells and scramble control cells were determined by ELISA specific for human TGF-β1(R&D Systems Inc., Minneapolis, MN, USA). The assay was performed following the manufacturer’s instruction. Optical density at 450 nm was measured and the final concentration of TGF-β1 (pg/ml) was calculated according to standard protein curve.

Firstly, Multiclass DESeq2 was used to identify 1744 cluster-specific altered genes with an FDR-adjusted P value of < 0.01. Secondly, Univariate cox analysis was performed for the 1744 genes and identify 709 genes with prognostic value in the TCGA SKCM dataset. Thirdly, LASSO (least absolute shrinkage and selection operator) logistic regression with 10-fold cross-validation was used to further reduce candidate genes using the ‘glmnet’ R package (37). Then, multivariate cox analysis (38) was used to further screen genes using the Cox multivariate proportional hazard regression model with a stepwise method (both). The risk score was calculated using the following formula: Risk score = 0.278×EI24 + 0.166×HEYL – 0.09×IFIT3 – 0.18×SNTB1 – 0.115×CSF1R. The risk score for each patient in the TCGA data set and the validation data set was calculated according to this formula. Patients were classified into high-risk and low-risk groups based on the median cut-off for risk scores. Survival analysis was performed using the ‘survival’ R package to assess differences in OS between high- and low-risk groups. To measure the specificity and sensitivity of the prognostic capability of this model, we calculated the area under the curve (AUC) using the R ‘timeRoc’ package (39).

Kaplan–Meier curve and log-rank test were used to evaluate differences in survival rate between the groups. Univariate and multivariate Cox regression analyses were used to determine prognostic factors. Pearson and Spearman correlation analyses were performed for calculating the correlation coefficient. The unpaired Student’s t-test and Mann–Whitney U-test were used for normally and non-normally distributed variables, respectively, so as to compare the two groups. To compare more than two groups, one-way analysis of variance (ANOVA) and Kruskal–Wallis test were applied as parametric and non-parametric methods, respectively. Chi-square and Fisher’s exact tests were used to examine differentially mutated genes and differential copy number gains and losses. The R ‘Maftools’ package (28) aided OncoPrint generation. The R ‘ggplot2’ (40) and ‘ComplexHeatmap’ packages (41) were used for data visualization. All survival curves were generated using the R ‘survival’ (42) and ‘survminer’ (43) packages. Statistical analysis was performed using the R software, and values represent mean ± standard deviation (SD). P < 0.05 indicated statistical significance.

First, we explored TCGA-SKCM dataset; Table S1 details pertaining to the clinical characteristics of TCGA cohort. In case of TCGA database, patients with ulceration had worse prognosis than those with non-ulceration ( Figure 1A ). Tumor samples consist of cancer and many non-tumor cells, such as infiltrating immune and stromal cells, as well as other non-cellular components, which crosstalk with each other and ultimately promote cancer progression (44). For TCGA-SKCM patients, immune cell infiltration was measured via ssGSEA using previously published immune cell signatures. The SKCM samples in TCGA were divided into two groups (high and low immunity) depending on the degree of immune cell infiltration ( Supplementary Figure 1 ). As with previous reports, the high immune cell infiltration group was associated with better overall survival (OS) than that the low immune cell infiltration group ( Figure 1B ). Based on our clinical-related data, patients were further divided into three groups: “ulcer_low-immunity,” “nonulcer_high-immunity,” and “others.” Survival analyses showed significant differences between the “ulcer_low-immunity” and “nonulcer_high-immunity” groups (P = 0.0014; Figure 1C ). Further, DESeq2 was used to identify ulcer-immunity-related DEGs on comparing the “ulcer_low-immunity” and “nonulcer_high-immunity” groups. Overall, 158 genes with a Benjamini and Hochberg-corrected P value < 0.05 and fold change > 2 were identified as primary DEGs ( Figure 1D ). We used these genes and OS data available in TCGA-SKCM dataset to perform a survival analysis using the univariate Cox proportional hazards model, and 53 genes were identified as ulcer-immunity-related DEGs. We used the R ‘ConsensusClusterPlus’ package for consistent clustering of genes in TCGA-SKCM dataset, and category identification based on the 53 ulcer-immunity-related DEGs was performed. Using the unsupervised clustering method, three different clusters were finally identified: 177 cases in Cluster 1, 200 in Cluster 2, and 90 in Cluster 3 ( Figure 1E and Supplementary Figures 2A, B ). The survival analysis revealed that Cluster 3 had the worst prognosis in comparison to the other clusters ( Figure 1F ). Previous study had shown that TCGA-SKCM sample can be further divided into three molecular subtypes, including ‘immune’, ‘keratin’ and ‘MITF-low’ (45). Furthermore, Alexander Bagaev et al. establish a classifier based on cancer microenvironment and classify tumor samples into four subtypes termed as (1) immune-enriched, fibrotic (IE/F); (2) immune-enriched, non-fibrotic (IE); (3) fibrotic (F); and (4) immune-depleted (D) (46). Interestingly, we found that the proportion of patients with low immune infiltration is highest in cluster 3 (51.2%) and most of the samples in cluster3 were defined as ulceration (78.6%). Moreover, most of the samples in cluster3 were defined as the keratin subtype (79.1%) and immune-depleted subtype (47.7%, Supplementary Figures 2E–H ). Multiclass DESeq2 was used to identify cluster-specific upregulated genes with an FDR-adjusted P value of < 0.05. Three GEO datasets with usable OS data and clinical information (GSE65904, GSE19234, and GSE59455; Table S1 ) were combined into one metacohort. The three clusters were validated in the merged GEO cohort using the NTP algorithm (FDR < 0.05). Cluster 3 was still associated with the worst prognosis ( Supplementary Figures 2C, D ).

To uncover the activation of the signaling pathways in each subtype, we performed GSVA and calculated the enrichment score for KEGG signaling pathways with MSigDB v7.4. We then selected a total of 30 most representative Cluster 1–3 gene sets and created a heatmap showing specific gene sets for each subtype ( Figure 2A ). As evident from the heatmap, in comparison with Clusters 1 and 2, Cluster 3 was characterized by the lack of signals related to the immune system and immune cells. To analyze mutations in each subtype, we created a waterfall chart to depict the top 20 significantly mutated genes (SMGs) in the three clusters ( Figure 2C and Supplementary Figures 3A, B ). This chart showed that 50% SMGs (TTN, MUC16, BRAF, PCLO, DNAH5, DNAH7, ADGRV1, LRP1B, ANK3, and CSMD1) were shared by Clusters 1–3. Clusters 1 and 3, but not Cluster 2, shared the top three SMGs: FLG (35% and 28%, respectively), XIRP2 (34% and 34%, respectively), and CSMD3 (33% and 26%, respectively). Clusters 2 and 3, but not Cluster 1, also shared the top three SMGs: MGAM (38% and 30%, respectively), HYDIN (34% and 30%, respectively), and USH2A (34% and 32%, respectively). Moreover, Cluster 3 showed four unique top SMGs: FAT4 (33%), FAM135B (28%), ZFHX4 (26%), and MUC17 (24%). In a recent prospective study, higher tumor mutation burden (TMB) was found to be associated with a better response to immunotherapy (47). The tumor mutations in Cluster 2 samples were slightly higher than those in Cluster 1 samples, and the severity of TMB in Cluster 1 and 2 samples was higher than that in Cluster 3 samples ( Figure 2B ). We then used GISTIC to analyze data related to somatic CNAs to identify areas that were repeatedly amplified and deleted among Clusters 1–3, and found an obvious similarity among chromosomal aberrations in the three clusters ( Figure 3A and Supplementary Figure 4A ). The similarity between Clusters 2 and 3 was stronger than that between Clusters 1 and 3. In total, 29 focal deletion and 17 focal amplification peaks were detected in Cluster 3, while 41 and 44 focal loss and 43 and 49 focal gain peaks were detected in Clusters 1 and 2, respectively ( Figure 3B and Supplementary Figures 4B, C ). Next, we examined the frequency of amplification or deletion for immune checkpoint genes in each subtype. Many immune checkpoint genes were deleted in Cluster 3 (PD-1, GZMA, CD276, CCL5, and VTCN1), while several were amplified in Cluster 1 (CD4, CD274, CD276, and HLA gene family) and Cluster 2 (CD274 and PDCD1LG2; Figure 3B and Supplementary Figures 4B, C ).

Our data indicated that Cluster 3 was negatively associated with immune-related pathways and mainly immune-depleted subtype. We investigated the heterogeneity of immune cell infiltration in the three clusters based on the EPIC, MCPcounter, Quantiseq, and CIBERSORT algorithms. Tumor-infiltrating immune cells in 467 TCGA melanoma patients was showed in the heatmap ( Figure 4 ). B cells infiltration in Cluster 3 showed a significantly decrease according to all the four algorithms; further, according to three algorithms, macrophages and CD8+ T cells infiltration in Cluster 3 was significantly lower. Collectively, these results revealed that Cluster 3 had a different immune phenotype than Clusters 1 and 2, with Cluster 3 showing less immune cell infiltration and less immune activation.

We investigated the expression of HLA family genes and immune checkpoint markers among the three subtypes in TCGA and GEO datasets. In case of HLA family genes, the expression levels of HLA-DQA2, HLA-DOB, HLA-DOA, and HLA-DMB in Cluster 3 of TCGA dataset were less than the other two clusters; there was no significant difference (ANOVA test, P < 0.05) in the expression level of other HLA family genes among the three clusters ( Figure 5A ). However, in the GEO dataset (GSE65904 + GSE19234 + GSE59455), all HLA family genes, except HLA-G, HLA-DQB2, exhibited downregulated expression levels (ANOVA, P < 0.05) in Cluster 3 relative to the other two clusters ( Figure 5C ). We then assessed immune checkpoint markers associated with antigen presentation, cell surface receptors, co-inhibition, ligands, and cell adhesion. In comparison with Clusters 1 and 2, Cluster 3 in TCGA dataset showed reduced expression of all immune checkpoint markers, except GZMB, TNF, CCL5, and CXCR3. Besides, in the GEO dataset, Cluster 3 showed reduced expression of all immune checkpoint markers ( Figures 5B, D ). Furthermore, we combined GSE78220 and GSE91061 with available response to immunotherapy and used the specific gene of three clusters for analyses with the NTP algorithm to predict the subtypes in the two GEO datasets (FDR < 0.05; Supplementary Figure 5C ). We also grouped the treatment response into a binary model and found that the percentage of patients with stable/progressive disease group in Cluster 3 was higher than the other two clusters ( Figure 5E ).

In comparison with Clusters 1 and 2, Cluster 3 showed the worst prognosis and responsiveness to immunotherapy. According to our results, Cluster 3 can be used as an independent prognostic factor in both TCGA and GEO databases ( Table S2 ). We then aimed to identify a hub gene in Cluster 3 with a key role in melanoma progression. To discover the main genes responsible for SKCM growth, CRISPR-based genome-wide loss-of-function screening was employed (DepMap, https://depmap.org/portal/download/), and 648 genes in the SKCM cell line were found to be important for survival. Among these 648 candidate genes, 60 were specific to Cluster 3 ( Figure 6A ). Supplementary Figures 5A, B shows the expression levels of these 60 genes in all three clusters. We then explored the PPI network based on these 60 Cluster 3-specific genes and the top 10 genes in the network were screened as hub genes by radiality method (35). ( Supplementary Figure 5D ). Survival analysis revealed that among these top 10 genes, only EIF3B had prognostic significance both in TCGA and GEO databases; higher expression level of EIF3B was associated with worse prognosis in patients with melanoma ( Figures 6B, C ). We also found that Cluster 3 showed the highest expression of EIF3B ( Figures 6D, E ). EIF3B has been previously reported to play a vital role in the progression of several types of cancers (17, 18, 48), but its effectiveness in melanoma remains to be reported. To further study the effects of EIF3B on the biofunction of melanoma cells in vitro, we constructed stable knockdown A375 and SK-MEL-28 cell lines that expressed EIF3B-specific siRNAs ( Figure 6F ). As showed in Figures 7A–D , the results of transwell migration assay and invasion assay showed that the relative number ratios of migrating and invasive cells were significantly decreased in EIF3B knockdown cells compared with si-NC cells. It indicated that EIF3B knockdown could suppress migration and invasion abilities in A375 and SK-MEL-28 cells.

To explore EIF3B-associated biological functions in melanoma, we performed GSVA using TCGA-SKCM dataset. High expression of EIF3B was associated with significantly downregulated immune-related pathways, such as the T cell receptor pathway, B cell receptor pathway, antigen processing and presentation, and leukocyte transendothelial migration ( Figure 8A ). Based on these finding, we believe that EIF3B plays a key role in the immune response in melanoma. We used the Spearman order correlation method to detect the correlation between immune system cells and EIF3B in TCGA dataset. EIF3B expression was found to be negatively correlated with immune cell levels, including those of effector cells, in immunotherapy ( Figure 8B ). Furthermore, we investigated the correlation between seven immune checkpoint markers (LAG3, PD-1, CD8A, GZMB, CTLA-4, BTLA, and IFNG) and EIF3B expression in TCGA and GEO datasets, and found that EIF3B expression was negatively correlated with all of them in both datasets ( Figure 8C and Supplementary Figure 6A ). These results suggested that high expression of EIF3B leads to deficient proinflammatory immune cell infiltration and might predict a worse response to immunotherapy. We then performed submap analyses to evaluate the response of melanoma patients with high and low EIF3B expression to anti-PD-1 immunotherapy. Interestingly, patients with low EIF3B expression showed partial response to anti-PD-1 immunotherapy ( Figure 8D ). The effectiveness of cancer immunotherapy is highly dependent on the development and activation of the steps in the cancer immunity cycle (49). Therefore, we further analyzed the correlation between EIF3B expression and anti-cancer immunity cycle and found EIF3B expression was negatively correlated with all the steps in the cycle ( Supplementary Figure 6B ). Furthermore, the ELISA results showed that the concentration of TGF-β1 in supernatant was significantly decreased in EIF3B knockdown cells compared with si-NC cells ( Figures 7E, F ). These results suggested that EIF3B plays a suppressive role in melanoma immune response and immunotherapy.

We firstly performed Lasso logistic regression and stepwise multivariate cox analysis to constructed a prognostic model based on cluster specific altered genes with prognostic value in the TCGA SKCM dataset ( Supplementary Figures 7A, B ). The risk score was calculated as follows: Risk score = 0.278×EI24 + 0.166×HEYL - 0.09×IFIT3 - 0.18× SNTB1 - 0.115×CSF1R. Then, patients in TCGA dataset were divided in high- and low-risk groups based on their risk score. We found that the expression of EIF3B was higher in high-risk groups ( Supplementary Figure 8A ) and the patients in cluster3 had a higher risk score ( Supplementary Figures 8B, C ). Kaplan-Meier survival analysis was performed on the training and testing datasets to assess the predictive power of our ulcer-immunity related signatures. In the TCGA dataset, high-risk patients exhibited worse OS compared to low-risk patients (P < 0.001, Figure 9A ). A similar result was observed in the testing datasets (GSE65904, GSE54467, GSE59455) ( Figures 9B, C and Supplementary Figure 8D ). Meantime, time-dependent AUC and AUC at 1, 2, 3, and 5 years suggested that the ulcer-immunity related score had a significant value in predicting the OS of melanoma patients in the TCGA and GEO datasets ( Figures 9D–F and Supplementary Figures 8E–I ). To examine the combined prognostic value of all training and testing datasets, we performed a prognostic meta-analysis. Results showed that the ulcer-immunity related score was a significant risk factor for overall survival of patients with melanoma (combined HR = 2.08, 95% CI = 1.73-2.51, P < 0.0001, Figure 9G ). Furthermore, we calculated the risk score of the patients in GSE78220 and GSE91061 datasets with available response to immunotherapy, and it was found that the percentage of stable/progressive disease in high-risk patients was higher than that in low-risk patients ( Figure 9H ).