The Overexpression of SLC25A13 Predicts Poor Prognosis and Is Correlated with Immune Cell Infiltration in Patients with Skin Cutaneous Melanoma

Purpose Skin cutaneous melanoma (SKCM) is one of the most malignant and aggressive cancers with poor prognosis due to its rapid progression towards metastasis. Thus, finding clinically relevant biomarkers for early diagnosis, prognosis, and therapy prediction is essential. This study focused on the identification of SLC25A13 as a novel biomarker for SKCM and is aimed at investigating the biological functions of solute carrier family 25 member 13 (SLC25A13) in the development of SKCM. Methods GEPIA was used to analyze the diagnostic and prognostic values of SLC25A13 in SKCM using the TCGA dataset. PrognoScan was used to validate the prognostic value of SLC25A13 and its coexpressed genes in SKCM. TISIDB was established to reveal the relationship between the expression of SLC25A13 and immune infiltration in SKCM. The protein expression of SLC25A13 in SKCM was evaluated by the Human Protein Atlas. The signaling pathways and biological functions of SLC25A13 in SKCM were analyzed by LinkOmics. Metascape was applied to analyze the functional enrichment analysis of SLC25A13. Protein-protein interaction analysis of SLC25A13 was performed by GeneMANIA. Results The mRNA and protein levels of SLC25A13 in the SKCM were much higher than those in the normal tissue. Furthermore, the overexpression of SLC25A13 predicts worse outcomes of SKCM patients. Moreover, the SLC25A13 expression was negatively correlated with the immune infiltration level of SKCM. The overexpression of SLC25A13 coexpressed genes, such as ACLY and AFG3L2, and SCL25A13 interacting genes also predicted the unfavorable prognosis of SKCM patients. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of SLC25A13 coexpressed genes showed that these genes are enriched in ATPase activity, cell cycle, mTOR, and VEGFA-VEGFR2 signaling pathways, which were relevant to tumor development and angiogenesis. Gene set enrichment analysis (GSEA) demonstrated that the SLC25A13 expression was related to infiltrating immune cells in SKCM. Conclusion Our findings revealed that SLC25A13 might be a potential prognostic and therapeutic biomarker for SKCM.


Introduction
The incidence of skin cutaneous melanoma (SKCM) is increasing in both males and females, and SKCM caused approximate 72% of deaths in skin carcinoma [1]. SKCM lacks effective treatment except early surgical resection, and it usually results in a poor prognosis and shockingly high mortality [2]. Meanwhile, SKCM is highly heterogeneous, which makes personalized treatment difficult [3]. Immune checkpoint inhibitors, such as cytotoxic T-lymphocyte antigen-4 (CTLA-4) and programmed cell death protein 1 (PD-1) inhibitors, are the main treatment measures for advanced-stage SKCM, but the efficacy of immunotherapy is limited due to different response rates [4]. In addition, primary and secondary resistance and the absence of predictive markers of response are biological challenges and clinical dilemmas in immunotherapy for SKCM patients [5]. Thus, it is urgent to identify a novel predictive biomarker for the immunotherapy of SKCM patients.
Solute carrier family 25 member 13 (SLC25A13), a calcium-binding/stimulated aspartate-glutamate carrier, plays a key role in metabolic pathways including aerobic glycolysis, gluconeogenesis, urea cycle, and synthesis of proteins and nucleotides [6]. The genetic deficiency of SLC25A13 causes adult-onset type II citrullinemia [7]. Furthermore, pathogenic germline variants correlated with increased tumor mutational burden have been found in 7 genes, including SLC25A13 and TP53, and germline variants can predict the response to immune checkpoint inhibitors [8].
In addition, SLC25A13 is correlated with tumor aggressiveness and poorer prognosis of colorectal cancer [9]. Thus, it is of importance to reveal the role of SLC25A13 in cancer development and evaluate the prognostic value of immunotherapy of SKCM patients.
This study compared the gene transcription levels between SKCM and Paracancerous tissue using RNA sequencing (data not shown), and the top 100 genes overexpressed in SKCM were considered as potential oncogenes involved in cancer process and verified by the cancer genome atlas (TCGA) data sets of SKCM. In this study, we investigated the diagnostic and prognostic value of SLC25A13 for SKCM patients and clarified the correlation between the expression of SLC25A13 and immune cell infiltration in SKCM patients; meanwhile, we also predicted the potential gene function of SLC25A13 by enrichment analysis of its neighbor genes in SKCM ( Figure 1).

Analysis of the Expression of SLC25A13 in SKCM and
Normal Skin Tissues. Gene Expression Profiling Interactive Analysis (GEPIA, http://gepia.cancer-pku.cn/) is a webbased tool, and The Cancer Genome Atlas (TCGA) and Genotype Tissue Expression (GTEx) database provide fast and customizable functionalities [10]. In this study, the expression of SLC25A13 in SKCM and normal skin tissues was compared and represented by box plots via the GEPIA. In addition, the overall survival and disease-free survival of SKCM patients based on the expression of SLC25A13 were analyzed using GEPIA.
Human Protein Atlas (https://www.proteinatlas.org) provides the expression profiles in human tissues of genes both on the mRNA and protein level [11]. The protein levels of SLC25A13 were compared in SKCM and normal skin tissues by Human Protein Atlas.

Survival Analyses in SKCM Patients
. PROGgene (http:// www.compbio.iupui.edu/proggene) serves as a tool to perform survival analysis on single genes, multiple genes as a signature, ratio of expression of two genes, and curated/published gene signatures [12]. The overall survival of SKCM patients based on the expression of SLC25A13 was analyzed using PROGgene. PrognoScan (http://dna00.bio.kyutech.ac.jp/PrognoScan/) shows the prognostic values of genes in multiple microarray datasets with meta-analysis [13]. PrognoScan was used to explore the prognostic value of interacting genes of SLC25A13 in SKCM patients.

Clinical Prediction Value of the Established Prognostic
Model. The TCGA database (https://portal.gdc.cancer.gov/) and the UCSC Xena database (https://xena.ucsc.edu/) were used to obtain SKCM data. All subsequent analyses were performed via R version 4.1.3. To test whether the expression of SLC25A13 was the effective prognostic indicator, the univariate and multivariate Cox regressions were performed [14,15]. Then, diagnostic ROC and PR analysis of the SLC25A13 expression in normal samples and SKCM samples was performed by R package "ROCR" [16]. Finally, to assess the risk of death for SKCM, a nomogram was constructed by combining the expression of SLC25A13 and clinical data. The patient "TCGA-EE-A2A2-06A" was illustrated.

Correlation
Analysis between the Expression of SLC25A13 and Immune Features in SKCM Patients. TISIDB (http://cis.hku.hk/TISIDB/index.php) is a website portal for interaction between tumor and immune systems in 28 types of tumor-infiltrating lymphocytes across human cancers [19]. The correlation between the abundance of tumorinfiltrating lymphocytes and the expression of SLC25A13 was investigated by TISIDB, and the relationship between abundance of immunomodulators and the expression of SLC25A13 was determined using TISIDB.

Analysis of Coexpressed Genes of SLC25A13 in SKCM.
The multiomics data from all 32 TCGA cancer types and 10 Clinical Proteomics Tumor Analysis Consortium (CPTAC) cancer cohorts are available at LinkedOmics (http://www.linkedomics.org) [20]. In this study, LinkFinder module was used to investigate the coexpression genes of SLC25A13, and the LinkInterpreter module was employed to perform enrichment analysis of co-expression genes associated with SLC25A13 in SKCM.
The STRING database (http://string-db.org/) aims to predict protein-protein associations including direct (physical) and indirect (functional) relations from knowledge transfer between organisms and interactions aggregated from other databases [21]. Gene Ontology (GO) and KEGG analysis of SLC25A13 coexpressed genes were evaluated by the STRING database to predict the biological functional of SLC25A13.
Cytoscape is an open source software project for integrating biomolecular interaction networks with highthroughput expression data and other molecular states into a unified conceptual framework [22]. Cytoscape and MCODE were established to estimate the most indispensable module in the protein-protein interaction networks.  3 Disease Markers annotation, and membership search to leverage more than 40 independent knowledge bases within one integrated portal [24]. In this study, Metascape was used to demonstrate the functional enrichment analysis of SLC25A13 and its neighboring genes.
GSEA is a powerful analytical method to perform genome-wide expression profiling of two phenotypes and evaluate cumulative changes in the expression of groups of multiple genes defined based on previous biological knowledge [25]. In this study, GSEA was applied to evaluate the biological functions and signaling pathways regulated by SLC25A13.

SLC25A13 Is Overexpressed in SKCM Compared with
Normal Skin Tissues. Firstly, the levels of SLC25A13 in SKCM and normal skin tissues were analyzed by immunohistochemistry data collected from Human Protein Atlas [11]. As shown in Figure 2(a), high levels of SLC25A13 were detected in SKCM samples, whereas the expression of SLC25A13 could not or barely be detected in normal skin tissues. Statistical analysis manifested that the expression of SLC25A13 in the SKCM was much higher than that in normal skin tissues (Figure 2(b)). In line with the protein expression, the mRNA levels of SLC25A13 were also upregulated in SKCM compared with normal skin tissues ( Figure 2(c)) [10].

High Levels of SLC25A13 Predict Poor Disease Outcome in SKCM.
To find whether the SLC25A13 expression is associated with overall survival (OS) and disease-free survival (DFS), we compared the OS and DFS of SKCM patients with high or low SLC25A13 expression. As shown in Figures 3(a) and 3(b), patients with low levels of SLC25A13 had longer OS (p = 0:00011) and DFS (p = 0:036) [10]. Furthermore, this was also validated using PROGgene based on GSE19234 and GSE22153 datasets. As shown in

Disease Markers
Figures 3(c) and 3(d), the overexpression of SLC25A13 was associated with worse OS in SKCM patients [13,26,27]. To explore the prognostic value of SLC25A13 gene in SKCM, univariate COX analysis was performed, and SLC25A13 was found to be an independent prognostic factor ( Figure 3(e), p < 0:05). SLC25A13 was also found to be an independent prognostic factor by multivariate COX analysis ( Figure 3(f), p < 0:05). To explore the diagnostic value of SLC25A13 gene in SKCM, diagnostic ROC analysis of SLC25A13 gene in 556 normal subjects and 468 SKCM patients was performed. ROC analysis showed that the area under the curve (AUC) was 0.849 (Figure 3(g)). Meanwhile, the area under the curve (AUC) was 0.894 according to PRC analysis. These results indicated that SLC25A13 might be a good diagnostic molecule for SKCM. To further determine the survival of SKCM patients, we drew a nomogram combining the risk value and clinical characteristics of the model [28]. As shown in Figure 4, we found that the 1-, 3-, and 5-year mortality rates of patient "TCGA-EE-A2A2-06" were 0.102, 0.478 and 0.64, respectively.

Analysis of Coexpressed Genes Correlated with SLC25A13 in SKCM.
To investigate the SLC25A13 coexpressed genes, we identified the potential candidate genes using LinkedOmics [20] and measured their Pearson coefficient of correlation with SLC25A13. The top 50 significant genes that were negatively (Figure 6(a)) and positively ( Figure 6(b)) correlated with SLC25A13 were presented in the heat map. As shown in Figure 6(c), 8959 genes (with red label) showed statistically significant positive correlations with SLC25A13, while 10769 genes (with green label) showed statistically significant negative correlation with SLC25A13.
3.5. The Prognostic Value of SLC25A13 Associated Genes for SKCM Patients. A total of 295 coexpressed genes with strong positive correlations with SLC25A13 were used to construct a protein-protein interaction network, and the vital module in the network was highlighted in yellow ( Figure 6(d)). Furthermore, the prognostic values of SLC25A13 coexpressed genes in SKCM were investigated using PrognoScan [13]. To better understand the biological functions of SLC25A13, the protein-protein interaction network of SLC25A13 and its interacting genes was constructed using GeneMANIA (Figure 8(a)) [23]. Meanwhile, the prognostic values of the interacting genes of SLC25A13 were analyzed using PrognoScan [13]. As shown in Figures 8

GO and KEGG Analysis of SLC25A13 Coexpressed
Genes. We input 649 coexpressed genes of SLC25A13 into the STRING database to perform functional enrichment analysis [29]. GO analysis of biological processes revealed that these genes were enriched in mitotic spindle midzone assembly, spindle midzone assembly, DNA replication checkpoint, and pigment biosynthetic process (Figure 9(a)). These genes were located in melanosome membrane, melanosome, and mitochondrial nucleoids (Figure 9(c)) and associated with the molecular functions of retromer complex binding, DNA replication origin binding, ATPase activity, etc. (Figure 9(b)).

Disease Markers
The results from KEGG analysis indicated that pathways including cell cycle and mTOR signaling pathway were particularly enriched, which were closely related to cancer progression (Figure 9(d)).

Functional Enrichment Analysis of SLC25A13 in Patients
with SKCM. Metascape analysis showed that SLC25A13 and its neighboring genes were enriched in glucose metabolism, amino acid metabolism, and cell cycle (Supplementary Figure 1A-1B) [24]. Moreover, they were also enriched in pathways such as VEGFA-VEGFR2 signaling pathway and urea cycle and associated pathways, which were relevant to tumor development and angiogenesis. Meanwhile, SLC25A13 and its neighboring genes were enriched in transcription factor target genes, including MEF2C, NCOA2, SNAI1, and SOX10 target genes, which were participated in the cancer progression ( Supplementary Figure 1 C-D). GSEA of SLC25A13 was utilized to further investigate the biological function of SLC25A13 in SKCM by Linke-dOmics [20]. GO analysis of molecular function showed that SLC25A13 positively participated in ATPase activity, while negatively regulated antigen binding and immunoglobulin binding (Supplementary Figure 2). GO analysis of the biological process showed that SLC25A13 positively participated in pigment metabolic and negatively regulated T cell activation, lymphocyte mediated immunity, humoral immune response, response to tumor necrosis factor, and leukocyte proliferation (Supplementary Figure 3). GO analysis of cell components displayed that SLC25A13 positively regulated pigment granules but negatively regulated immunological synapse and MHC protein complexes in SKCM (Supplementary Figure 4). KEGG pathway analysis showed that SLC25A13 positively regulated Fanconi anemia pathway, cell cycle, DNA replication, basal transcription factors and negatively regulated natural killer cellmediated cytotoxicity and antigen processing and presentation (Supplementary Figure 5).

Discussion
SKCM is a highly aggressive skin tumor worldwide due its metastatic nature, and the incidence rate of SKCM keeps increasing in recent years [30]. Once SKCM spreads through the dermis, it easily spreads to vital organs and lymph nodes and causes poor prognosis [31]. Therefore, it is urgent to identify biomarkers to predict the prognosis of SKCM. This study indicated that the mRNA level and protein expression of SLC25A13 was significantly higher in SKCM than that in normal skin tissues. Moreover, Cox regression analysis showed that the SLC25A13 overexpression was an independent risk factor for poor outcome of SKCM patients. Furthermore, the overexpression of SCL25A13 associated genes including ACLY, AFG3L2, ASL, CHCHD3, CTH, GFM1, MDH2, NDUFS1, SLC25A15, GOT1, GOT2, RNF31, SIRT7, TUBG1, and ELAVL1 were also related to poor prognosis of SKCM patients. These findings revealed that SLC25A13 is a promising and prognostic biomarker in SKCM.
Multiple members of the mitochondrial carrier family (SLC25) are involved in the invasion and metastasis of many cancers. For example, Wong et al. reported that SLC25A22 promotes the proliferation and migration of colorectal cancer cells, and high levels of SLC25A22 predict poor prognosis in colorectal cancer [32]. SLC25A13 was also found to be pivotal in colorectal cancer aggressiveness [9]. Although the function of SCL25A13 in SKCM development may need further investigation, the function enrichment analysis revealed that SLC25A13 and its neighboring genes were enriched in cell cycle, VEGFA-VEGFR2 signaling, and mTOR pathway, which were related to the development of SKCM. SKCM develops from mutation of pigment-producing cells in the skin, and there is uncontrolled growth of pigmentproducing cells during SKCM development [33]. In this study, the function enrichment analysis of SLC25A13 and its neighboring genes revealed that SLC25A13 was involved in pigment granules, pigment biosynthetic, and metabolic processes. Our results suggested that SLC25A13 might be critically involved in the progression of SKCM. It is well known that cancer cells have a higher metabolic rate than normal cells as a result of environmental selection and the accumulation of mutations [34]. Then, we downloaded the mutation data of SKCM through mafTools package and showed the landscape map of the top 10 genes with the highest mutation frequency, among which TTN and MUC16 genes are more prone to mutation (Supplementary Figure 6A). However, there is no relationship between SLC25A13 and tumor mutation load (TMB) through correlation analysis (Supplementary Figure 6B, COR = −0:07, p > 0:05). Thus, SLC25A13 was not associated with TMB of SKCM.
SLC25A13 plays a pivotal role under glucose-deprived conditions and is correlated with tumor aggressiveness [9]. The presence and amount of glucose have a great impact on the immune cell functions of the immune system, resulting in the tumor pathological stage [35]. Thus, we were interested in investigating the relationship between the expression of SLC25A13 and immune cell infiltration of SKCM. The data of the current study revealed that the SLC25A13 expression was negatively correlated with immune cell infiltration, including effector memory CD8+ and CD4+ T cells, Th1 cells, Th2 cells, and Treg cells. Furthermore, the function enrichment analysis revealed that SLC25A13 negatively regulates T cell activation, lymphocyte mediated immunity, humoral immune response, natural killer cell-mediated cytotoxicity, etc. These data indicated that the SLC25A13 expression was negatively correlated with immune cell infiltration in SKCM. Furthermore, SKCM treatment has been revolutionized with the approval of immune checkpoint inhibitors, which have a significant impact on the prognosis of patients with SKCM [36]. Nevertheless, only a small percentage of SKCM patients can benefit from immunotherapy [37]. Tumor immune cell infiltration can affect the sensitivity of immunotherapy, and higher immune cell infiltration scores showed a significant immune therapeutic advantage and clinical benefit [38]. Thus, the expression of SLC25A13 might be a potential biomarker useful for screening suitable SKCM patients for immunotherapy. Comprehensive psychological intervention can remarkably improve the negative emotions, social support degree, sleep quality, and immune functions of cancer patients receiving conventional nursing intervention during the perioperative period [39]. As SLC25A13 could be a 13 Disease Markers biomarker for tumor immune cell infiltration in SKCM, the feasibility of exploring SLC25A13 as a biomarker of psychological intervention in SKCM treatment might need further investigation.

Conclusions
This work demonstrated that the SLC25A13 expression was increased in SKCM compared with normal skin tissues, and the high expression of SLC25A13 and its neighbor genes were associated with a poor prognosis in SKCM patients. In addition, the expression of SLC25A13 was negatively correlated with the tumor-infiltrating cells in the tumor immune microenvironment. Functional enrichment analysis of SLC25A13 showed that it was associated with tumor development and progression. These findings suggest that SLC25A13 may be a potential prognostic biomarker for the diagnosis and prognosis of SKCM patients.

Data Availability
All data generated or analyzed used to support the findings of this study are included within the article and supplementary materials.

Conflicts of Interest
The authors report no conflicts of interest in this work.