SLC20A1 is a prospective prognostic and therapy response predictive biomarker in head and neck squamous cell carcinoma

Data collection

The Cancer Genome Atlas (TCGA) database was used as the source from which data on SLC20A1 expression, prognostic information and clinical features of HNSCC cases was extracted (https://xena.ucsc.edu/). Normalized gene transcript data from Fragments Per Kilobase of transcript per Million mapped reads (FPKM) were involved in this study. The clinicopathological features of the included specimens, like gender, age, grade, node (N) and tumor (T) statuses, and stage were used in the analyses. Another 24 paired tissues were stored in our laboratory. The Ethics Committee of the Affiliated Lihuili Hospital has granted authorization for the conduct of this research, under the approval number 2022SL415-01.

Expression difference of the SLC20A1 and its association with DNA methylation in HNSCC

Perl software was used to extract the SLC20A1 mRNA expression level in HNSCC from the HTSeq data. R software (version 3.6.0; limma package) is applied for analysis of the differential expression of the SLC20A1 in HNSCC samples by comparing it to normal tissues. In addition, the methylation of these cg sites (cg05422897, cg05672265, cg06121808, cg10762132, cg10898730, cg23886783, cg26703507) in SLC20A1 gene’s promoter regions was obtained by downloading DNA methylation data from Illumina's Human Methylation 450K data. Annotation data used for cg sites were downloaded from Illumina's help desk (https://support.illumina.com). Utilizing the Pearson correlation with the R package corrplot, our team analyzed the correlation between DNA methylation and SLC20A1 mRNA expression.

Prognosis value of SLC20A1 in HNSCC

Using the optimal threshold for SLC20A1 expression identified by the iterative algorithm of the survival package, we constructed the overall survival (OS) curves using the Kaplan-Meier method for a more precise analysis. A risk assessment model for prognosis prediction was generated to bind the expression level with survival outcome of each patient. TCGA datasets were used to perform all analyses. The independent prognostic value of SLC20A1 in HNSCC was assessed using uni- and multivariate Cox regression analyses. Two R packages, survminer and survival, were used for all these analyses and the ggplot was used for the results represented by forest plot.

Nomogram construction for prognosis prediction

A nomogram was constructed to provide a quantitative approach for predicting survival probability in HNSCC, which comprised the SLC20A1 clinical, pathological, and expression factors, including gender, age, grade, stage, T status and N status, by using R package rms. In the nomogram, we used the OS rates at 1, 3, and 5 years as endpoints. Area under the curve (AUC) and receiver operating characteristic (ROC) curve were developed to evaluate these rates’ predictive ability. Furthermore, calibration plots were drawn to assess the nomogram accuracy for internal validation [19].

Prediction of therapeutic sensitivity in HNSCC with different expression of SLC20A1

To evaluate the chemoresistance in HNSCC, the 50% inhibiting concentration (IC50) value of the 3 medicine which was mainly used in chemotherapy for HNSCC was inferred using the pRRophetic algorithm. The TIDE score uses gene expression patterns to predict clinical outcome for immune checkpoint inhibitors (ICIs) [20]. A lower TIDE score is associated with a greater immunotherapy response.

Expression of SLC20A1 and the association with tumor infiltrating immune cells (TIICs)

The relationship between SLC20A1 and 22 subsets of tumor-infiltrating lymphocytes (mostly neutrophils, macrophages, dendritic cells, CD8+ and CD4+ T cells, and B cells) were evaluated using the CIBERSORT database (https://cibersort.stanford.edu/). For this purpose, we utilized the Tumor Immune Estimation Resource (TIMER) database (https://cistrome.shinyapps.io/timer/) to investigate TIICs and SLC20A1 expression in depth.

Functional enrichment analysis

Enrichment studies of the thematic expression of SLC20A1 were performed to identify the possible molecular processes and roles of SLC20A1 in HNSCC. Two further studies, using the Kyoto Encyclopedia of Genes and Genomes (KEGG) and the Gene Ontology (GO), were conducted. The R programs cluster profiler and enrichplot were applied to run KEGG pathway and GO functional enrichment studies. To investigate potential signaling pathways in which SLC20A1 participates, a gene set enrichment analysis (GSEA) was run. Remarkably variable expression levels were defined as having a false discovery rate (FDR) of < 0.05 and a | log2FC | ≥ 1.

Cell culture and transfection

TU686 (laryngeal squamous cell carcinoma cell line) and CAL-27 (tongue squamous cell carcinoma cell line) human HNSCC cells were employed in this investigation. HNSCC cells were cultured in 0.05 mg/mL gentamycin (Schering-Plough Europe, Belgium) and RPMI 1640 with 10% FBS. Saturated humidity with 5% CO₂, and 37° C was used to cultivate cells. GeneChem created SLC20A1-specific siRNAs (si-SLC20A1#1/2/3) and si-NC (China). These plasmids were introduced into CAL-27 and TU686 cell lines using Lipotransfectamine 3000 (Thermo Fisher Scientific, USA).

CCK-8 experiment and colony formation assay

Cell lines CAL-27 and TU686 were grown in 96-well plates at a density of 5×10³ cells/well for 0, 24, 48, 72, and 96 hours. After two hours in the dark at 37° C, cells were counted using a cell counting kit 8 (CCK-8) (Dojindo, Japan). Cell absorbance was measured at a wavelength of 450 nm. Six-well plates were seeded with CAL-27 and TU686 cells (at a density of 500 cells per well), and the media were replaced every four days. Colonies were fixed in 4% paraformaldehyde for 30 minutes and dyed with 0.1% crystal violet (Sigma-Aldrich, USA) for 15 minutes after being incubated for two weeks. The colony count was performed using Image J program.

Wound healing assay

A six-well plate was used for the cellular seeding. After the cells had adhered together to create a monolayer, a scratch was made in it using the tip of a plastic pipette with 200 L. Cells should be washed with PBS three times. After discarding the whole DMEM, 2 mL of DMEM devoid of serum was applied to each well. Using an IX71 inverted microscope (Olympus, Japan), we photographed the drawn distance of each hole at the stated periods (0 h, 48 h) and analyzed the photos using Image J. There were three rounds of each test.

Transwell assay

By a 24-well Transwell chamber, a Transwell experiment was done to test how well the cells could invade. TU686 and CAL-27 cell lines in media without serum were put in top chambers that had been covered with 2% Matrigel (BD Biosciences, USA). In the bottom chambers, medium with 20% FBS was added. After two days of incubation, 4% paraformaldehyde was used to fix the cells that had moved to the bottom of the membrane. The cells were fixed, and then crystal violet was used to stain them for 15 minutes. With the use of Image-Pro Insight software and an Olympus IX71 microscope from Japan, the findings were seen and recorded (Olympus, Japan).

Quantitative real-time polymerase chain reaction (qRT-PCR)

Using TRIZOL (Invitrogen, USA), total RNA was extracted from 24 HNSCC and normal samples, and then cDNA was synthesized using a reverse transcription kit (Novoprotein, China). The qRT-PCR reactions were performed using a LightCycler 480 (Roche, USA) instrument and NovoStart SYBR qPCR (Novoprotein, China), as per the manufacturers' instructions. During each of the 35 cycles, the reaction was subjected to the following conditions: heating, denaturation, annealing, and extension at 95° C, 60° C, 72° C, and 95° C for 30 s, 10 s, 30 s, and 30 s, respectively. The SLC20A1 relative expression was determined using the two-fold comparative Ct technique (2Ct). Each sample has three independent tests conducted on it. These were the primer sequences that were provided: SLC20A1, F: 5’-TGGCAACGCTGATTACCAGT-3’; R: 5’-CAGCCCTTGAGTCGAGTTGT-3’. GAPDH, F: 5’-TCAAGATCATCAGCAATGCC-3’; R: 5’-CGATACCAAAGTTGTCATGGA-3’.

Western blotting analysis

Tissue was collected, lysed, and protein content was determined using a BCA protein assay kit (Beyotime, China) per the manufacturer's instructions. Using a 10% sodium dodecyl sulfatepolyacrylamide gel, we separated the protein that had been extracted using the improved RIPA buffer (Beyotime, China). The membrane was then transferred onto a PVDF membrane (Millipore, USA) at a continuous current of 240 mA for 2 hours, followed by blocking in 5% skim milk for 2 hours and incubation with anti-SLC20A1 primary antibody (Proteintech, China) at 4 degrees Celsius for 1 night. Following a 3-time TBST wash, the membrane was incubated with a secondary antibody from Proteintech (China) for an hour. Eventually, chemiluminescence (Advansta, USA) was used to identify and visualize the target protein with Image Lab software.

Statistical analysis

Documents for GESA preparation, DNA methylation, and HTSeq FPKM were parsed using Perl software version 5.32 to glean the relevant data. Differential expression, nomogram generation, prognostic value assessment, and Pearson correlation were all performed using R 4.0.3 program with specialized packages. Appropriate testing methods were selected for comparison based on the type of samples, including the T-test, Chi-square test, and Wilcoxon test. A significance level of 0.05 was used.

Data availability

The data that supported the findings of this study are openly available from The Cancer Genome Atlas (TCGA) program (https://portal.gdc.cancer.gov/).

HNSCC prognosis and SLC20A1 expression

We analyzed 383 cases to see whether or not there was a correlation between SLC20A1 expression and patient features. It was found that patients that expressed a lot of SLC20A1 had a higher-than-average T stage, as shown in Table 1 (P < 0.05). Our findings demonstrated the significantly higher expression levels of SLC20A1 mRNA in 499 HNSCC samples compared to 44 adjacent normal ones, obtained from the TCGA database (Figure 1A, P<0.05). qRT-PCR analyses of 24 pair-matched samples stored in our laboratory consistently showed a remarkable increase of SLC20A1 mRNA in HNSCC tissues (Figure 1B, P<0.05). Western blot of five paired tissues showed that SLC20A1 protein was highly expressed in HNSCC tissues (Figure 1C). Furthermore, SLC20A1 methylation status in promoter cg sites was also investigated, as was the correlation with gene expression. From Pearson correlation, we could obtain that the abnormal DNA methylation might lead to the high expression of SLC20A1 (Figure 1D).

Figure 1. Exaggerated SLC20A1 expression levels in HNSCC. (A, B) SLC20A1 mRNA expression between the HNSCC specimens and nearby sound tissues from TCGA database and qRT-PCR analysis. (C) SLC20A1 protein expression in the HNSCC specimens and nearby sound ones from Western blotting analysis. (D) DNA methylation modification associated with SLC20A1 expression in HNSCC. *, P<0.05.

Additional investigation of SLC20A1's relevance to survival in HNSCC was performed using Kaplan-Meier analysis. We determined the threshold for separating patients into those with high or low SLC20A1 expression based on SLC20A1 expression level, survival status, and OS time. The findings showed that individuals having a high SLC20A1 expression had a higher risk score, and a shorter OS time (Figure 2A, 2B). In addition, stage, age, and SLC20A1 expression were shown to be independent predictive variables for HNSCC in both uni- and multivariate analyses using Cox regression (Figure 2C, 2D).

Figure 2. HNSCC prognosis via SLC20A1. (A) Kaplan-Meier analysis of poor survival outcomes with SLC20A1 overexpression in HNSCC from TCGA database. (B) Risk score with survival time in TCGA database. (C, D) Forest plots representing the uni- and multivariate analysis of the significantly marked prognostics.

Silencing SLC20A1 inhibited cell invasion, migration, and proliferation in HNSCC

In our research, we used the CAL-27 and TU686 cell lines to examine SLC20A1's role in HNSCC cells. Transfection of siRNA for nucleotides cofactors (NC) and SLC20A1 were performed in CAL-27 and TU686 cell lines, respectively. Cell proliferation was reduced in both cell lines treated with si-SLC20A1 compared to controls, as measured by cell counting tests and colony formation (Figure 3, P<0.05). The wound healing experiments were utilized to measure the migratory potential of HNSCC cells, and the findings indicated a substantially smaller wound area in the controls compared to cases in the si-SLC20A1 groups (Figure 4A, P<0.05). Furthermore, Transwell assay findings showed that SLC20A1-down-expressed groups were less invasive than controls (Figure 4B, P<0.05). Finally, we hypothesized that knocking down SLC20A1 in HNSCC cells would reduce their ability to proliferate, migrate, and invade.

Figure 3. Silencing SLC20A1 inhibited cell proliferation in TU686 and CAL-27 cells lines. (A) CCK8 and (B) colony formation experiments showed that silencing SLC20A1 inhibited cell proliferation. *, P<0.05. **, P<0.01. ***, P<0.001.

Figure 4. Silencing SLC20A1 inhibited migration and invasion in HNSCC. (A) Wound-healing experiments showed that silencing SLC20A1 inhibited cell migration. (B) Silencing SLC20A1 inhibited the ability to invade CAL-27 and TU686 cell lines. *, P<0.05. **, P<0.01. ***, P<0.001.

Nomogram development

A nomogram has been established based on SLC20A1 expression and clinical characteristics (Figure 5A), as a quantitative approach for HNSCC prognosis prediction. Calibration curve indicates good agreements between the observed outcome and predicted probability (Figure 5B). For further validation of the prognostic power of combined clinical features and SLC20A1 expression, ROC curve has been created and the AUC for survival rates at 1, 3, and 5 years were 0.676, 0.750 and 0.716, with moderate prediction accuracy (Figure 5C), indicating the significant ability of the developed monogram for predicting HNSCC prognosis.

Figure 5. A nomogram based on SLC20A1 expression and clinical features. (A) Nomogram to predict the overall survival (OS) of the 1^st, 3^rd, and 5^th years for HNSCC cases. (B) Calibration curve of the nomogram prognostic model. (C) ROC curve of the nomogram for predicting the overall survival (OS) of the 1^st, 3^rd, and 5^th years.

The association between SLC20A1 overexpression with chemotherapy and immunotherapy response in patients with HNSCC

As shown in Figure 6A, SLC20A1 expression was remarkably associated with cisplatin, gemcitabine, and paclitaxel. Moreover, SLC20A1 overexpression was more sensitive to cisplatin and gemcitabine, while low expression of SLC20A1 were more sensitive to paclitaxel.

Figure 6. Chemotherapy and immunotherapy response were associated with SLC20A1 expression in HNSCC. (A) Chemotherapeutic responses with differential SLC20A1 expression in HNSCC. (B) Correlation between SLC20A1 expression and immunotherapy responses. *, P<0.05. **, P<0.01. ***, P<0.001. IC50, 50% inhibiting concentration.

The results demonstrated that SLC20A1 overexpression was associated with more immune exclusion and dysfunction, and a lower TIDE score, indicating that patients with a higher SLC20A1 expression level are less sensitive to immunotherapy (Figure 6B). Specifically, we looked at how SLC20A1 expression relates to TIICs in HNSCC. SLC20A1 expression variation was strongly correlated with a wide range of immune cell types (Supplementary Figure 1). In our analysis of the TIMER and the CIBERSORT databases, we found that high SLC20A1 expression was negatively correlated with resting mast cells, dendritic cells, resting NK cells, T cells regulatory (Tregs), and resting CD8+ T cells, and positively correlated with Eosinophils, activated mast cells, macrophages M0, resting NK cells, and CD4+ memory T cells (Figure 7A, 7B). This led us to hypothesize that SLC20A1 expression has a role in modulating the immunological microenvironment in HNSCC, which is important for further studies.

Figure 7. Association between SLC20A1 expression and TIICs in HNSCC. (A) Abundance of TIICs infiltration between low and high SLC20A1 expression in HNSCC from CIBERSORT database. (B) Correlation between SLC20A1 expression and TIICs from TIMER database.

SLC20A1-related underlying molecular mechanism in HNSCC

A protein-protein interaction (PPI) network was constructed using the GENEMANIA database, followed by a Pearson correlation analysis based on gene expression data from TCGA. This analysis revealed a strong correlation between SLC20A1 and SLC20A2, especially in shared protein structural domains (Supplementary Figure 2). Analyses of GO functions and KEGG pathways were used to uncover a possible molecular mechanism involving SLC20A1 in HNSCC. The findings of the GO analysis revealed the correlation between cellular component, biological process and molecular function, and SLC20A1 expression (Supplementary Figure 3A, 3B). Meanwhile, the results from KEGG enrichment analysis revealed a remarkably enrichment of SLC20A1 in immune process and several pathways, including “IL−17 signaling pathway”, “cytokine−cytokine receptor interaction”, “TNF signaling pathway”, “NOD-like receptor signaling pathway”, “NF-κB signaling pathway”, “JAK-STAT signaling pathway”, “Chemokine signaling pathway” (Supplementary Figure 3C, 3D). According to the initial screening of pathways, the TNF signaling pathway has garnered our attention. Our bioinformatics analysis revealed significant overexpression of TNF and TNFR2 in head and neck squamous cell carcinoma (Supplementary Figure 4). Furthermore, there was a notable positive correlation between the expression of SLC20A1 and both TNF and TNFR2. In vitro experiments with the head and neck squamous carcinoma cell line TU686 showed that knocking out SLC20A1 significantly reduced TNF and TNFR2 expression levels compared to the control group (Supplementary Figure 5). Hence, we hypothesize that SLC20A1 may be involved in the regulation of the TNF/TNFR2 signaling pathway.

In addition, GSEA interpretation indicated that SLC20A1 overexpression had markedly positive correlations with numerous enrichment pathways (Supplementary Figure 6A, 6B), including “bladder cancer”, “colorectal cancer”, “renal cell cancer”, “small cell lung cancer”, “cell cycle”, “natural killer cell mediated cytotoxicity”, “P53 signaling pathway” and “MAKP signaling pathway”, In contrast, SLC20A1 low expression was significantly related with several metabolism and biochemical processes, such as “alpha linolenic acid metabolism”, “drug metabolism cytochrome P450”, “arachidonic acid metabolism” and “oxidative phosphorylation” (Supplementary Figure 6C). Accordingly, SLC20A1 might act as a novel biomarker in tumor progression and immune environment in HNSCC.