The Cancer Genome Atlas (TCGA) skin cutaneous melanoma (SKCM) dataset, which included 98 primary and 356 metastatic melanoma samples with survival information, was retrieved from UCSC Xena (https://xenabrowser.net/datapages/). The GSE46517 (GPL96) dataset, which contained 31 samples of primary melanoma, 73 samples of metastatic melanoma, and seven control samples, was mined from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). TCGA-SKCM and GSE46517 datasets were used as training sets 1 and 2, respectively. The validation set GSE65904 (GPL10558), which contained survival information of 210 melanoma tumor samples, and the single-cell dataset GSE72056 (GPL18573), which contained gene expression data for 4,645 quality-controlled (QC) cells, were obtained from the GEO database. Gln metabolism-related genes (GRGs) were mined from the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb).

The differentially expressed genes (DEGs) between the melanoma and control groups in training set 2 were selected using the R language limma package (v 3.52.4) (15), with adj. p < 0.05 and |log₂FC| ≥ 0.5. Moreover, the enrichment analysis of Gene Ontology (GO) [including cellular component (CC), molecular function (MF), and biological process (BP) analyses] and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of the DEGs were completed using the clusterProfiler package (v 4.4.4) (16). The single-sample gene set enrichment analysis (ssGSEA) algorithm of the GSVA package (v 1.44.5) (17) was applied to compute the Gln metabolism score (GMS) in all samples of training set 1, with GRGs serving as the background gene set. Then, all genes of training set 1 were assigned to modules utilizing weighted gene coexpression network analysis (WGCNA) (v 1.71) (18). Modules relevant to GMS (p < 0.05) were confirmed as key modules, and key module genes were utilized for subsequent analyses.

Data from the single-cell dataset were integrated using the Seurat package for R (v 4.3.0) (19), with QC for the number of genes contained in the cells >1,700, housekeeping expression (corrected) >3, and all genes expressed >2% in at least five cells. The vst method was selected to screen the top 2,000 highly variable genes for downstream analysis. Subsequently, Uniform Manifold Approximation and Projection (UMAP) was utilized to reduce the dimensions. Using GRGs as the background gene set, the RcentageFeatureSet was utilized to calculate the percentage of GRG expression levels in each cell, and all cells were categorized into high- and low-expression groups to select intercellular differentially expressed GRGs (DE-GRGs) according to the median percentages.

Overlapping DEGs and key module genes were obtained from intersecting genes. Based on the intersecting genes, a protein–protein interaction (PPI) network was constructed. Then, the correlations between the intersecting genes and the DE-GRGs were computed according to Pearson’s correlation in training set 2, and the genes with |r| > 0.3, at least one differential GRG, and p < 0.05 were retained as candidate genes.

Based on the candidate genes via the glmnet package (v 4.1-4), key genes were selected via univariate Cox regression and least absolute shrinkage and selection operator (LASSO) regression analyses. Then, according to the correlation between expression of key genes and overall survival (OS), a risk model was constructed via LASSO (19). Risk scores were assessed utilizing the following formula:

where coef and X indicate coefficients and gene expression, respectively. Moreover, the samples of training set 1 and the validation set were sorted into the high- and low-risk groups based on the median risk score. The Kaplan–Meier (K–M) survival curves were drawn using the survival package (v 3.4-0) for both risk groups in training set 1 and the validation set (20). To further assure the validity of the risk model, receiver operating characteristic (ROC) curves were generated at 3, 5, and 7 years, and the area under the curve (AUC) values were computed using survivalROC (v 0.4) (21).

Risk scores and seven clinical characteristics (age, sex, Clark stage, metastasis status, Breslow status, and ulcer status) were entered into the risk model for univariate and multivariate Cox regression analyses. Then, independent prognostic factors were selected to construct a nomogram via rms (v 6.3-0) (22). The 3-, 5-, and 7-year survival rates were predicted depending on the total points (the higher the points, the lower the survival rate). The predictive ability of the nomogram was assessed using calibration curves. Correlations between the risk score and eight clinical characteristics were analyzed via correlation analysis.

To understand prognostic gene-related biological functions and signaling pathways, the correlations between prognostic genes and other genes were calculated and sequenced separately in training set 1. Based on the C2:KEGG gene set downloaded from the msigdbr package in R (v 7.5.1), the sequenced genes were enriched using the GSEA function in R (adj. p < 0.05).

The immune-related genes identified in the literature were used as background gene sets (23), and the samples in the training set were analyzed using ssGSEA to obtain enrichment scores for 28 immune cell types. Differences in enrichment scores for each immune cell between the melanoma and control groups were analyzed via the Wilcoxon test. The stromal score, immune score, and ESTIMATE score (summed over the first two) of the samples in training set 1 were estimated using the estimate package in R (v 1.0.13). Moreover, the expression of common immune checkpoints, including PD-L1, CTLA-4, LAG-3, GAL9, TIM-3, PD-1, PD-1LG2, and TIGIT, was compared between the high- and low-risk groups. Prognostic gene-related miRNAs were predicted using the starBase database.

The 50% inhibitory concentration (IC₅₀) values of 198 chemotherapeutic drugs were computed and compared according to the Genomics of Drug Sensitivity in Cancer database using the OncoPredict package (v 0.2).

The cell lines A375 (primary cutaneous melanoma) and A2058 (metastatic melanoma), purchased from ATCC (cat. nos. CRL-1619 and CRL-11147), were cultured in Dulbecco’s modified Eagle’s medium (DMEM; #12634-010, USA) supplemented with 10% fetal bovine serum (FBS; No. SH30070.02, HyClone, Utah, USA) in incubators with 5% CO₂ at 37°C. The human immortalized keratinocyte HaCaT cell line (CVCL-0038) as the control cell line was purchased from the Kunming Institute of Zoology and cultured in Dulbecco’s modified Eagle’s medium.

Total RNA was extracted from melanoma cell lines by TRIzol (15596018, Thermo, Beijing, China). A PrimeScript™ RT kit (R232-01, Vazyme, Nanjing, China) was applied to synthesize cDNA. Real-time polymerase chain reaction (RT-PCR) was achieved using SYBR Green Master Mix (Q111-02, Vazyme), and the expression levels were confirmed via the 2^−ΔΔCt method. The expression of each mRNA was standardized to the expression of GAPDH mRNA. All primers, as shown in Supplementary Table 1 , were purchased from Tsingke Biotech (Beijing, China).

A total of 3,216 DEGs (1,784 upregulated and 1,432 downregulated) between the melanoma and control groups were selected from training set 2 ( Figures 1A, B ). Functional enrichment analysis indicated that DEGs were related to 1,238 GO terms ( Figure 1C ), including 83 CCs (e.g., “collagen-containing”, “extracellular matrix”, and “cornified envelope membrane raft”), 56 MFs (e.g., “cytokine binding”, “cytokine activity”, and “actin binding”), and 1,099 BPs (e.g., “epidermis development”, “skin development”, and “epithelial cell proliferation”). In addition, 76 functional pathways were enriched according to the KEGG analysis. The KEGG enrichment analysis results displayed that chemical carcinogenesis-receptor activation, focal adhesion, and apoptosis were the pathways enriched in the DEGs ( Figure 1D ).

First, samples from the training set were used to construct a clustering tree ( Figure 2A ). A soft threshold of 7 (R² = 0.98) was applied to construct a scale-free network ( Figure 2B ). Then, an adjacency matrix and topological overlap matrix were constructed ( Figure 2C ). Finally, 13 modules were obtained based on average hierarchical clustering and dynamic tree clipping ( Figure 2D ). The blue module (|cor| = −0.34, p < 0.05, containing 2,165 module genes) was associated with GMS, which was identified as a key module ( Figure 2E ). Subsequently, 75 intersecting genes were obtained by overlapping DEGs and module genes ( Figure 2F ). A PPI network of 75 intersecting genes showed multiple pairs of relationships for intersecting genes. For example, CYP2E1 was associated with multiple genes, such as GPT and NR1I2 ( Figure 2G ).

After QC, 2,887 cells and 23,684 genes remained in the single-cell dataset ( Figure 3A ). Then, 2,000 highly variable genes were selected for subsequent analysis ( Figure 3B ). The cells were clustered into 13 clusters based on distance and were annotated to six cell types [T cells, B cells, Cancer-Associated Fibroblasts (CAFs), macrophages, natural killer (NK) cells, and endothelial cells] via marker genes ( Figures 3C, D ). Furthermore, 81 DEGs were screened between the high- and low-expression groups (min.pct = 0.25, logfc.threshold = 0.25) and were crossed with 80 GRGs to obtain 14 DE-GRGs (ALDH18A1, ASL, ASNSD1, ATP2B4, ALDH18A1, ASL, ASNSD1, ATP2B4, CLN3, FPGS, GLS, GLUD1, GMPS, GOT2, MTHFS, NIT2, OAT, and UCP2) ( Figures 3E, F ). In addition, 65 candidate genes were obtained based on Pearson’s correlation analysis of the intersecting genes and DE-GRGs ( Figure 3G ).

A total of nine genes were identified via univariate Cox regression analysis ( Figure 4A ), and further eight key genes (CHMP4A, IFFO1, ANKRD10, ZDHHC11, CLPB, ANKMY1, TCAP, and POLG2) were identified via LASSO based on 65 candidate genes ( Figure 4B ). Subsequently, a risk model was constructed according to the expression of eight key genes, and risk scores were also computed. Risk curves ( Figure 4C ) and gene expression data of the two risk groups were plotted based on risk scores ( Figure 4D ). It was observed from the K–M curves that the difference in the survival of melanoma patients was highly significant (p < 0.005) ( Figure 4E ). The AUC values exceeded 0.6 at 3, 5, and 7 years for melanoma patients. It suggested that the eight key genes could reliably predict survival status ( Figure 4F ).

To screen independent prognostic factors, clinical characteristics and risk scores were subjected to univariate and multivariate Cox analyses. The risk score, Clark stage, and ulcer status were identified as independent prognostic factors, which were used to construct a nomogram ( Figures 5A–C ). The slope of each calibration curve was close to 1, indicating favorable prediction accuracy of the nomogram ( Figure 5D ). In addition, correlation analysis of seven clinical characteristics and CD274 expression demonstrated that major differences existed in Clark stage and CD274 expression in melanomas (p = 0.019 and p < 0.001, respectively; Figure 6A ). There was a marked difference in the survival status of the samples in the Clark subgroups ( Figure 6B ).

GSEA demonstrated that eight key genes were enriched in KEGG pathways ( Figure 7 ).

There were significantly different enrichment scores of the 16 immune cells between the high- and low-risk groups (p < 0.05). The expression of immune cells in the high-risk group, except for CD56dim NK cells, was higher than that of the low-risk group ( Figure 8A ). A heatmap of the correlation between the eight key genes and risk scores of the 16 immune cells was also drawn in Figure 8B . The stromal, immune, and ESTIMATE scores were substantially different between the two risk groups (p < 0.05), and all were low in the high-risk group ( Figure 8C ). In addition, all immune checkpoint genes exhibited low expression in the high-risk group (p < 0.001) ( Figure 8D ).

Drug sensitivity analysis revealed that the high-risk group was strongly sensitive to rapamycin (p < 0.0001, Figures 9A, B ). A total of 81 miRNAs (e.g., hsa-miR-421, hsa-miR-449a, and hsa-miR-375) associated with the key genes were predicted using the starBase database ( Figure 9C ).

In GSE46517, IFFO1, ANKRD10, CLPB, TCAP, and POLG2 displayed high expression, but the expression levels of CHMP4A, ZDHHC11, and ANKMY1 were low in the melanoma group ( Figure 10A ). The expression of these genes was further examined via RT-PCR in HACAT and melanoma cells (A375 and A2058). It was interesting that IFFO1, ANKRD10, and POLG2 were markedly upregulated and that CHMP4A was also markedly downregulated in A375 cells ( Figure 10B ), which was partially consistent with the results of biological analysis.