Identification of SLC2A3 as a prognostic indicator correlated with the NF-κB/EMT axis and immune response in head and neck squamous cell carcinoma

ABSTRACT

SLC2A3 is an important member of the glucose transporter superfamily. It has been recently suggested that upregulation of SLC2A3 is associated with poor survival and acts as a prognostic marker in a variety of tumors. Unfortunately, the prognostic role of SLC2A3 in head and neck squamous cell carcinoma (HNSC) is less known. In the present study, we analyzed SLC2A3 expression in HNSC and its correlation with prognosis using TCGA and GEO databases. The results showed that SLC2A3 mRNA expression was higher in HNSC compared with adjacent normal tissues, which was validated with our 9 pairs of HNSC specimens. Moreover, high SLC2A3 expression predicted poor prognosis in HNSC patients. Mechanistically, GSEA revealed that high expression of SLC2A3 was enriched in epithelial-mesenchymal transition (EMT) and NF-κB signaling. In HNSC cell lines, SLC2A3 knockdown inhibited cell proliferation and migration. In addition, NF-κB P65 and EMT-related gene expression was suppressed upon SLC2A3 knockdown, indicating that SLC2A3 may play a preeminent role in the progression of HNSC through the NF-κB/EMT axis. Meanwhile, the expression of SLC2A3 was negatively correlated with immune cells, suggesting that SLC2A3 may be involved in the immune response in HNSC. The correlation between SLC2A3 expression and drug sensitivity was further assessed. In conclusion, our study demonstrated that SLC2A3 could predict the prognosis of HNSC patients and mediate the progression of HNSC via the NF-κB/EMT axis and immune responses.

KEYWORDS:

• SLC2A3
• head and neck squamous cell carcinoma (HNSC)
• epithelial–mesenchymal transition (EMT)
• NF-κB, prognosis

Introduction

Head and neck squamous cell carcinoma (HNSC) accounts for more than 90% of head and neck cancers, and the recurrence rate and distant metastasis are relatively high [1]. Unfortunately, HNSC usually lacks apparent early clinical symptoms, resulting in late-stage diagnosis for most patients [2]. There are more than 600,000 new cases and 380,000 HNSC-related deaths worldwide every year [3]. The current 5-year survival rate for HNSC is lower than 65% [4]. HNSCs originate from mucosal squamous cells, including tumors occurring in the oral cavity, sinus, pharynx and larynx. Numerous genetic and epigenetic aberrations have been described in different individual types, although the molecular mechanisms are not fully understood. Studies have shown that many phenotypic and molecular characteristics of epithelial cells are shared in squamous cell carcinoma [5,6]. Therefore, the molecular characteristics and prognostic markers that could be applied in HNSC need to be further explored.

SLC2A3 belongs to the SLC2A family, consisting of 14 members designated SLC2A1 through SLC2A14. SLC2A (also known as GLUT) genes are the major facilitator superfamily of membrane glucose transporters, which play a central role in glucose uptake along a concentration gradient [7]. In the brain, SLC2A3 is the most abundant glucose transporter and plays a critical role in energy metabolism [8]. The upregulation of SLC2A3 has also been reported in numerous cancer types and acts as a poor outcome indicator in patient-derived prostate xenograft tumors [9], bladder cancer [10], papillary thyroid carcinoma [11], colorectal carcinoma [12] and gastric cancer [13]. However, the role of SLC2A3 in HNSC has not been well clarified.

In this study, we delineated the expression of SLC2A3 in HNSC through The Cancer Genome Atlas (TCGA) database and the role of SLC2A3 in prognosis. Moreover, we validated SLC2A3 expression in our HNSC samples and investigated its impact on the proliferation and migration of HNSC cell lines. In addition, SLC2A3 expression-related signaling pathways, the immune response and drug susceptibility were further explored. Overall, our study demonstrated that upregulation of SLC2A3 may be a potential prognostic biomarker and therapeutic target for HNSC.

Materials and methods

Data collection

We downloaded RNA-sequencing data for 546 HNSC cases, including 56,536 genes and corresponding clinical information, from The Cancer Genome Atlas (TCGA) (http://cancergenome.nih.gov). RNA sequencing data were converted into transcripts per million reads (TPM) and normalized to log2 (TPM+1). Moreover, 270 tumor arrays and 31,330 genes from GSE65858 were applied for prognosis prediction.

Expression of SLC2A3 in HNSC

The TIMER online database (https://cistrome.shinyapps.io/timer) was used to evaluate the expression of SLC2A3 and SLC2A1 across different tumors in TCGA data sets. Then, we used the R package “ggplot2” in R studio to explore the difference in SLC2A3 expression between HNSC tissues and normal samples from the TCGA database. The clinical data for HNSC on the TCGA platform and GEO (www.ncbi.nlm.nih.gov/geo; GSE65858, n = 270) data sets were used to study the relationship between the clinicopathological characteristics and the expression level of SLC2A3, which were classified into high and low groups based on the median value of the SLC2A3 transcript. Using the UALCAN database (http://ualcan.path.uab.edu/index.html), we analyzed the expression levels of the SLC2A3 gene in different tumor stages of HNSC and normal tissues in The Cancer Genome Atlas (TCGA) HNSC data set.

Prognosis and survival analysis

To evaluate the potential independent prognostic value of SLC2A3 in HNSC patients, we obtained clinicopathological data from the TCGA platform, including gender, age, T/N/M stage, clinical stage and lymphatic vascular invasion. The R package “glmnet” was used for univariate and multivariate Cox regression analyses. The parameters with P < 0.1 were retained for multivariate analysis and adjusted according to other standards (N/M stage, lymphovascular invasion). Statistical significance was set as p < 0.05. The R packages “ggrisk,” “survival,” and “survminer” were utilized for survival analysis. Receiver operating characteristic curve (ROC) analysis was carried out using the R package “pROC” and “ggplot2.” The patients with a relative value above or below the median were considered the high or low group for each gene expression level. Kaplan-Meier survival analysis with the log-rank test was then used to investigate the overall survival (OS) difference between the abovementioned groups.

Enrichment analysis of the SLC2A3 gene coexpression network in HNSC

The TCGA-HNSC data sets were analyzed using the R package “stat” to study the coexpressed genes related to SLC2A3 expression. Statistical correlations were obtained by calculating Pearson’s correlation coefficients, and the R package “ggplot2” was used to draw a heatmap. Gene Ontology (GO) function and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of coexpressed genes were performed by the “clusterProfifiler” package (version: 3.18.0) of R software, and visual analysis of data was performed by the “ggplot2” software package.

Gene Set Enrichment Analysis

To explore the possible pathway through which SLC2A3 may be involved in the regulation of HNSC, we divided the TCGA HNSC samples into two groups based on the median expression of SLC2A3 and performed gene set enrichment analysis (GSEA) with the differentially expressed genes (DEGs) between the low and high SLC2A3 expression groups. GSEA software (v4.1.0) was used to conduct the enrichment analysis. The MSigDB database V7.2 was used to obtain the hallmark gene set. False discovery rate (FDR)<0.25 and p.adjust<0.05 were considered statistically significant.

Immune cell analysis

The CIBERSORTx website (https://cibersortx.stanford.edu/) was utilized to analyze the relative abundance of 22 immune cell subtypes among TCGA HNSC. CIBERSORTx used a deconvolution algorithm to calculate the frequency of immune cells in each sample based on the expression levels and the matrix of gene expression of 22 immune cell-associated genes [14]. We used the somatic copy number alteration (SCNA) module in the TIMER platform to analyze the relationship between SLC2A3 copy number variation (CNV) and immune cells. Kaplan‒Meier curve analysis was further performed to study the survival status based on different levels of immune cells.

Drug sensitivity analysis

We downloaded NCI-60 compound activity data from CellMiner™ (http://discover.nci.nih.gov/cellminer/). The data were analyzed by the R packages “impute,” “limma,” “ggplot2,” and “ggpubr” to identify the correlation between SLC2A3 expression and drug IC50.

Clinical samples

A total of 9 HNSC and paired adjacent normal tissues were obtained from our Department of ENT, the Affiliated Hospital of Qingdao University, Qingdao, Shandong, China. The study was approved by the Research Ethics Committee of our hospital. The histopathology of the specimens was determined by pathologists in the pathology department. The patients had not received any chemotherapy and/or radiotherapy before surgery.

Cell culture and shRNA construction

FaDu and CNE-2 cells were obtained from Shanghai Cell Bank (Shanghai, China) and cultured in high glucose DMEM with 10% fetal bovine serum. The two SLC2A3 shRNA target sequences were 5-CTTGGTCTTTGTAGCCTTCTT-3 and 5-AGTAGCTAAGTCGGTTGAAAT-3. The target sequence for the negative control was 5-CGCTGAGTACTTCGAAATGTC-3 against firefly luciferase. The primers were annealed and ligated into the tet-PLKO vector digested with AgeI and EcoRI. The shRNA plasmids were transfected into cells, which were then treated with 2 μg/ml doxycycline.

Western blot

The cells were lysed with 2% SDS buffer. After transfer, the membrane was incubated with primary antibodies against SLC2A3 (R24453, Zen Bio), E-cad (bs -10,009 R, Bioss), N-cad (bs-1172 R, Bioss), Vimentin (bs-0756 R, Bioss), Snai1 (bs-1371 R, Bioss), and NF-κB P65 (AP0123, ABclonal). GAPDH (E-AB-20059, Elabscience) and α-tubulin (AC007, ABclonal) were used as loading controls. The bands were quantified with ImageJ software and normalized to the loading control.

Quantitative PCR

Total RNA was extracted with RNA-easy Isolation Reagent (#R701, Vazyme). After quantification, the same amount of RNA was reverse transcribed with an RT Reagent Kit (RR037A, Takara). Quantitative PCR was performed using SYBR qPCR reagent (Q111, Vazyme). GAPDH for cells and 18S rRNA for tissues were used as reference genes. The following primers were used: SLC2A3-F: GTCGGACTCTTCGTCAACCG, SLC2A3-R: AGAGGCCAATAACCAAGCGA; CXCR4-F: CTCCTCTTTGTCATCACGCTTCC, CXCR4-R: GGATGAGGACACTGCTGTAGAG; CCR1-F: GCATCATCATTTGGGCCCTG, CCR1-R: AAAGTGAAGGCTGCAGGTGT; CCR8-F: TGGCTGTTGTCCATGCCGTGTA, CCR8-R: TGGGATGGTAGCCATAATGGCG; CCL2-F: GCAATCAATGCCCCAGTCAC, CCL2-R: CTTTGGGACACTTGCTGCTG; CCL13-F:CAGCAGGTGTCCCCAGAAG, CCL13-R: TTCAGGGTGTGAGCTTTCCG; 18 s rRNA-F: CTGGATACCGCAGCTAGGAA, 18 s rRNA-R: GAATTTCACCTCTAGCGGCG; GAPDH-F: TCGACAGTCAGCCGCATCTT, GAPDH-R: GAGTTAAAAGCAGCCCTGGTG.

Cell Counting Kit-8 (CCK8) assay

A total of 5000 cells were seeded into 96-well plates with six technical replicates. The next day, the cells were treated with 2 μg/ml doxycycline after transfection. For drug sensitivity validation, cells were treated with 750 nM obatoclax (Beyotime, SC4381) and 10 μM triciribine (Beyotime, SF2721) for 24 h after 48 h of transfection. At the indicated time point, the mixed solution of 100 μl DMEM and 10 μl CCK8 (A311, Vazyme, china) was added into each well and the cells were further incubated at 37°C for 30–60 minutes. The absorbance was measured at 450 nm after a half-hour incubation at 37°C.

In vitro scratch assay

The cells were seeded into 6-well plates and transfected with plasmids. At 24 hours after transfection, the image was taken as time point zero after scratching. Then, the medium was changed to serum-free medium with 2 μg/ml doxycycline. Images were taken at the indicated time points. The wound area was measured with ImageJ software. The wound recovery percentage was normalized to time point zero.

Statistical analysis

The data are presented as the mean ± SD. The experiments were repeated three times. Student’s t test and Mann – Whitney U test were used to determine the differences between the two groups. R (version 4.03) and Rstudio software were used for statistical analysis. Univariate and multivariate Cox regression analyses were used to identify the correlation between SLC2A3 expression levels and clinical characteristics. Survival analysis was assessed by the Kaplan – Meier method and the log-rank test. p < 0.05 was set as statistically significant.

Results

Elevated expression of SLC2A3 predicts poor survival in HNSC

First, we explored the expression of SLC2A3 mRNA in human tumors in the TCGA database. Compared to the normal tissues, SLC2A3 was found to be notably higher in tumor tissues of cholangiocarcinoma (CHOL), stomach adenocarcinoma (STAD), kidney renal clear cell carcinoma (KIRC) and HNSC but markedly decreased in bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), colon adenocarcinoma (COAD), glioblastoma multiforme (GBM), kidney chromophobe (KICH), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC) and uterine corpus endometrial carcinoma (UCEC) (Figure 1a). In addition, the expression of another important family member, SLC2A1, was also analyzed (Supplementary Figure S1).

Figure 1. SLC2A3 expression and correlation with survival status in HNSC. (a) SLC2A3 mRNA expression levels in pan-cancer data from TCGA. (b) Comparison of SLC2A3 mRNA expression levels between all normal tissues and HNSC samples from TCGA. (c) SLC2A3 mRNA expression in matched HNSC samples from TCGA. (d) the relationship between SLC2A3 mRNA expression and tumor stages of HNSC patients in the UALCAN database. The correlation between SLC2A3 expression and (e) disease-specific survival, (f) overall survival (OS) in TCGA and (G) OS in GSE65858. (h) ROC curve analysis of SLC2A3 diagnosis. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001.

To study the expression level and prognostic impact of SLC2A3 in HNSC, we next analyzed the HNSC RNA-seq data in the TCGA database. The results showed that SLC2A3 expression was significantly upregulated in all HNSC tissues compared with normal tissues (Figure 1b). In line with the results, SLC2A3 was also increased in paired HNSC tissues (Figure 1c). Further exploration with the UALCAN database revealed that SLC2A3 expression was higher in late-stage HNSC than in early-stage HNSC (Figure 1d). We next investigated the relationship between SLC2A3 expression and survival status in TCGA and GSE65858. HNSC patients with high SLC2A3 expression had poor disease-specific survival (Figure 1e) and overall survival (Figures 1f, 1g). These results suggest that SLC2A3 overexpression is associated with HNSC progression. In addition, ROC curve demonstrated that SLC2A3 could distinguish malignancy from adjacent normal tissues with AUC 0.727 (95% CI: 0.635–0.819) (Figure 1h), indicating that SLC2A3 could be used as a diagnostic marker in HNSC.

Association of SLC2A3 expression and clinical features of HNSC

Next, we deciphered the relationship between SLC2A3 expression and the clinical characteristics of HNSC. The results showed that SLC2A3 expression was significantly correlated with T classification, M classification and clinical stage (Table S1). N type, M type, and lymphatic infiltration were prognostic risk factors in univariate Cox regression analysis (Figure 2a). In multivariate Cox regression analysis, all three prognostic risk factors were identified as independent significant predictors (Figure 2b). Moreover, SLC2A3 expression was compared among different clinical TNM stages and pathological degrees in HNSC patients. The results showed that higher SLC2A3 expression predicted poor survival in the same TNM stages and clinical stages and the differences look more apparent for the T stage stratification (Figure 2c-f).

Figure 2. SLC2A3 expression was associated with the clinicopathological features of HNSC based on TCGA. (a) Univariate and (b) multivariate Cox analyses of SLC2A3 expression and clinicopathological variables. Kaplan – Meier analysis of overall survival for patients showed significant differences between high and low SLC2A3 expression groups with (c) T stage (d), N stage (e), M stage and (f) clinical stages.

Enrichment analysis of genes coexpressed with SLC2A3 in HNSC

The coexpressed genes related to SLC2A3 expression in the TCGA-HNSC database were obtained by the deseq2 package of R software. A heatmap showed the top 50 genes positively or negatively correlated with SLC2A3 expression (Figures 3a, 3b). We further performed GO and KEGG analyses for the top upregulated and downregulated genes in Figure 3a,b. The results demonstrated that the top 50 upregulated genes were mainly related to ECM and focal adhesion and the PI3K-AKT pathway (supplementary figure S2A). Interestingly, the top 50 downregulated genes were strongly related to epidermal cell differentiation, keratinization and cornification (supplementary figure S2B). The results suggested that SLC2A3 may play an important role in cell proliferation, differentiation and migration. In addition, we performed OS univariate analysis of the top 50 up/downregulated genes associated with SLC2A3 expression, which is shown in Supplementary Table 1 and Table 2. GO function and KEGG pathway enrichment of the top 200 coexpressed genes positively correlated with the expression of SLC2A3 were performed by the “clusterProfifiler” package. With the threshold of p.adj<0.05 and qvalue<0.2, 318 biological processes (GO-BP), 20 cellular components (GO-CC), 29 molecular functions (GO-MF) and 10 KEGGs were obtained. KEGG pathway analysis demonstrated that SLC2A3 was mainly related to the PI3K−AKT signaling pathway (Figure 3c). The GO analysis data showed that most of the genes were associated with extracellular matrix structure (Figure 3d-f). The top 10 are displayed in the bubble chart, among which the PI3K-AKT signaling pathway and extracellular matrix structure were the most prominent features (Figure 3c-f).

Figure 3. Enrichment analysis of the SLC2A3 gene coexpression network in HNSC. Heatmaps show the top 50 coexpressed genes (a) positively and (b) negatively correlated with SLC2A3 expression in HNSC. Enrichment analysis of (c) KEGG and (d – f) GO for SLC2A3 coexpression genes.

SLC2A3 was positively correlated with the EMT signaling pathway in HNSC

To further elucidate the underlying mechanisms of SLC2A3 in the progression of HNSC, the TCGA HNSC data sets were divided into two groups based on the median expression level of SLC2A3. GSEA data set genes were sorted according to their differential expression between the two groups. The results demonstrated that the epithelial-mesenchymal transition (EMT) pathway was significantly enriched in the activated pathways, which contained the most abundant genes among the top 3 pathways (Figure 4a). We further performed GSEA of the EMT pathway subset for DEGs, which ranked in order as shown in the lower panel of Figure 4b. GSEA demonstrated that the EMT pathway was enriched in the upregulated genes (NES = 2.376) (Figure 4b). Furthermore, a significant positive correlation was observed between the expression of SLC2A3 and EMT-related gene signatures, including N-cad, vimentin, MMP9, Snai1, Snai2, Twist1 and Twist2 (Figure 4d-j). There was a negative correlation between SLC2A3 and E-cad, although it was not statistically significant (Figure 4c). Our results suggested that SLC2A3 may promote the development of HNSC by regulating the EMT pathway.

Figure 4. SLC2A3 expression was correlated with EMT gene signatures in TCGA-HNSC. GSEA of DEGs between the low and high SLC2A3 expression groups from TCGA-HNSC for (a) activated pathways and (b) EMT gene signature. NES: normalized enrichment score. FDR: false discovery rate. (c-j) Correlation between SLC2A3 expression and EMT-related genes.

SLC2A3 was positively correlated with NF-κB in HNSC

NF-κB activation is essential for the EMT pathway in many cancers [15,16]. GSEA revealed that NF-κB was enriched in the upregulated genes in HNSC (Figure 5a). SLC2A3 expression was positively correlated with five NF-κB family members, NF-κB1, NF-κB2, REL, RELA and RELB (Figure 5b). Moreover, we screened all five NF-κB family members and only RELB and NF-κB2 gene expression levels were positively correlated with N-cad, MMP9, Snai1 and Vimentin and negatively correlated with E-cad (Figure 5c).

Figure 5. SLC2A3 expression was associated with NF-κB, and NF-κB was associated with EMT gene signatures in TCGA-HNSC. (a) GSEA of DEGs between the low and high SLC2A3 expression groups from TCGA-HNSC for NF-κB signaling. NES: normalized enrichment score. FDR: false discovery rate. (b) Correlation between SLC2A3 and the expression of NF-κB1, NF-κB2, REL, RELA and RELB. (c) Both RELB and NF-κB2 were positively correlated with the expression of N-cad, vimentin, Snai1 and MMP9 and negatively correlated with the expression of E-cad.

SLC2A3 regulated HNSC progression through the immune response

SLC2A3 may affect the immune response [17]. GSEA revealed that IL6/JAK/STAT3 signaling and interferon-γ response pathways were enriched in the upregulated genes in HNSC (Figures 6a, 6b). In addition, the enrichment analysis of immune characteristics was performed for all the CXCR, CCR and CCL family members in HNSC samples and HNSC cell lines. The SLC2A3 expression was only positively correlated with the immune response genes CXCR4, CCR1, CCR8, CCL2 and CCL13 (Figure 6c-g). The qPCR results validated that SLC2A3 knockdown inhibited the mRNA expression of CXCR4, CCR1, CCR8, CCL2 and CCL13 (Figure 6h).

Figure 6. SLC2A3 expression was correlated with IL6/JAK/STAT3 signaling and interferon-γ hallmark gene sets. GSEA of DEGs between the low and high SLC2A3 expression groups from TCGA-HNSC for the (a) IL6/JAK/STAT3 signaling pathway and (b) interferon-$\gamma$ response pathway. NES: normalized enrichment score. FDR: false discovery rate. (c-g) the positive correlation between SLC2A3 and the expression of CXCR4, CCR1, CCR, CCL2 and CCL13 based on TCGA-HNSC. (h) SLC2A3 knockdown inhibited the mRNA expression of CXCR4, CCR1, CCR8, CCL2, and CCL13. Each value represents the mean ± SD of three separate experiments. *P < 0.05; **P < 0.01.

Relationship between SLC2A3 expression and tumor immune cells

Tumor immune cells play an important role in the occurrence and development of tumors [18–21]. To explore the relationship between SLC2A3 and tumor immune cells, we analyzed the correlation between SLC2A3 expression and HNSC immune cells by CIBERSORTx. Compared with the low SLC2A3 expression group, in the high SLC2A3 expression group, the proportions of resting memory CD4 T cells, macrophages (M0 and M2), activated mast cells and neutrophils were increased, whereas decreased proportions of CD8 T cells, activated NK cells, monocytes, resting DCs and resting mast cells were observed (Figure 7a). Meanwhile, compared with the diploid/normal group, Arm-level gain of SLC2A3 resulted in decreased infiltration levels of B cells, CD8+ T cells, neutrophils and dendritic cells (Figure 7b). However, TIMER analysis demonstrated that high infiltration levels of B cells, CD8+ T cells and dendritic cells still predicted better OS of HNSC patients with high SLC2A3 expression (P < 0.05) (Figure 7c).

Figure 7. Correlation between SLC2A3 expression and tumor immune infiltrating cells. (a) Infiltration of 22 immune infiltrating cells between HNSC patients with low and high SLC2A3 expression obtained by CIBERSORTx. (b) Comparison of immune infiltrating cell levels in HNSC with different SLC2A3 CNVs. (c) Clinical survival outcome of HNSC patients with high B cell, CD8+ T-cell and myeloid dendritic cell infiltration. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001.

Drug sensitivity analysis

Based on the CellMiner database, we next investigated the potential correlation between drug sensitivity and SLC2A3 expression in HNSC. As a result, SLC2A3 expression was negatively associated with the IC50 of obatoclax, triciribine phosphat, carfilzomib, BX − 912, SCH−1473759, acetalax, BML−277, CH − 5132799, 7−hydroxystaurosporine and palbociclib (Figure 8a, b, D-H, J, K) and positively associated with XL − 147 and trametinib (Figures 8c, 8I), which indicated that SLC2A3 expression levels could be an effective prognostic indicator for HNSC patients in chemotherapy. We further validated the effect of obatoclax and triciribine with a CCK8 assay. The FaDu and CNE-2 cells showed slightly better viability upon SLC2A3 knockdown (Figure 8M).

Figure 8. Scatter plot of the relationship between SLC2A3 expression and drug sensitivity. SLC2A3 expression was (a, b, d-h, j-l) negatively or (c, i) positively correlated with the IC50 of drugs approved by the FDA. (M) Effect of obatoclax and triciribine on cell viability upon SLC2A3 knockdown. *P < 0.05.

Validation of SLC2A3 expression in HNSC samples

To validate the SLC2A3 expression pattern in HNSC, we collected 9 pairs of HNSC originating in the nasopharynx, hypopharynx, larynx, oropharynx and skull base (Table S1). The results showed that SLC2A3 mRNA expression was significantly elevated in cancer samples compared to adjacent normal tissues (Figure 9a).

Figure 9. SLC2A3 knockdown inhibited proliferation and migration and NF-κB/EMT-related gene expression in FaDu and CNE-2 cell lines. (a) SLC2A3 mRNA expression in adjacent normal tissues (ANT) and HNSC samples. (b) Cell proliferation and (c) migration were suppressed upon SLC2A3 knockdown. (d) SLC2A3 knockdown inhibited N-cad, Vimentin, and Snai1 expression but increased E-cad expression. (e) NF-κB P65 expression was decreased upon SLC2A3 knockdown. Each value represents the mean ± SD of three separate experiments. * p < 0.05, # p < 0.01.

SLC2A3 silencing attenuated EMT signaling in FaDu and CNE-2 cells

Next, we uncovered the underlying mechanisms of SLC2A3 in FaDu and CNE-2 cells. Cell proliferation (Figure 9b) and migration (Figure 9c) were inhibited upon SLC2A3 knockdown. Our results of data mining in HNSC revealed that SLC2A3 expression was strongly associated with the NF-κB pathway and EMT signaling. We found that SLC2A3 knockdown resulted in enhanced expression of E-cad and downregulation of N-cad, Vimentin and Snai1 in both cell lines (Figure 9d). It has been reported that the NF-κB pathway is involved in the induction and maintenance of EMT. Our GSEA also indicated that NF-κB expression is positively correlated with the EMT signature. In HNSC cell lines, we observed that NF-κB P65 protein was suppressed by SLC2A3 knockdown. Together, our results suggested that SLC2A3 may regulate EMT signaling via the NF-κB pathway (Figure 9e).

Discussion

HNSC is one of the most common malignant cancers worldwide and arises from the epithelium in the head and neck, leading to as many as 50% of deaths [22]. However, there is a lack of effective biomarkers for the prognosis of HNSC. Our results suggested that SLC2A3 could be used as a biomarker in HNSC. Studies have demonstrated that SLC2A3 is a primary neuronal glucose transporter. It also plays an important role in glucose absorption in sperm, preimplantation embryos and cultured leukocytes [23]. In the context of tumors, SLC2A3 expression is increased in various types of tumors, including colon cancer [24], gastric cancer [25], breast cancer [26], and lung cancer [27]. It has been well documented that tumor cells increase the uptake and metabolism of glucose to generate energy even under normal oxygen conditions. The low K_M and high V_MAX of SLC2A3 on glucose uptake suggested a preeminent role for SLC2A3 in energy production [28]. We searched all the SLC2A members and found that SLC2A3 could be a novel prognostic biomarker in HNSC.

In our study, we found that the expression level of SLC2A3 in HNSC tissues was significantly higher than that in normal tissues through bioinformatic analysis. Moreover, we also observed the enhanced expression of SLC2A3 in our clinical samples. A large number of studies have shown that SLC2A3 promotes the growth and proliferation of tumor cells in different cancers [24,27,29]. Our results showed that cell proliferation and migration were suppressed upon SLC2A3 knockdown in FaDu and CNE-2 cells. Recent studies have indicated the predictive value of the SLC2A family. A 3-gene prognostic index including BRCA1, QKI, and SLC2A1 indicated a high risk of recurrence and death in early-stage lung adenocarcinoma patients [30], and SLC2A1 and SLC2A3 were significant indicators of poor prognosis in stage I non-small cell lung cancer (NSCLC) [31]. In laryngeal squamous cell carcinoma, patients with positive expression of a 3-gene signature including SLC2A1, vimentin and N-cad have much shorter survival times than those with negative expression [32], whereas another study demonstrated that SLC2A3 but not SLC2A1 acted as a significant marker of poor prognosis [33]. In mammary epithelial cells, SLC2A1 knockdown led to increased Vimentin and decreased E-cad [34]. A recent meta-analysis revealed that SLC2A1 was associated with poorer OS among HNSC patients [35]. The expression of SLC2A1 and SLC2A3 was significantly elevated in biopsies of primary HNSC compared to the corresponding normal gingiva [36]. Botha et al. revealed that only SLC2A1 and SLC2A3 have potential value in the prognosis of oral squamous cell carcinoma [37]. High expression of SLC2A3 was correlated with poor survival in different cancer types, including gastric cancer [25], glioma [38], colorectal cancer (CRC) [24] and NSCLC [27]. Unfortunately, the predictive role of SLC2A3 in HNSC is less known. Here, our results demonstrated that SLC2A3 predicted poor survival in HNSC. The SLC2A3 level was correlated with increased clinical stage and lymph node metastasis, indicating the potential prognostic value of SLC2A3 in HNSC. The correlation with stratification of T stage is more apparent, indicating SLC2A3 expression may represent an early event in HNSC progression. Nevertheless, no obvious prediction value of SLC2A1 was observed. This may be attributed to different sample sources in these studies.

Current studies on the role of SLC2A3 have mainly focused on glucose metabolism [24,39,40]. In this study, the coexpressed genes associated with SLC2A3 in HNSC were analyzed by the R software package. The results showed that the PI3K/AKT signaling and extracellular matrix pathways were enriched in the KEGG and GO analyses. It has been well documented that the PI3K/AKT signaling pathway is involved in GLUT1/4 (SLC2A1/4) translocation to the plasma membrane in normal tissues and tumors [41,42]. However, the link between PI3K/AKT signaling and GLUT3 (SLC2A3) trafficking is less addressed. AKT could regulate GLUT3 surface localization during Mycobacterium tuberculosis infection in macrophages [43]. In addition, PI3K/AKT signaling could promote GLUT3 expression in the SH-SY5Y neuroblastoma cell line [44]. Conversely, in nonmalignant human breast cells, GLUT3 overexpression activated AKT signaling, and SLC2A3 knockdown led to AKT inactivation [45]. Interestingly, the extracellular matrix pathway ranked the top in all three GO subontologies (BP, CC, MF). It has been reported that extracellular matrix organization is one of the top three enriched GO terms in SLC2A3 compared to SLC2A1-overexpressing glioblastoma cells [46]. On the other hand, the breakdown of extracellular matrix hyaluronan (HA) promoted SLC2A1 expression, resulting in increased tumor cell migration via upregulation of glycolysis. SLC2A3 was also increased in LiSa-2 cells upon HAase treatment, although it did not show similar consistency in some other cell lines, indicating cell type dependency [47]. It seems that the regulation between the extracellular matrix and glucose transporters might be reciprocal. However, the link between extracellular matrix remodeling and glucose metabolism needs to be further explored in HNSC.

We also found that the EMT pathway was correlated with SLC2A3 in HNSC. EMT is an important pathological process that converts epithelial cells to mesenchymal-like cells to increase invasiveness, metastasis and drug sensitivity and therefore promotes tumor progression. In laryngeal squamous cell carcinoma, SLC2A1 expression was positively correlated with Vimentin and N-cad expression but negatively correlated with E-cad expression [32]. SLC2A3 expression has been suggested to be correlated with EMT signatures in CRC and NSCLC [27,48]. In HNSC, we also found that SLC2A3 expression was highly associated with EMT signaling. Furthermore, by western blotting, we validated that upon SLC2A3 silencing, the expression of EMT-related genes, including N-cad, vimentin and Snai1, decreased, whereas that of E-cad increased. Interestingly, the important EMT gene ZEB1 could directly bind to the SLC2A3 gene promoter to activate transcription [27]. Our results demonstrated that the strong positive correlation between NF-κB and SLC2A3 expression and NF-κB P65 protein expression was inhibited by SLC2A3 knockdown in FaDu and CNE-2 cells. It has been reported that the NF-κB pathway mediates EMT-related gene expression, including Snail, E-cad, ZEB-1/2, and Twist1 [49–52], resulting in EMT phenotypes. Together, our study indicated that the SLC2A3/NF-κB/EMT axis is involved in HNSC progression. In addition, it has also been shown in the top 3 pathways that TGF-β and hypoxia are correlated with SLC2A3 expression. Hypoxia is a crucial factor in tumor development. Hypoxia could enhance SLC2A3 expression in normal and tumor cells. As a hypoxia marker, SLC2A3 expression was associated with poor prognosis in esophageal adenocarcinoma [53], intestinal tumors [54,55], clear cell renal cell carcinoma [56], oral squamous cell carcinoma [57] and hepatocellular carcinoma [58]. TGF-β signaling pathway plays a pivotal role in inducing EMT. In lung cancer cells TGF-β could upregulate SLC2A3 expression through the EMT gene ZEB1 [59]. However, the link between TGF-β signaling and SLC2A3 expression is much less addressed and needs to be further investigated.

SLC2A3, as one of the combined biomarkers of the high-risk group, was significantly correlated with immune cells and lower expression of immune checkpoint molecules in esophageal adenocarcinoma, HNSC and CRC [53,60,61]. In CRC, the expression of SLC2A3 was positively correlated with CD4+ and CD8+ T-cell infiltration [48]. SLC2A3 is involved in inflammatory gene expression and M1 tumor-associated macrophages partially by regulating CXCL8 expression in triple-negative breast cancer [26], yet SLC2A3 mediates the M2 subtype transition of macrophage infiltration in the gastric cancer microenvironment [13]. In oral squamous cell carcinoma, a high level of CD8+ T-cell infiltration is associated with an improved patient outcome [62]. Our results showed that CD8+ T cells and resting dendritic cells were decreased in high SLC2A3 expression group. Moreover, the arm-level gain of SLC2A3 was related to low immune cell infiltration, including B cells, CD8+ T cells, neutrophils and dendritic cells. We screened the correlation of SLC2A3 expression with all the CXCR, CCR and CCL family members. A significant correlation was observed only with CXCR4, CCR1, CCR8, CCL2 and CCL13. Interestingly, CXCR4 has been identified as an important immune suppressor, and inhibition of CXCR4 relieves this suppression, resulting in enhanced T and NK cell infiltration [63,64]. These results suggest that SLC2A3 may inhibit the immune response of the HNSC tumor microenvironment, which may explain the worse outcome of HNSC patients with high SLC2A3 expression. However, the underlying mechanisms need to be further investigated. It should be noted that the high infiltration of B cells, CD8+ T cells and dendritic cells could still indicate better outcomes in HNSC patients with high SLC2A3 expression.

In conclusion, our study demonstrated that the upregulation of SLC2A3 was strongly associated with clinical prognosis in HNSC. SLC2A3 plays an important role in promoting HNSC cell proliferation and migration. The SLC2A3-induced phenotype may be attributed to the activation of the NF-κB/EMT pathway and immune response. Together, these findings suggest that SLC2A3 could be used as a biomarker for HNSC prognosis and as a treatment target.

Abbreviation

CNV: copy number variation, DEGs: differentially expressed genes, EMT: epithelial-mesenchymal transition, GEO: Gene Expression Omnibus, GO: Gene ontology, HNSC: head and neck squamous cell carcinoma, GSEA: gene set enrichment analysis, KEGG: Kyoto Encyclopedia of Genes and Genomes, OS: overall survival, ROC: receiver operating characteristic, SCNA: somatic copy number alteration, TCGA: The Cancer Genome Atlas

Author contributions

FC and JZ wrote the paper and conceived and designed the experiments; TF, PJ and YH analyzed the data; LW, SY, XY, LY, ZX, RW, BX and XD collected and provided the sample for this study. All authors approved the final version of the manuscript.

Ethical approval

This study was performed in line with the principles of the Declaration of Helsinki. Approval was granted by the ethics committee of the Affiliated Hospital of Qingdao University.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The data that support the findings of this study are available on request from the corresponding author.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/19336950.2023.2208928.

Additional information

Funding

This study was supported by the National Natural Science Foundation of China (81672662), Natural Science Foundation of Shandong Province (ZR2022MH073), Qingdao Medical and Health Science and Technology Program (2021-WJZD184) and Clinical Medicine +X Scientific Research Project of Affiliated Hospital of Qingdao University (QDFY+X202101012).