Combining single-cell and transcriptomic analysis revealed the immunomodulatory effect of GOT2 on a glutamine-dependent manner in cutaneous melanoma

Background: Reprogramming in glutamine metabolism is a hallmark of cancers, while its role in cutaneous melanoma has not been studied at great length. Methods: Here, we constructed a glutamine metabolism-related prognostic signature in cutaneous melanoma with a variety of bioinformatics methods according to the glutamine metabolism regulatory molecules. Moreover, experimental verification was carried out for the key gene. Results: We have identified two subgroups of cutaneous melanoma patients, each with different prognoses, immune characteristics, and genetic mutations. GOT2 was the most concerned key gene among the model genes. We verified its role in promoting tumor cell proliferation by CCK-8 and clone formation assays. Conclusion: Our study cast new light on the prognosis of cutaneous melanoma, and the internal mechanism regulating glutamine metabolism of GOT2 may provide a new avenue for treating the cutaneous melanoma disease precisely.


Introduction
Melanoma is a potentially fatal malignancy resulting from the transformation and uncontrolled proliferation of pigment-producing melanocytes in the human epidermis (Curti and Faries, 2021).CM is the most aggressive form of melanoma due to its high incidence and metastasis rate.Only 14% of CM patients who develop metastasis survive beyond 5 years (Siegel et al., 2023).For advanced/metastatic CM, chemotherapy is usually ineffective, while radiotherapy is only effective in cases where brain metastases have occurred.A current guideline only recommends targeted therapy for patients with BRAF-mutated CM, whereas other types of mutations do not have sufficient evidence to support treatment recommendations (Seth et al., 2020).Recently, a series of important advances have been made in immunotherapy, such as anti-PD-1 and anti-CTLA-4 antibodies, which greatly reduce the mortality of CM (Seth et al., 2020;Huang and Zappasodi, 2022).Notably, due to specific environmental stresses induced by the treatment, many CM cases can develop genetic progression even under targeted or immunotherapies (Yen et al., 2021).
There is increasing evidence that during cancer development, not only do cancer cells themselves undergo metabolic alterations, but the TME also undergoes metabolic reprogramming due to metabolic competition (Pavlova and Thompson, 2016;Ruocco et al., 2019).Metabolic reprogramming of the TME would lead to invasive migration, metastasis, and even drug resistance of tumor cells, which is an important driver of cancer progression and immune impairment (Viale et al., 2014).For example, glycolysis is one of the key metabolic reprogramming processes in tumor cells, which provides immediate and sufficient energy for rapid proliferation, adaptation to hypoxic environments, and construction of an immunosuppressive TME.Lactate, a major metabolite of glycolysis involved in the construction of an acidic TME, promotes tumor progression and suppresses anti-tumor immunity by inhibiting T cells and PD-L1 upregulation (Feng et al., 2017).Therefore, metabolic reprogramming of tumor cells and immune cells is one of the key barriers to tumor immunotherapy (Kouidhi et al., 2018).
In hyperproliferative cells, glutamine plays an essential role in providing intermediates to the TCA cycle, which contributes to cell proliferation and biosynthesis (Yoo et al., 2020).As well as serving as a nitrogen donor for purines and pyrimidines, glutamine also serves as a precursor for GSH synthesis (Zhu et al., 2022).Glutamine is transported across membranes with the assistance of amino acid transporters and hydrolyzed into glutamate and ammonia through the catalysis of glutaminase.Then glutamate becomes α-KG catalyzed by dehydrogenase or aminotransferase and enters the TCA cycle (Yang et al., 2017).Meanwhile, glutamate can be oxidized into GSH via the glutamate-cysteine ligase, which participates in the immune defense, nutritional metabolism, and antioxidant function of cells by neutralizing mitochondrial ROS (Wang et al., 2019).Previous studies have reported that CM upregulates mitochondrial oxidative phosphorylation by enhancing glutamine anaplerotic metabolism, thereby promoting tumor cell resistance to targeted therapy (Haq et al., 2013;Hernandez-Davies et al., 2015;Baenke et al., 2016;Zhang et al., 2016;Vashisht Gopal et al., 2019).Several small molecule inhibitors have been utilized to take advantage of the "glutamine fragility" in cancer treatment due to the crucial role it plays in metabolic reprogramming (Li et al., 2019).Nonetheless, the clinical development of these drugs as monotherapies has been limited by toxicity or ineffectiveness (Wu et al., 2018;Leone et al., 2019).
It is important to note that glutamine is not just an important nutritional supplement for tumor cells, but also an essential component of immune cell function, including its role in activating T lymphocytes (Wang and Green, 2012;Cruzat et al., 2018).A recent study found that glutamine metabolizing enzyme inhibitors could promote Th1 and cytotoxic T-cell proliferation and viability by altering the epigenome (Johnson et al., 2018).The consumption and metabolism of glutamine may differ between tumor cells and TILs.Targeting the fragility of glutamine metabolism in tumor cells may change the fate of the TME and strengthen the anti-tumor ability of the TILs (Halama and Suhre, 2022;Ma et al., 2022).
In the current study, we established a glutamine metabolism risk score model based on GMRGs and evaluated the immune characteristics and tumor burden characteristics of the high-and low-GMRS groups under this characterization.The mechanism by which GOT2 regulates glutamine metabolism in tumor cells may provide new insights into targeted therapy for CM.Abbreviations and corresponding words and phrases used in this article can refer to Supplementary Table S1.

Candidate data sources
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.

Construction of the prognostic GMRG signature
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.

Subgroup analysis of the gene signature
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).

Clustering analysis
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).

Biological and mutational analysis
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 immunerelated 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.

Single-cell RNA-seq data processing
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.

Cell culture and proliferation assay
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 2 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.

Statistical analysis
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.

Establishment of the GMRGs related prognostic signature
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: 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).

The relationship between GMRS and clinicopathologic factors
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).

Mutation burden between high-GMRS and low-GMRS groups
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).

Immune microenvironment difference between High-GMRS and Low-GMRS groups
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.

Identification of glutamine-metabolismrelated molecular subtypes
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 "collagencontaining 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.

Identification of the potential key gene in 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).

Potential function of GOT2 in regulating TME
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.

Mapping GOT2 in single-cell data
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 GOT2malignant cells.Substantially high metabolic activity was observed in GOT2+ malignant cells compared with GOT2-malignant cells and normal cells.

Verification of GOT2 expression in CM
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.

Pan-cancer analysis of GOT2 expression and prognostic significance
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.CM shares the oncogenic transformation features of glutamine addiction, like many other solid malignancies (Filipp et al., 2012).In other words, CM cells depend on glutamine for growth, irrespective of their oncogenic background.In a study quantitatively assessing glutamine metabolism in melanoma cell lines, it was shown that energy-producing anaplerosis and asparagine biosynthesis are responsible for CM cell growth (Ratnikov et al., 2015).The glutamine metabolism in CM cells is characterized by overexpression of genes involved in the production of proline from glutamate, thus increasing proline production by tumor cells in comparison with melanocytes (De Ingeniis et al., 2012).It is these metabolic changes that underlie the development of CM and that affect its response to treatment.Most tumors, in the meantime, trigger their own growth and development by remodeling the TME and recruiting corresponding cells, like TAMs and Tregs, which usually lead to poor outcomes (Watson et al., 2021;Shi et al., 2022).In CM, the TME consists of adjacent cells such as keratinocytes, adipocytes, immune cells, CAFs, and extracellular matrix elements (Mazurkiewicz et al., 2021).There is a strong correlation between immune cell infiltration in the TME and glutamine metabolism, according to several recent studies (Oh et al., 2020;Yang et al., 2021).
In our study, genes involved in glutamine metabolism in CM were based to generate a prognostic signature and evaluate its predictive value.The model possessed independent predictive power and was strongly correlated with other commonly used clinicopathological factors.In terms of immune function and infiltration, we found that patients in the high-GMRS group were developing an immune-silent microenvironment, which partly explained their poorer prognosis.By clustering CM cases according to the genes used for modeling, we defined two glutamine metabolism-related clusters.Among them, the CM patients in cluster A owned a high degree of overlap with the patients in the high-GMRS group, and both possessed poor prognoses.Mostly tumor-related signaling pathways were enriched in DEGs between the two clusters.Our findings further provided new insights into glutamine metabolism in melanoma.
Among these model genes, GOT2 stood out as an important enzyme for cellular redox homeostasis and aspartate production.As a crucial component of protein, purine, and pyrimidine nucleotide synthesis, aspartate also plays a crucial role in cell growth (Garcia-Bermudez et al., 2018).In our study, CM patients who expressed high levels of GOT2 had poor prognoses and low immune infiltration.
According to previous studies, GOT2 promoted the growth and proliferation of malignant tumors mainly through three pathways (Yang et al., 2015;Yang et al., 2018;Hong et al., 2019;Guan et al., 2021;Abrego et al., 2022).The first one is the supply of cellular building materials.Hong et al. reported that the protein BRCA1 and ZBRK1 can together form a complex that inhibits GOT2 transcription and translation.Due to the impairment of this complex, the production of GOT2 increased, which resulted in the rapid proliferation of breast cancer cells (Hong et al., 2019).The second one is the protection from oxidative damage.It is indispensable for cells to create NADPH and maintain the cellular redox state by glutamine carbon flow through GOT2.The knockdown of GOT2 in pancreatic ductal adenocarcinoma resulted in an increased production of ROS, which led to the cyclin-dependent kinase inhibitor p27-dependent senescence (Yang et al., 2018).The third one is the suppression of immune function.PPARδ is a lipid metabolism-related transcription factor.According to Abrego et al., GOT2 can shuttle to the nucleus and enhance PPARδ activity (Abrego et al., 2022).The translation products of many genes activated by PPARδ possess immunosuppressive properties, resulting in low immune infiltration as well as high tumor burden.We speculated that GOT2 could also play a cancer-promoting role through the above three pathways in CM.
However, GOT2 could exhibit a cancer-suppressive profile in some cancers.For instance, HCC cells exhibit downregulation of GOT2, and low GOT2 expression is associated with advanced disease progression (Li et al., 2022).Mechanistically, the reduction of GOT2 in HCC mediated the reprogramming of glutamine metabolism towards the synthesis of reduced GSH, which maintained redox homeostasis in tumor cells by resisting the ROS damage in HCC progression and stimulating the PI3K/AKT/mTOR signaling pathway thereby contributing to the malignant progression of HCC.The mechanisms behind GOT2 are still being explored, but it appears it is involved in reprogramming glutamine metabolism in order to promote cancer progression.This can be used as a therapeutic and diagnostic target for CM.
Admittedly, our research has some limitations.Our risk model included 5 genes, which may increase the cost of testing.We may be able to involve multiple clinicopathological parameters in the construction of the model to increase the practicality.The effect of GOT2 on immune function and immune cell infiltration may need more experiments to verify.

Conclusion
To summarize, we built up a brand-new prognostic model and stratified CM patients according to the model scores (GMRSs).Significant differences were found in prognosis, immune characteristics, and genomic mutation between these subgroups.Then, we identified the glutamine-metabolism-related molecular subtypes with different biological features.Furthermore, GMRGs like GOT2 could contribute to an in-depth understanding of the underlying mechanisms of CM and may become a new independent biomarker and target for the diagnosis and treatment of CM.

FIGURE 1
FIGURE 1 Development of the GMRS model.(A) Nine prognosis-related genes were selected by the univariate Cox regression analysis.(B,C) Nine prognosisrelated genes underwent the LASSO regression analysis.(D) Five genes were selected with the best value of coefficient by the multivariate Cox regression analysis.(E) The K-M survival curve shows the poorer survival probability of CM patients in the high GMRS group (p < 0.001).(F) The AUC value of the ROC curve for predicting CM patients' prognosis.(G,H) Predictive independence assessment of classical clinical predictors and our GMRS score by the (G) univariate and (H) multivariate Cox regression analysis.

FIGURE 2
FIGURE 2 Interpretation of the gene mutation landscape between GMRS groups.(A,B) Oncoplots showing the top 20 mutation profiles in the (A) high-GMRS and (B) low-GMRS groups.The types of mutation in each gene and the corresponding percentage in the samples were also included.(C-E) The differential expression of the model genes among the subgroups stratified by the three top gene mutations (C) TTN; (D) MUC16; (E) BRAF).(F) The K-M survival curve shows the poorer survival probability of CM patients in the low TMB group.(G) The K-M survival curves show the survival differences between subgroups stratified by both TMB and GMRS scores in the CM cohort.* represents statistical p-value < 0.05, ** represents statistical p-value < 0.01, *** represents statistical p-value < 0.001.

FIGURE 3
FIGURE 3 Comparison of the immune characteristics between risk groups.(A) Heatmap showing the distribution of tumor-infiltrating immune cells and GMRS scores by 7 mainstream algorithms.(B) Spearman correlation analysis of the tumor-infiltrating immune cells under GMRS score in CM. (C) Differential abundance analysis of tumor-infiltrating immune cells between GMRS groups using ssGSEA algorithm.(D) Differential immune function analysis between GMRS groups.(E) Differential analysis of immune-related scores between GMRS groups using the ESTIMATE algorithm.(F) Differential analysis of the anticancer immune functions of the cancer-immunity-cycle between GMRS groups.* represents statistical p-value < 0.05, ** represents statistical pvalue < 0.01, *** represents statistical p-value < 0.001.

FIGURE 4
FIGURE 4 Identification of the potential key gene in CM by Xiantao tool.(A) Differential expression analysis of the model genes between CM and normal tissues.(B) ROC curves of diagnostic efficiency of the model genes.(C-G) K-M survival curves of the CM patients stratified by expression of model genes (C) GLUD2; (D) DGLUCY; (E) SLC6A14; (F) LAT; (G) GOT2), respectively.(H-I) Enrichment analysis of (H) KEGG and (I) Retcome pathways in CM patients according to the GOT2 expression.* represents statistical p-value < 0.05, ** represents statistical p-value < 0.01, *** represents statistical p-value < 0.001.NES, normalized enrichment score.P.adj, adjusted p value.

FIGURE 5
FIGURE 5 Comparison of the immune characteristics in GOT2 subgroups by Xiantao tool.(A) Differential abundance analysis of tumor-infiltrating immune cells between GOT2 groups using the ssGSEA algorithm.(B) Spearman correlation analysis of the tumor-infiltrating immune cells under GOT2 expression in CM. (C) Differential analysis of immune-related scores between GOT2 subgroups using the ESTIMATE algorithm.* represents statistical p-value < 0.05, ** represents statistical p-value < 0.01, *** represents statistical p-value < 0.001.

FIGURE 6
FIGURE 6 Mapping GOT2 in single-cell data.(A) A tSNE cluster view of 11701 cells was obtained from three GEO datasets.The identified clusters are annotated according to their origin.(B) The expression of marker genes for each cell type.(C) The density and distribution of GOT2 expression in different clusters.(D) Copy number analysis of the epithelial cells with the R package inferCNV.(E,F) The distribution of malignant cells and normal epithelial cells is based on the CNA values.(G) The multilineage interactome network among different cell clusters.

FIGURE 7
FIGURE 7 Verification of GOT2 expression in CM. (A) An accurately predicted complex structure of GOT2 protein by AlphaFold.(B) IHC staining of GOT2 in clinical CM tissues and normal skin tissues.(C) Verification of the knockdown efficiency of si-GOT2s by evaluating GOT2 mRNA expression using RT-qPCR.(D-E) Cell proliferation curve following transfection of si-GOT2 in (D) A375 and (E) SK-28 cells.(F) Comparison of the clone formation ability between negative control groups and GOT2 knockdown groups in A375 and SK-28 cells.NC, negative control; KD, knockdown; OD, optical density.* represents statistical p-value < 0.05, ** represents statistical p-value < 0.01, *** represents statistical p-value < 0.001.

FIGURE 8
FIGURE 8 Pan-cancer analysis of GOT2 expression and prognostic significance.(A) Differential expression analysis of GOT2 between tumor tissues and corresponding adjacent normal tissues in 33 cancer types.(B) K-M survival curves of patients of 8 cancer types stratified by expression of GOT2.* represents statistical p-value < 0.05, ** represents statistical p-value < 0.01, *** represents statistical p-value < 0.001.