The mRNA-seq datasets and clinical information of CM patients and normal controls were acquired from the TCGA-SKCM project and the GTEx website. The key enzymes involved in regulating glutamine metabolism were summarized from one previous study (Altman et al., 2016). We defined the genes encoding these 22 proteins as GMRGs. The single-cell count matrix was obtained from the GSE72056 (Tirosh et al., 2016), GSE115978 (Jerby-Arnon et al., 2018), and GSE120575 (Sade-Feldman et al., 2018) datasets from the GEO database. The gene expression matrix and clinical information used to validate the expression of the key gene and its correlation with clinicopathologic parameters were also obtained from the GSE3189 (Talantov et al., 2005), GSE15605 (Raskin et al., 2013), GSE114445 (Yan et al., 2019), and GSE46517 (Kabbarah et al., 2010) datasets from the GEO database. The online Xiantao tool (https://www.xiantao.love/products) is a comprehensive bioinformatics analysis website. Differential analysis, survival analysis, biofunctional enrichment analysis, and immune-related analyses of target genes were performed by this tool partly. The predicted 3D structure and immunohistochemical staining of the protein encoded by the key gene were obtained from the HPA database. The main R packages used in the analysis process were displayed in Supplementary Table S2.

The prognostic GMRGs correlated significantly with the OS of CM patients were screened out by univariate Cox analysis. The threshold of screening was p < 0.05. The composition of the optimal multi-gene signature was determined by successively applying the LASSO and multivariate Cox regression analyses (Gui and Li, 2005). The standardized expression value of each signature gene is multiplied by the corresponding regression coefficient and the results are summed to obtain the risk score of this signature, which is defined as GMRS. CM patients were stratified into two subgroups according to the median GMRSs. Survival analysis was performed to evaluate survival differences. The Time-dependent ROC curves were used to evaluate the predictive performance of our signature. The predictive independence of classical prognostic factors and our signature was identified by regression analyses.

To confirm the consistency between our gene signature and conventional predictors, we performed a hierarchical difference analysis. Similarly, we also discussed the distribution of clinical parameters between GMRS subgroups. Moreover, the subgroup survival analysis was conducted based on different clinical parameters to verify the predictive ability of our model in different clinical signs.

Eight algorithms were used to compare the immune infiltration status between GMRS subgroups, including the ssGSEA algorithm and the corresponding gene sets used (Charoentong et al., 2017). The ESTIMATE algorithm allowed for a quantitative comparison of immune correlation status among GMRS subgroups using 4 kinds of scores (Yoshihara et al., 2013). The TIP tool was used to evaluate 23 standardized immune activity scores between GMRS subgroups during the cancer immune cycle (Xu et al., 2018).

The molecular subtypes of CM samples were identified by the consistent clustering method (Wilkerson and Hayes, 2010). The clustering algorithm is set to K-means and the coefficients k are set from 2 to 9. The optimal clustering result was presented as a heat map including the expression of modeled genes and the distribution of classical clinical parameters across cases. The distribution of different subtypes was explored by the t-SNE method (Laurens and Hinton, 2008).

The Wilcoxon test with an adjusted p < 0.05 and a |log2 FC| > 1 threshold was used to screen out DEGs between clusters. The pathways enriched in each cluster were evaluated by GSVA. The GO and KEGG pathway enrichment analyses based on the above DEGs were also performed. The correlation between the expression of our key gene GOT2 and immune-related pathways from the KEGG database and Reactome database was evaluated by GSEA. The top 20 mutated genes were evaluated and the oncoplots were generated to provide a visual representation of the mutational landscape in GMRS subgroups. The survival analysis of the subgroups stratified by the TMB state of each sample was conducted.

The raw data for single-cell analysis first underwent a process of normalization, dimensionality reduction, and clustering (Stuart et al., 2019). The criteria for cell filtering were more than 200 gene expressions and less than 20% mitochondrial gene expression. The criteria for gene filtering were expressed in more than 3 cells. The top 30 PCs for UMAP construction and a resolution of 0.4 for graph-based clustering were utilized to identify each cell cluster. The InferCNV method was utilized to estimate the CNV signal for individual cells (Patel et al., 2014). To define the reference cell-inferred copy number profiles, B cells, fibroblast, and T cells were utilized. Epithelial cells were used for the observations. According to the flow described by Sunny Z. Wu et al., (Wu et al., 2021), Epithelial cells were classified as either normal epithelial or malignant. The potential communication molecules between various cell subtypes were investigated using the CellPhone DB database (Vento-Tormo et al., 2018). Following that, we calculated the mean and the significance of cell communication and visualize the network plot according to the interaction matrix and the cell count matrix.

Human melanoma cell lines A375 and SK-28 were purchased from the University of Colorado Cancer Center Cell Bank and cultured in DMEM medium supplemented with 10% FBS (Invitrogen, Carlsbad, CA, USA) at 37°C in a 5% CO₂ atmosphere. siRNAs were used to knock down the expression level of GOT2 in cell lines. The siRNA sequences were as follows: siRNA#1-GCTACAAGGTTATCGGTATTA; siRNA#2-GCCTTCACTATGGTCTGCAAA. We used RT-qPCR to determine the knockdown efficiency of siRNA-GOT2. CCK8 and clone formation were applied to detect the proliferation activity of melanoma cells.

All analyses were performed using R software (version 4.2.2). All statistical tests were two-sided. p-value <0.05 or Spearman correlation coefficient >0.3 was considered statistically significant unless otherwise noted.

The prognostic analysis was first carried out to filtrate GMRGs significantly related to outcomes in CM patients. As shown in the forest plot (Figure 1A), 9 GMRGs were found to be associated with prognosis (p < 0.05). 5 of the above 9 genes were regarded as protective genes and 6 were considered as risk genes. Then the multivariate Cox and LASSO regression analyses were applied in succession to single out the best model with the TCGA dataset (Figures 1B, C). Furthermore, A glutamine-metabolism-related prognostic signature for CM was constructed given the best coefficient value (Figure 1D). The risk score for our signature, which was defined as GMRS, was worked out through the established formula method: G M R S = − 0.3806 * exp D G L U C Y ⁡ + − 0.1598 * exp G L U D 2 ⁡ + 0.4241 * exp G O T 2 ⁡ + − 0.4854 * exp L A T ⁡ + 0.2668 * exp S L C 6 A 14

The GMRS was then calculated for each case. Its mean threshold divided the CM patients into two subgroups (low-GMRS (n = 229) and high-GMRS (n = 228)). KM analysis confirmed a lower likelihood of survival in the high-GMRS group (p < 0.05, Figure 1E). An evaluation of the predictive performance was conducted with the time-dependent ROC curve of the model and in Figure 1F, the AUC reached 0.696. Univariate and multivariate Cox analyses were also performed to further explore the independent predictive ability of our signature. The univariate Cox regression analysis showed that Breslow depth, Clark level, age, pathologic T stage, pathologic N stage, and our GMRS were all significantly associated with survival (Figure 1G). However, when applying the corresponding multivariate Cox regression analysis, we found that only age (p = 0.007), pathologic T stage (p = 0.025), pathologic N stage (P < 0.001), and our GMRS (p = 0.001) were clinically independent (Figure 1H).

In order to evaluate the association between GMRS and clinicopathologic factors, we performed differential analysis and evaluated the distribution characteristics of classic clinicopathologic factors between two risk groups. The bar plots in Supplementary Figure S1A showed that patients with advanced T stage (p < 0.001), overall stage (p < 0.05), and Clark level significantly presented higher GMRSs. Consistent with this result, the high-GMRS group exhibited notably higher proportions of advanced T stage (p = 0.002), N stage (p = 0.043), overall stage (p = 0.004), and Clark level (p = 0.007) (Supplementary Figures S1A, B). Further, we performed a stratification analysis for survival probability to determine the prognostic power of our prognostic signature in subgroups of CM. As was shown in Supplementary Figure S1C, except for the subgroup of Clark level I/II (p = 0.088), patients in the high-GMRS group showed a more abysmal outcome than those in the low-GMRS group (all p < 0.05).

Since a large variety of genetic alterations occur in the common progression trajectories of CM, the gene mutation status was comprehensively assayed in patients with different GMRSs (Shain and Bastian, 2016). The top 20 most mutative genes in each group were displayed (Figures 2A, B). Some of them overlapped, with a higher mutation frequency in the low-GMRS group. Figures 2C–E showed the differential expression of the model genes among the subgroups stratified by the three top gene mutations. Previous evidence revealed that the TMB severely affected the response to cancer immunotherapy and the long-term prognosis of tumor patients (Kang et al., 2020; Moldoveanu et al., 2022). According to our analysis, patients with high TMB experienced improved OS, relative to low TMB (p = 0.023, Figure 2F). Next, we evaluated the collaborative impact of our signature and TMB in predicting prognosis. The results showed that the high-GMRS group had worse OS than the low-GMRS group regardless of the TMB status, which implied that the prognostic value of our signature was not intervened by the TMB status of individuals (p < 0.001, Figure 2G).

Considering CM as immune-related, and immunotherapy as preferable (Coit et al., 2013), we explored the difference in the immune infiltration situation between the two risk subgroups. The heatmap in Figure 3A showed visually some immune cells demonstrated significant differences according to most of the algorithms, such as B cells, CD8⁺ T cells, neutrophils, macrophages, and NK cells. The specific correlation analysis was shown in Figure 3B. And we can conclude that most immune cells have a significant negative correlation with the value of GMRS (r < −0.3).

Whereafter, the enrichments of 23 immune cell types and 13 immune-related functions were assessed between the low-GMRS and high-GMRS groups. The significant between-group difference showed up in most of the active immune cells and immunocompetent (p < 0.05), with a higher-level immune cell infiltration in the low-GMRS group (Figure 3C). Both groups demonstrated significant differences in all immune-related functions, such as checkpoint, inflammation-promoting, Type I IFN response, etc. (Figure 3D). By the ESTIMATE algorithm, we explored differences in immune-related scores between GMRS groups. Lower immune, stromal, and ESTIMATE scores were observed in the high-GMRS group, while its tumor purity scores were significantly higher (all p < 0.05, Figure 3E).

Chen et al. referred to seven sequential processes of antitumor immunity as the “cancer-immunity cycle” (Chen and Mellman, 2013). TIP (a web service for determining tumor immunophenotype profiling) was utilized to evaluate the anticancer immunological functions of the seven-step cancer-immunity cycle between GMRS groups (Xu et al., 2018). Our results in Figure 3F revealed that patients in low-GMRS groups owned significant functional enhancement in step 1 (release of cancer cell antigens), step 3 (priming and activation of effector T cell responses), step 4 (trafficking of immune cells to tumors), and step 5 (infiltration of immune cells into tumors). Specifically in step 4, the low-GMRS group presented the high activity of recruiting of 10 main immune cells, including T cells, macrophages, NK cells, Th17 cells, and so on. Only the process of MDSC recruiting was enhanced in high-GMRS groups. Taken together, the patients in the high-GMRS group were predisposed to an immune-silent microenvironment featuring full down-modulation of immune cell infiltration and immune functions, which may be followed by a dismal prognosis.

According to the expression levels of five GMRGs from our prognostic signature, CM samples from TCGA were divided by clustering analysis. When k = 2, the consensus matrix heatmap in Supplementary Figure S2A exhibited high consistency between clusters. And the K-M survival curve in Supplementary Figure S2B showed a significant difference in survival between clusters (p = 0.015). The area under the CDF curve tended to be stable starting from k = 2 (Supplementary Figures S3A, B). The PCA plot in Supplementary Figure S2C also agreed with the clustering outcome of the two subtypes. The alluvial diagram demonstrated the distribution overlap of the CM samples between GMRS and molecular subtypes (Supplementary Figure S3C). The distribution of cluster A subtype samples was mainly comprised of high GMRSs, while cluster B subtype achieved a high ratio of the low-GMRS group in the samples. As shown in Supplementary Figure S2D, the heatmap illustrated the clinicopathologic features and expression patterns of five GMRGs.

To explore the between-group differences in biological features, GSVA was carried out to seek the enriched pathways in each cluster. Results with adjusted p < 0.05 were shown in Supplementary Figure S2E. Cluster A subtype was enriched in more tumor-related pathways such as “MYC_TARGETS_V2” (adj.p = 0.014), “P53_PATHWAY” (adj.p = 0.006), “WNT_BETA_CATENIN_SIGNALING” (adj.p = 0.031) and so on, which may account for the poor prognosis of patients in cluster A. More details can be obtained from Supplementary Table S3. Then, the enrichment analyses of the GO functions and KEGG pathways were performed. Enriched BPs of the DEGs mainly concluded “negative regulation of protein phosphorylation” and “cell-substrate adhesion”. Enriched of CCs were mainly in “cell-cell junction” and “collagen-containing extracellular matrix”. Enriched MFs consisted of “integrin binding” and “insulin-like growth factor I binding” (Supplementary Figure S2F). Enriched KEGG pathways were listed as “PI3K-Akt signaling pathway” and “Human papillomavirus infection” (Supplementary Figure S2G). These significantly enriched GO terms and KEGG pathways can offer a better understanding of the roles of these DEGs in the initiation and progression of CM.

The differences in prognostic value and immune characteristics between the two clusters suggested that specific molecular biomarkers may contribute to CM patients’ prognostication. To find these key genes, we conducted the differential analysis and plotted the ROC curve of diagnostic efficiency. We found differential expressions of all five model genes between tumor and normal tissues. DGLUCY and SLC6A14 were lowly expressed, while those of LAT, GOT2, and GLUD2 were highly expressed in CM tissues (all p < 0.001, Figure 4A). The AUCs of DGLUCY, SLC6A14, GOT2, and GLUD2 were all greater than 0.7 (Figure 4B), implying that every one of them owned so adequate diagnostic efficiency to become an independent diagnostic biomarker. After a validation and survival analysis of these five genes in the TCGA cohort (Figures 4C–G), only GOT2 with a high expression was found to be associated with poor prognosis in CM patients (p = 0.015, Figure 4G). We also validated GOT2 expression in CM and its correlation with classical prognostic factors in external datasets. Among the GSE3189, GSE15605, and GSE114445 cohorts, the expression level of GOT2 was significantly higher in tumor tissues of CM than in adjacent normal skin tissues (all p < 0.05, Supplementary Figure S4A–C). Based on data from the GSE46517 cohort, we found that the expression of GOT2 was significantly higher in subgroups with more severe clinical predictors, including pathological T stage, distant metastasis, and overall stage (all p < 0.05, Supplementary Figure S4D–F).

Hence, we hypothesized that GOT2 might have a modest effect on tumor biological behavior in cancer cells. We then performed GSEA in CM to evaluate the gene enrichment according to GOT2 expression levels. In terms of KEGG pathways, overexpressed GOT2 was negatively associated with immune-related pathways (Figure 4H), such as “Primary Immunodeficiency” (NES = −2.308, adj. p < 0.001), “Cytokine Cytokine Receptor Interaction” (NES = −2.213, adj. p < 0.001), “Natural Killer Cell Mediated Cytotoxicity” (NES = −2.204, adj. p < 0.001). However, in terms of Retcome pathways, overexpressed GOT2 was positively associated with tumor-related pathways (Figure 4I), such as “MTOR Signaling” (NES = 1.44, adj. p < 0.001), “Cellular Response To Hypoxia” (NES = 1.647, adj. p = 0.004), “Signaling By NOTCH4” (NES = 1.661, adj. p = 0.007).

Since immune-related pathways were primarily enriched via GSEA, we decided to explore the regulatory effect of GOT2 on the TME. We first compared the differences in immune infiltration between high- and low-GOT2 expression groups. The infiltration levels of “NK CD56 bright cells”, “TFH”, “aDC”, “B cells”, “CD8 T cells”, “Cytotoxic cells”, “iDC”, “Macrophages”, “pDC”, “T helper cells”, “T cells”, “Tgd”, “Th1 cells”, and “Treg” were significantly lower in the high-GOT2 expression group (Figure 5A, all p < 0.05). The correlation between GOT2 and immune cells was also identified. The results in Figure 5B revealed that GOT2 was significantly negatively linked with the infiltration of most immune cells (|R|> 0.08, p < 0.05), such as CD8 T cells, B cells, Tgd, and T cells. Then, the ESTIMATE algorithm was conducted to calculate the immune and stromal scores. The high-GOT2 expression group showed lower immune, stromal, and ESTIMATE scores than the low-GOT2 expression group (all p < 0.001, Figure 5C). Our results demonstrated that GOT2 had the potential to become the diagnostic biomarker of CM and may pose as an immune-silent modulator in the TME of CM.

In addition, we also explored the ability of GOT2 as a predictor of immunotherapeutic response in CM. In three GEO cohorts that included immunotherapy-related information, we found that GOT2 expression did not appear to be significantly different between immunotherapy-responsive and non-responsive groups (all p > 0.05, Supplementary Figure S5A–C). Also, in the IMvigor210 cohort, there was no significant difference between the high- and low-GOT2 subgroups in the percentage of cases responding to immunotherapy (p > 0.05, Supplementary Figure S5D). Subsequently, we obtained the IPS scores of immunotherapy response to predict the efficacy of PD-1 and CTLA4 immune checkpoint blockers from the TCIA database. We found that only in the CTLA4-positive-PD-1-positive cohort, the IPS score was significantly higher in the low-GOT2 subgroup than in the high-GOT2 subgroup (Supplementary Figure S5E–H). The results suggested that the efficacy of immunotherapy could be enhanced by the combination of PD-1 inhibitor, CTLA4 inhibitor, and GOT2 inhibitor in CM.

Because of the tumor heterogeneity, we sought to explore further at single-cell resolution. The scRNA-seq data of CM samples were obtained from three GEO datasets. Among the remaining cells after mass filtration, a total of 11701 single cells and 7 cell clusters were derived from tumors and non-malignant samples (Figure 6A). The fundamental markers of each major cell type were shown in the form of a heatmap in Figure 6B. The GOT2 expression was upregulated in several epithelial and stromal cells but downregulated in immune cells (Figure 6C). We further explored the genomic CNA in various cell subtypes to subdivide and identify malignant epithelial cells (Figure 6D). The tSNE plots in Figures 6E, F displayed the distribution of normal and malignant epithelial cells based on the CNA values. According to the expression abundance of GOT2, malignant epithelial cells were divided into GOT2+ and GOT2-malignant cells. Substantially high metabolic activity was observed in GOT2+ malignant cells compared with GOT2-malignant cells and normal cells.

The simulated visualization structure of GOT2 protein was predicted by AlphaFold and was displayed in Figure 7A (Tunyasuvunakool et al., 2021). The HPA database provided the immunohistochemical staining images of GOT2 protein. It was consistent that GOT2 was highly expressed in the CM tissues (Figure 7B). Subsequently, the RT-qPCR was performed and the si-GOT2 molecules were employed to effectively downregulate GOT2 mRNA expression in both A375 and SK-28 cell lines (all p < 0.01, Figure 7C). Furthermore, we used two methods to examine the effect of knocking down GOT2 on cell proliferation potential. Obviously, the cell proliferation was significantly restrained once A375 and SK-28 cells were transfected with si-GOT2 (all p < 0.01, Figures 7D, E). Additionally, the knockdown of GOT2 significantly inhibited the clone formation in A375 and SK-28 cells (Figure 7F). These results implied that GOT2 could promote the growth and proliferation of CM cells.

Based on TCGA data, we conducted the pan-cancer analysis to assess the expression and prognostic significance of GOT2 in various cancers. Differential expression analysis showed that GOT2 expression was significantly increased in 22 cancer types, such as ACC, BRCA, and CESC, compared with normal tissues (Figure 8A). In contrast, GOT2 expression was significantly decreased in 3 cancers, including CHOL, KIRC, and LIHC. Survival analyses showed that beyond CM, GOT2 expression also had prognostic significance in other 8 cancer types (all p < 0.05, Figure). GOT2 may be protective in ACC, LIHC, KIRP, and UCEC, while it is a risk factor for disease progression in LGG, HNSC, MESO, and LAML. All these data indicated that the dysregulation of this glutamic-oxaloacetic transaminase was common across several tumors, but the role of GOT2 in different tumors cannot be generalized, and further exploration was still needed.