SLC13A4 Might Serve as a Prognostic Biomarker and be Correlated with Immune Infiltration into Head and Neck Squamous Cell Carcinoma

SLC13A4 is a sodium sulfate co-transporter, which is expressed in brains, placentas, thymes and other tissues, plays an essential role in maintaining the metabolic balance of sulfate in vivo. The TCGA database shows that it is differentially expressed in a variety of tumors, but its prognostic value in tumors has not been clarified. TCGA, Oncomine and Timer databases were used to analyze SLC13A4 mRNA expression in cancer tissues and normal tissues, and its correlation with clinical prognosis in head and neck tumor. The CIBERSORT database was used to analyze the correlation between SLC13A4 expression and the infiltration of immune cells. SLC13A4 enrichment analysis was carried out by GSEA. SLC13A4 mRNA levels were significantly lower in head and neck tumors than in paracancer tissues. SLC13A4 expression in Head and neck squamous cell carcinoma (HNSCC) was closely related to tumor pathological grade and clinical stage. Decreased SLC13A4 expression was associated with poor overall survival (OS), progression free survival (PFS), disease specific survival (DSS) and recurrence free survival (RFS) in HNSCC patients. The expression of SLC13A4 was negatively correlated with Monocytes, M1 macrophages, M2 macrophages, resting CD4+ memory T cells, resting NK cells and activated NK cells, but positively correlated with neutrophils, plasma cells, T follicular helper cells, gamma delta T cells, regulatory T cells and naive B cells. In addition, the genes in SLC13A4 low-expression group were mainly concentrated in immunity-related activities, viral diseases, typical tumor pathways and metabolism. The SLC13A4 high expression group was mainly enriched in metabolic pathways. These suggest that SLC13A4 may be a potential prognostic biomarker in HNSC and correlated with immune infiltrates.


INTRODUCTION
Most head and neck cancers are squamous cell carcinomas that originate in the upper respiratory epithelia (located in oral cavity, oropharynx, larynx or hypopharynx). Head and neck squamous cell carcinoma (HNSCC) tend to develop in the mucous membrane of the mouth, nose, and throat. It is the sixth most common cancer around the globe, with a 5-year mortality rate of 50% [1,2]. Although the application of various treatment modalities, and combination use of surgical, radiotherapeutic and chemotherapeutic treatments can significantly reduce the short-term and long-term survival of the malignancy, the situation remains grim [3]. The immune cells and stromal infiltration contribute to the tumor microenvironment (TME), in which various immune cells act as either promoters or inhibitors of cancers, and they are closely related to the prognosis of patients [4]. Intensive researches demonstrated that PD-1 inhibitors exert a significant therapeutic effect on HNSCC [5,6]. In the clinical practice, early diagnosis and accurate prognostic prediction of HNSCC patient remain a challenge due to the lack of specific biomarkers. To improve the long-term survival of tumor patients, it is of great importance to identify the biomarkers that are conductive to early diagnosis, precise outcome prediction of cancers and to the identification of people who well respond to immunotherapy.
A feature shared by all members of solute carrier family 13 (SLC13) is the multi-bonded proteins encoding 8-13 transmembrane domains, and is found in a wide array of tissues. SLC13 family has 5 members. In terms of gene codes and function, they fall into two categories: sodium sulfate cotransporter (SLC13A1, SLC13A4) and sodium carboxylic acid co-transporter (SLC13A2, SLC13A3, SLC13A5) [7]. SLC13A4 is one of the major transporters of sulfate metabolism and plays an anti-apoptotic part by mediating the transmembrane transportion of thiosulfate on cell membrane, regulating the concentration of intracellular sulfate and modulate caspase-3 sulfurization [8]. SLC13A4 is very different from other sulfate transporters, and its functions have been preliminarily explored in tonsil, placenta, testis and central nervous system, but the physiological implication of its low-level expression in heart, thymus and liver has yet to be elucidated [9][10][11][12]. Moreover, its roles in tumors remains poorly understood. The mutated sulfate transporter gene and/or disordered sulfate metabolism are closely related to multiple human diseases, such as aberrant early embryonic development, chondrodysplasias and Hunter syndrome, which may be the predisposing factors of neurodevelopmental diseases, such as mental retardation and autism spectrum disorder [11,[13][14][15]. Thiosulfate is one of the important inorganic sulfates in the body and the fourth most abundant anion in mammalian plasma, and is crucial for a great many physiological functions [16,17]. In prostate cancer, the content of thiosulfate is used as an auxiliary diagnostic marker [18]. All these highlight the important role that sulfate transporter protein plays in human physiology. In addition, thiosulfate is an important product of H2S metabolism. H2S was reported to induce DNA damage and affect the cell signal transduction and cell cycle, and participate in a variety of physiological and pathological processes, including tumorigenesis and development [19][20][21]. These further indicate that the sulfate process of metabolic pathway may be closely involved in TME.
This study examined the differential expression, clinical relevance and prognostic value of SLC13A4 in HNSCC by searching RNA-seq and clinical data in the TCGA public database. On the basis of its expression in thymus, we further analyzed the potential pathways and possible molecular mechanisms by which SLC13A4 is involved in HNSCC tumor immune infiltration by searching cibersort and GSEA database.

TCGA Database
All level 3 RNA expression data (Workflow Type: HTSeq-FPKM) from 502 HNSCC cancer tissues and 44 adjacentcancer tissues were retrieved from the Cancer Genome Atlas (TCGA) database (https://tcga-data.nci.nih.gov/tcga/) [22]. The clinical information of HNSCC patients, including age at diagnosis, gender, HPV status, radiation therapy, pathological stage, T stage, N stage, M stage, grades, overall survival time, overall survival status, disease free survival time, disease free survival status, disease specific survival time, disease specific survival status, progression-free survival time, progression-free survival status were also downloaded from the TCGA. TCGA-HNSCC has 500 HNSCC patients with clinical information, but the information for each patient is not complete. In detail, TCGA-HNSCC provides gender information for 500 patients, age information for 499 patients, HPV status information for 481 patients, T stage information for 444 patients, N stage information for 407 patients, and M stage information for 187 patients, pathological stage information for 432 patients, grade information for 481 patients and radiation therapy information for 429 patients. Furthermore, 500 patients' overall survival status information, 499 patients' overall survival time information, progression free survival time and status, 474 patients' disease specific survival time and status, 125 patients' disease specific survival time and status were provided. All data are included in the study. According to the publication guidelines, datasets can be used for publication without restrictions or limitations.

GEO Database
RNA-seq data of SLC13A4 in GSE2837, GSE27020 and GSE41613 were obtained from Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) [23] to study the relationship between SLC13A4 expression level and prognosis of HNSCC patients. Normalized expression matrix files and sequencing platform annotations of the gene sets were downloaded.

Oncomine Database Analysis
Oncomine database (https://www.oncomine.org/resource/main. html) [24] was used for analyzing SLC13A4 expression in various cancer types. The cut-off value of p value and fold change were 0. 01 and 2, respectively.

TIMER Database Analysis
TIMER database (https://cistrome.shinyapps.io/timer/) [25] contained 10,897 samples, including 32 different kinds of cancer types covered by the TCGA dataset. We used the database to evaluate the SLC13A4 expression in tumors of various types and its correlation with T cell exhaustion markers.

The Human Protein Atlas
The human protein atlas (HPA) contains the detailed RNA sequences of 37 major normal tissue types and immunohistochemical results of tissue microarray involving 44 different tissue types. HPA (https://www. proteinatlas.org) [26] has the data about the expression and localization of human proteins in tissues and organs. All the data are openly accessible. We analyzed the expression of SLC13A4 protein in normal oral tissues and oral squamous cell carcinoma by on the basis of immunohistochemistry image in HPA. The antibody used is HPA048582 (SIGMA-ALDRICH).

Survival Analysis
Survival analysis was performed by using survival and survminer packages in R. 499 patients had detailed survival time data and progression-free survival time data, 474 patients had disease specific survival time data, and 124 patients had disease free survival time data in TCGA. 97 Patients had overall survival time data in GSE41613. 28 Patients had recurrence-free survival time data in GSE2837. 109 Patients had disease specific survival time data in GSE27020. All the time data were employed for survival analysis. The median value of survival time from TCGA database and mean value of survival time from GEO database was used as the cut-off value for the survival analysis. The survival curve was plotted by using Kaplan-Meier method, and log rank was used for significance assessment. A p < 0.05 was considered statistically significant.

Immune Infiltration Analysis
To estimate the abundance profile of the tumor-infiltrating immune cells (TIICs) in HNSCC, CIBERSORT (http:// cibersort.stanford.edu/) [26] computational method was applied for analyzing mRNA sequence data of all tumor samples. Based on the validated leukocyte gene signature matrix (termed LM22), the proportion of infiltration of 22 kinds of immune cells in each sample was statistically analyzed. In each sample, the ratio of all TIICs was equal to 1.

Barplot, Corrplot, Vioplot and Scatter Plot
Barplot and corrplot were visualized by using R software with package corrplot, while vioplot was produced with packages BiocManager and vioplot. Except for package BiocManager, scatter plotting was performed by employing R software with packages ggplot2, ggpubr and ggExtra.

Gene Set Enrichment Analysis
GSEA is supported by the Broad Institute Website (http:// software.broadinstitute.org/gsea/index.jsp) [27]. We got HNSCC RNA-sequencing data from TCGA, and performed gene set enrichment analysis by using GSEA software (V4.0.2). We downloaded kegg, hallmark and immune signatures gene sets from Molecular Signatures Database to investigate the possible biological functions of SLC13A4. Gene sets with nominal p value < 0.05 and FDR ≤ 25% were considered statistically significant.

Statistical Analysis
All statistical analyses were conducted by using SPSS 26.0 statistical software (SPSS, Chicago, IL), GraphPad prism 7.0, R X64 4.0.0 and Bioconductor packages (https://www. bioconductor.org/). Wilcoxon test was used to compare the different proportions of TIICs between high SLC13A4 expression group and low expression group. Fisher test was employed to analyze the relationship between SLC13A4 mRNA expression and clinical characteristics. Univariate Cox regression models were utilized to study the prognostic value of clinicopathological parameters for HNSCC. The multivariate Cox regression models were used to assess the association between SLC13A4 mRNA expression and survival adjusted for age, gender, T stage, N stage, pathological stage, grade, HPV status and radiation therapy. Overall survival (OS), disease free survival (DFS), progression free survival (PFS), disease specific survival (DSS) and recurrence free survival (RFS) curves were calculated by Kaplan-Meier method. GraphPad Prism 7.0 was applied for statistical analysis of clinical information. Normal distribution data were subjected to parametric test, and non-normal distribution data to Mann-Whitney test. Significant difference was set at p < 0.05.

Expression of SLC13A4 in Various Human Cancer Types
Using Oncomine database, we analyzed SLC13A4 gene expression in different tumors and the adjacent tissues. The results are shown in Figure 1A. SLC13A4 expression was higher in lymphoma tissues than in adjacent tissues, and was lower in esophageal cancer, HNSCC and sarcoma. To further evaluate SLC13A4 expression in various tumors, we used Timer database to analyze the transcriptional sequence data of 22 TCGA tumors. The results are shown in Figures 1B,C. SLC13A4 expression was lower in breast cancer, head and neck squamous cell carcinoma, kidney chromophobe, thyroid cancer, and higher in renal clear cell carcinoma, hepatocellular carcinoma, cholangiocarcinoma, and rectal adenocarcinoma. In addition, immunohistochemical analyses available from the HPA are shown in Figure 1D, SLC13A4 expression was lower in tumor tissues when compared to normal tissues.

Correlation Between SLC13A4 Expression and Clinical Features in HNSCC Patients
We obtained the data on transcription sequences and clinical characteristics of HNSCC patients from the TCGA database, and divided HNSCC patients into the high-expression group and the low-expression group in terms of the median mRNA expression of SLC13A4. The results are shown in Table 1. HPV status and overall survival in HNSCC patients were significantly correlated with SLC13A4 mRNA expression. HPV-positivity rate (17.7%) was higher in high SLC13A4 expression group than in its low-expression counterpart (10.9%). For the overall survival status, the proportion of deaths in HNSCC patients with high SLC13A4 expression (34.4%) was significantly lower than that with low SLC13A4 expression (52.4%).
We grouped the clinical characteristics of the patients and compared the differences in mRNA expression of SLC13A4 in each group. The results are shown in Figure 2 and Table 1. SLC13A4 expression showed statistically significant differences among T stage, N stage and pathological grades, exhibiting an overall tendency that the more malignant the tumor, the lower the mRNA expression of SLC13A4. The mRNA expression level of SLC13A4 at T1 stage was significantly higher than those at T2 stage (p 0.0342) and T3 stage (p 0.0175). SLC13A4 expression level at N0 stage was significantly higher than that at N2 stage (p 0.0022). SLC13A4 expression at grade 1 was significantly higher than those at grade 2 (p 0.0287) and grade 3 (p 0.003). Since there were only two samples that were in grade 4, no statistical comparison was made.

Univariate and Multivariate Cox Regression Analysis of Factors Affecting Survival in HNSCC Patients
To assess the risk factors influencing the prognosis of HNSCC patients, Cox regression models were used to analyze the clinical data of the patients. Univariate analysis revealed that SLC13A4 expression level, T Stage, N Stage and pathological stage are high risk factors affecting overall survival, disease specific survival and progression free survival in patients with HNSCC. In addition, HPV negative status is also a higher risk factor affecting both overall survival and progression free survival, and radiation therapy is a risk factor only affecting overall survival ( Table 2). Moreover, after adjusted for age, gender, T Stage, N Stage, pathological stage, grade, HPV status and radiation therapy,  the relationship between SLC13A4 mRNA expression and overall survival, disease specific survival and progression free survival were still significant ( Tables 3-5).

Prognostic Value of SLC13A4 in HNSCC
We theorize that SLC13A4 expression is related to the prognosis of cancer patients. We used TCGA data to evaluate the effects of SLC13A4 expression on the total survival time, disease free survival time, disease specific survival time, and progression-free survival time of HNSCC patients. According to the median mRNA value of SLC13A4 expression, tumor samples were divided into a high expression group and a low expression group. Kaplan-Meier plot analysis was performed and results are shown in Figure 3. Compared with the low SLC13A4 expression group, the patients in the high-expression group had longer overall survival time, progression-free survival time and disease specific survival time (p < 0.001), and there was no significant difference in disease free survival time (p 0.217). GEO validation cohort results are shown in Supplementary Figure S1. Compared with the low SLC13A4 expression group, the patients in the high-expression group had longer overall survival time (GSE41613, p 0.010) and recurrence-free survival time (GSE2837, p 0.011). Data from GSE27020 showed that there was no statistical difference in disease free survival time between SLC13A4 low expression group and SLC13A4 high expression group (p 0.092). The findings suggest that SLC13A4 expression level is closely related to clinical characteristics and patients' survival time, and is potentially of value in the assessment of the prognosis of HNSCC. M2 macrophages ( p < 0.001)were significantly higher in the lowexpression group than in the high-expression group. Correlation analysis showed that the expression of SLC13A4 was negatively correlated with monocytes, M1 macrophages, M2 macrophages, resting CD4+ memory T cells, resting NK cells and activated NK cells, but was positively correlated with neutrophils, plasma cells, T follicular helper cells, gamma delta T cells, regulatory T cells and naive B cells. Furthermore, the expression of SLC13A4 in HNSCC is negatively correlated with a variety of T cell exhaustion markers, including LAG3 ( p < 0.001), TIM-3( p < 0.001) and CTLA-4 ( p < 0.001), indicating its potential mechanism in regulating tumor immunity (Supplementary Table S1). These findings showed that SLC13A4 expression was related to the proportion of TIICs and distribution of TIICs in HNSCC, which further demonstrated that SLC13A4 expression affected the immune activity of HNSCC.

Correlation Between SLC13A4 Expression and Modulation of Immunity and Tumor Pathways
In order to further understand the potential mechanism of SLC13A4 affecting the prognosis of HNSCC, patients were divided into a high-and a low-expression groups according to the median expression level of SLC13A4 for GSEA enrichment analysis. The results are shown in Figure 6. The genes in low SLC13A4 expression group were principally involved in immunity-related activities, viral diseases, typical tumor pathways and metabolism, including antigen processing and presentation, regulatory T cells, interleukin-family members, viral myocarditis, epithelial mesenchymal transition (EMT), glycosaminoglycan biosynthesis chondroitin sulfate, among others. In the high SLC13A4 expression group, mainly the metabolism-related pathways were enriched, including metabolisms of linoleic acid, retinol, arachidonic acid and alpha linolenic acid. However, there were no enrichment of genome sets involving immune-related activities and typical tumor pathways in the high SLC13A4 expression group, which further suggested that low SLC13A4 expression might be implicated in the process of cancer development and immune activities.

DISCUSSION
SLC13A4 is a sulfate transporter and is also responsible for the anion transport of selenium and chromium [29]. At present, most studies focused on its roles in severe fetal abnormalities [11], nerve system development [15], and so on. However, sulfate metabolism disorders can cause a variety of diseases and might serve as a predictor for tumor diagnosis [18]. Squamous cell carcinoma is a cancer that arises from particular cells known as squamous cells. Squamous cells are found in the outer layer of skin and in the mucous membranes, which are the moist tissues that form the lining layer of body cavities, such as the airways and intestines. HNSCC is the most common tumor of the head and neck regions, developing in the mucous membrane of the mouth, nose and throat. This study showed that the expression of SLC13A4 was intimately related to the clinical features and prognosis of HNSCC. In HNSCC patients, the expression of SLC13A4 decreased with the increase of tumor pathological grade and clinical stage. The expression level of SLC13A4 in HPV-positive patients was higher than in HPVnegative patients, and the finding was consistent with the phenomenon that the survival rate of HPV-positive patients was higher than that of HPV-negative patients. In this study, only one patient was in M1 stage, so no further analysis was made. Kaplan Meier analysis exhibited that the high expressing SLC13A4 group had longer overall survival time, progression free survival time and disease specific survival time as compared with the low SLC13A4 expression group. On the basis of the aforementioned findings, we are led to believe that that SLC13A4 expression can serve as a predictor of tumor prognosis.
The infiltration of immune cells is one of the characteristic host immune responses to tumor cells and is strongly associated with the development and progression of cancer. Our results also indicated that SLC13A4 expression was significantly correlated with the infiltration of immune cells into HNSCC. SLC13A4 expression was negatively correlated with monocytes, M1 macrophages, M2 macrophages, resting CD4+ memory T cells, resting NK cells and activated NK cells, but was positively correlated with neutrophils, plasma cells, T follicular helper cells, gamma delta T cells, regulatory T cells and naive B cells. These results further confirmed that SLC13A4 level affected the immune activity of TME. SLC13A4 was closely related to the counts of macrophages, M1 and M2, which suggests that the expression of SLC13A4 might be a regulator of macrophage activation and polarization. SLC13A4 was positively correlated with the expression of Tregs, TFH and γδT cells, suggesting that SLC13A4 might regulate the T cell function of HNSCC. Our analysis also showed that SLC13A4 was associated with the degree of immune infiltration  and the number of immune cells in HNSCC. It has been reported that the low ratio of M0 macrophages and high expression of Tregs contribute to the favorable prognosis of OS and DFS in HNSCC patients, which is consistent with our conclusion [30]. These findings suggest that SLC13A4 plays an important role in lymphocyte recruitment and control of immune infiltration in HNSCC.
Our GSEA enrichment analysis demonstrated that the genes of in low SLC13A4 expression group were mainly enriched in immunity, viral diseases, typical tumor pathways and metabolism, while genes in high SLC13A4 high expression group were mainly enriched in metabolismrelated activities, such as metabolisms of linoleic acid, retinol, arachidonic acid and alpha linolenic acid. However, there was no genomic enrichment involving immunity and typical tumor pathways, suggesting that SLC13A4 might act as an indicator of TME transformation from metabolic dominance to immune dominance. The change of infiltrating immune cell components might be triggered by the growth of tumor tissue, and also affect the development of tumors, such as the tumor immune phenotype could predict their prognosis in patients with squamous cell carcinoma of the oral tongue and with squamous cell carcinoma of the head and neck [31,32]. The change of tumor immune microenvironment provides a new target for studying the development and progression of HNSCC. Immunotherapies aim to enhance the activity of the

Pathology & Oncology Research
November 2021 | Volume 27 | Article 1609967 immune system, thereby eradicating cancer cells. Targeting immune checkpoints is a strategy extensively used in cancer therapy. The expression of SLC13A4 in HNSCC is negatively correlated with T cell exhaustion markers, including LAG3, TIM-3 and CTLA-4. Immune checkpoint inhibitors can block these T cell exhaustion markers signals and have proven their good anti-tumor effects, reverse depletion of T cells, and restore tumor microenvironment and tumor-infiltrating T cell functions [33][34][35]. These suggest that SLC13A4 might reverse T cell immunosuppression and inhibit tumor immune escape.
In conclusion, this study demonstrated that SLC13A4 might serve as a biomarker of HNSCC. Our study helps clarify the role of SLC13A4 in the prognosis of HNSCC patients and the transformation from metabolic dominance to immune one in TME, and pave the way to future application of the biomarker in the diagnosis and prognostic prediction of HNSCC [28].

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.