Non-coding RNA-mediated high expression of SFXN3 as a prognostic biomarker associated with paclitaxel resistance and immunosuppressive microenvironment in head and neck cancer

Article type: Research Article

Keywords: sideroflexin 3, paclitaxel resistance, immune infiltration, targeted therapy, head and neck cancer

Authors: Kailin Chen, Sha Gong, Xueliang Fang, Qian Li, Mingliang Ye, Junyan Li, Shengyan Huang, Yuheng Zhao, Na Liu, Yingqin Li, Jun Ma

Affiliations: Department of Radiation Oncology, State Key Laboratory of Oncology in South China, Collaborative Innovation Center of Cancer Medicine, Guangdong Key Laboratory of Nasopharyngeal Carcinoma Diagnosis and Therapy, Sun Yat-sen University Cancer Center, Guangzhou, China; Department of Experimental Research, State Key Laboratory of Oncology in South China, Collaborative Innovation Center of Cancer Medicine, Guangdong Key Laboratory of Nasopharyngeal Carcinoma Diagnosis and Therapy, Sun Yat-sen University Cancer Center, Guangzhou, China

License: Copyright © 2022 Chen, Gong, Fang, Li, Ye, Li, Huang, Zhao, Liu, Li and Ma CC BY 4.0 This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Article links: DOI: 10.3389/fimmu.2022.920136  |  PubMed: 36159813  |  PMC: PMC9493355

Relevance: Moderate: mentioned 3+ times in text

Full text: PDF (10.0 MB)

Introduction

Head and neck cancer is the seventh most common cancer globally, and squamous cell carcinoma is the most common histological type (ref. 1, ref. 2). The majority of patients (>50%) with head and neck squamous cell carcinoma (HNSC) present with a locally advanced stage or distant metastasis at diagnosis (ref. 3). Despite multimodal therapeutic interventions, the outcome of patients with advanced-stage HNSC remains unsatisfactory: 40%–60% of these patients ultimately experience recurrence or local progression, with poor prognosis (ref. 4). Therapeutic resistance, especially chemoresistance, is the leading cause of poor prognosis (ref. 2, ref. 5). Therefore, identifying promising biomarkers for stratified chemotherapy and prognosis, and novel therapeutic targets for the effective management of HNSC, are urgently required.

Sideroflexin 3 (SFXN3) is a member of the sideroflexin family, of which there are five homologs of these mitochondrial proteins in humans. The expression of SFXN homologs varies across tissues (ref. 6). The SFXN3 protein level is relatively high in the mature brain and neurons (ref. 7–ref. 9). Consistent with these finding, SFXN3 was reported to be involved in neurodegenerative diseases, including Parkinson’s disease, Alzheimer’s disease, and other neurological dysfunction diseases (ref. 10–ref. 12). It is widely accepted that many cancer cells are highly dependent on serine metabolism for their growth and proliferation (ref. 10, ref. 13, ref. 14). SFXN3, one of the main mitochondrial serine transporters during one-carbon metabolism, might play a crucial role in the development of cancer. However, few studies have focused on the role of SFXN3 in cancer. Murase et al. reported that serum autoantibodies to SFXN3 might be a novel tumor marker for oral squamous cell carcinoma (ref. 15). These findings prompted us to conduct a comprehensive study of the expression, prognosis, and possible mechanism of SFXN3 in HNSC. Moreover, although SFXN3 was reported to be expressed in the stromal fibroblasts around cancer nests (ref. 15), the correlation of SFXN3 with drug resistance and tumor immune infiltration in HNSC has not been determined.

In this study, we analyzed the expression and prognostic value of SFXN3 in pan-cancer. Next, we determined the effect of SFXN3 expression on drug sensitivity. Importantly, the roles of SFXN3 in HNSC immune infiltration were further explored. Mechanistically, non-coding RNAs (ncRNAs) related to the regulation of SFXN3, such as long non-coding RNAs (lncRNAs) and microRNAs (miRNAs), were also investigated in HNSC. Taken together, our findings suggested that the ncRNA-mediated upregulation of SFXN3 was associated with poor prognosis, paclitaxel resistance, and the immunosuppressive microenvironment in HNSC.

Materials and methods

TCGA data preprocessing

The normalized pan-cancer dataset, TCGA Pan-Cancer (PANCAN, N = 10,535, G = 60,499), was downloaded from the UCSC (https://xenabrowser.net/) database. The samples of common cancer types were used for further analysis, including glioblastoma multiforme (GBM), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), lung adenocarcinoma (LUAD), colon adenocarcinoma (COAD), breast invasive carcinoma (BRCA), esophageal carcinoma (ESCA), stomach and esophageal carcinoma (STES), kidney renal papillary cell carcinoma (KIRP), stomach adenocarcinoma (STAD), prostate adenocarcinoma (PRAD), uterine corpus endometrial carcinoma (UCEC), head and neck squamous cell carcinoma (HNSC), kidney renal clear cell carcinoma (KIRC), lung squamous cell carcinoma (LUSC), liver hepatocellular carcinoma (LIHC), thyroid carcinoma (THCA), rectum adenocarcinoma (READ), pancreatic adenocarcinoma (PAAD), pheochromocytoma and paraganglioma (PCPG), bladder urothelial carcinoma (BLCA), kidney chromophobe (KICH), and cholangiocarcinoma (CHOL). Then, we extracted the expression data of the ENSG00000107819 (SFXN3) gene in each sample with log2(x+0.001) transformation for each expression value.

In addition, transcriptional profiles of count values and clinicopathological characteristics for HNSC (n = 546) were also obtained from the UCSC (https://xenabrowser.net/) database. Normalization was conducted by the “limma” package. The univariate and multivariate Cox regression analyses for survival in HNSC and CIBERSORTx and xCell analysis were conducted.

Differential expression and survival analysis

The Wilcoxon rank-sum test was used to compare expression values between normal and tumor samples in pan-cancer analysis (ref. 16). Cox regression was used to analyze the relationship between gene expression and prognosis in each tumor. The log-rank test was used to compare two survival curves. The Sangerbox online tool (ref. 17) was used for statistical analysis and graphing. The UALCAN and Gene Expression Profiling Interactive Analysis (GEPIA) databases were used to further confirm these results. UALCAN (http://ualcan.path.uab.edu/index.html) is a web resource providing a comprehensive analysis of cancer transcriptome data, such as that from TCGA (ref. 18). Transcripts per million (TPM) values were employed to estimate the significance of the differences in gene expression levels between groups using a t test. GEPIA (http://gepia.cancer-pku.cn/) is a web tool for analyzing RNA sequencing expression and correlation data of cancer and normal samples from TCGA and the Genotype-Tissue Expression (GTEx) data (ref. 19). SFXN3 and lncRNA expressions on box plots were determined by GEPIA in various types of human cancer. The method used for differential analysis was one-way ANOVA, and the p value cutoff was 0.05. Survival analysis for SFXN3 expression and candidate lncRNAs in cancers was also performed at GEPIA using the log-rank test. The statistical analyses of the online databases in this study were calculated automatically. Median values were employed as group cutoffs, and a log-rank p value < 0.05 was identified as statistically significant.

Tumor immune estimation resource, CIBERSORTx, and xCell analysis

Tumor Immune Estimation Resource (TIMER) is a comprehensive resource for the analysis of tumor-infiltrating immune cells across more than 30 types of cancer (https://cistrome.shinyapps.io/timer/) (ref. 20). In the present study, we compared tumor infiltration levels among tumors with different somatic copy number alterations for SFXN3 and the correlation of SFXN3 expression with tumor immune infiltration using TIMER in the HNSC dataset.

Moreover, we applied the CIBERSORTx computational method to calculate the proportion of 22 subtypes of immune cells in case of HNSC from TCGA database (ref. 21). According to the cell type description in the LM22 signature matrix, 11 major leukocyte types were divided into 22 leukocyte subsets (ref. 22). The classification of immune cells was unified based on the description of comparative data from other databases more conveniently. We analyzed the fractions of the 10 main types of immune cells, including B cells (naïve B cells and memory B cells), plasma cells, CD8 T cells, CD4 T cells (CD4 naïve T cells, CD4 memory resting T cells, CD4 memory activated T cells, T follicular helper cells (Tfhs), and regulatory T cells (Tregs)), gamma delta T cells, natural killer (NK) cells (resting NK cells, activated NK cells), monocytes, macrophages (M0 macrophages, M1 macrophages, M2 macrophages), dendritic cells (resting dendritic cells, activated dendritic cells), and neutrophils. xCell (https://xcell.ucsf.edu/) is another signature-based method for quantifying the infiltration of immune cells in tissues. We used the “xCell” R package to generate immune estimates to validate the CIBERSORT results by R software (version 3.6.4) (ref. 23).

Assessment of the prognostic role of SFXN3 in relation with the immune infiltration

We employed a univariate Cox regression to identify the variables with a close relationship to OS in HNSC patients. The variables in the analysis included SFXN3 expression, the levels of immune cell infiltration by CIBERSORT, clinical data (age, T stage, N stage, and stage), and tumor purity. The high and low groups of immune cell infiltration levels were determined based on the median values. Tumor purity estimates on HNSC were obtained from a previous study (ref. 24). Four methods were used for tumor purity estimates: ESTIMATE, ABSOLUTE, LUMP, and IHC. ESTIMATE applied gene expression profiles of 141 immune genes and 141 stromal genes to evaluate (ref. 25). Somatic copy-number data were used for ABSOLUTE (ref. 26). LUMP is the abbreviation of leukocytes unmethylation for purity; it averages 44 non-methylated immune-specific CpG sites. IHC was estimated by image analysis of hematoxylin and eosin stain slides (ref. 24).

All covariates with a p value less than 0.1 based on univariate analysis were included in the multivariate model. The backward Wald method was used for the multivariate Cox proportional hazards regression model by SPSS 24.0 software (IBM, Armonk, NY, USA).

Drug sensitivity analysis

To clarify the association between SFXN3 expression and drug sensitivity and tolerance, we downloaded transcriptome data and drug sensitivity half-maximal inhibitory concentration (IC50) values of the NCI-60 panel of human cancer cell lines from the CellMiner database (https://discover.nci.nih.gov/cellminer/) (ref. 27). The therapeutic effects of 574 kinds of drugs or compounds in clinical trials and 218 US FDA-approved drugs in the database were selected for further analysis. The effect of SFXN3 expression on the drug sensitivity was analyzed using the “impute” (ref. 28) and “limma” R packages (ref. 29). Missing data for some drugs were imputed using impute. knn from the impute R package. There are 1,275 missing values, 5% of the total. Spearman’s correlation test was used to analyze the relationship between the activity of drugs or compounds and the expression of SFXN3. Linear models were used for analysis of drug sensitivity. The Wilcoxon test was utilized to analyze the difference in the activity of drugs between the high- and low-risk SFXN3 expression groups, and then a box plot was drawn using the “ggplot2” and “ggpubr” functions of R.

To assess the relationship between SFXN3 expression and the sensitivity of commonly used chemotherapy drugs for HNSC, we downloaded the HNSC cell line half-maximal inhibitory concentration (IC50) values of commonly used chemotherapy drugs from the Genomics of Drug Sensitivity in Cancer (GDSC) database (ref. 30). These values can be linked with the SFXN3 expressions of the corresponding cell lines reported in the Cancer Cell Line Encyclopedia (https://sites.broadinstitute.org/ccle/). In the analysis of each drug, the expression of SFNX3 was grouped high and low according to the median value. The Wilcoxon rank-sum test was used to compare IC50 values between the high and low SFXN3 expression groups using the Sangerbox online tool.

Prediction of candidate miRNAs and lncRNAs

The upstream miRNAs of SFXN3 were predicted using eight target gene prediction programs, including DIANA TOOLS (ref. 31), miRDB (ref. 32), miRmap (ref. 33), miWalk (ref. 34), RNA22 (ref. 35), TargetScan (ref. 36), microT-CDS (ref. 37), and miRTarBASE (ref. 38). Only the miRNAs predicted in more than two programs were selected for subsequent analysis. CancerMIRNome (http://bioinfo.jialab-ucr.org/CancerMIRNome/) is a comprehensive database for the interactive analysis of miRNA expression profiles from TCGA projects (ref. 39). The datasets of significant downregulated miRNAs and positive prognostic miRNAs in HNSC were downloaded from CancerMIRNome. Pearson’s correlation analysis was conducted between the expression of hsa-miR-29c-3p and the expression of SFXN3.

The upstream candidate lncRNAs of hsa-miR-29c-3p were predicted using the ENCORI database (https://starbase.sysu.edu.cn/index.php) (ref. 40). Correlation analyses for lncRNA–hsa-miR-29c-3p or lncRNA–SFXN3 and the expression level of lncRNAs in HNSC were performed at ENCORI using Pearson’s correlation coefficient.

Results

Pan-cancer analysis of SFXN3 expression and its prognostic value

To explore the possible roles of SFXN3 in pan-cancer, we first analyzed its expression in 22 common cancer types. Compared with that in the normal samples, the expression of SFXN3 was significantly upregulated in 12 cancer types, namely, COAD, ESCA, STES, KIRP, STAD, HNSC, KIRC, LIHC, THCA, READ, PCPG, and CHOL (all p < 0.01, f1), while the expression of SFXN3 was significantly downregulated in GBM, LUAD, PRAD, UCEC, LUSC, and KICH. Next, the UALACAN and GEPIA databases were used to validate the expression of SFXN3 in pan-cancer. As shown in Figure S1A, the expression of SFXN3 was significantly upregulated in five cancer types by the UALACAN database, including CHOL, COAD, HNSC, KIRC, KIRP, LIHC, READ, THCA, and STAD (all p < 0.001). SFXN3 expression levels in CHOL, HNSC, KIRC, acute myeloid leukemia (LAML), PAAD, PCPG, and sarcoma (SARC) were upregulated compared with those in the corresponding normal controls in the GEPIA database (all p < 0.05, Figure S1B). Taken together, these data demonstrated that SFXN3 was upregulated in CHOL, HNSC, and KIRC, indicating that SFXN3 might act as a crucial regulator in the carcinogenesis of these three cancer types.

Furthermore, Cox survival analysis was performed on SFXN3 expression in pan-cancer. Both overall survival (OS) and disease-specific survival (DSS) analyses were conducted. For OS, a higher expression of SFXN3 in patients was associated with significantly worse prognosis in COAD, HNSC, LIHC, PAAD, and KICH (p = 0.008, 0.01, 0.005, 0.02, and 0.02, respectively). COAD, HNSC, and KICH were consistent with a high expression in the tumor (f1). However, patients with ESCA, STES, KIRP, STAD, KIRC, THCA, READ, PCPG, and CHOL with a higher expression of SFXN3 showed no significance in terms of OS (all p > 0.05).

For DSS, Cox regression analysis found that a higher expression of SFXN3 demonstrated poorer prognosis in COAD and HNSC (p = 0.002 and 0.005, respectively, f1). In addition, survival analysis was also performed in GEPIA, including CESC, CHOL, HNSC, KICH, KIRC, LUSC, PAAD, PCPG, PRAD, and UCEC. A higher expression of SFXN3 in patients was associated with significantly worse prognosis in HNSC for overall survival (p = 0.006) (Figure S1C). For DFS, a higher SFXN3 expression demonstrated poorer prognosis in HNSC (p = 0.037, Figure S1D).

These results demonstrated that SFXN3 expression may serve as a negative and unfavorable prognostic marker for patients with HNSC.

Poor prognosis of patients with HNSC with a high SFXN3 expression was attributed to enriched tumor-infiltrating macrophages

Head and neck cancer is a profoundly immunosuppressive disease with high immune infiltration (ref. 41, ref. 42). Therefore, we evaluated the relationship between SFXN3 expression and the level of immune infiltration in HNSC. As shown in f2, SFXN3 expression was significantly and positively associated with the infiltration level of CD4⁺ T cells (r = 0.304, p < 0.001), macrophages (r = 0.239, p < 0.001), neutrophils (r = 0.257, p < 0.001), and dendritic cells (r = 0.323, p < 0.001), but not B cells (r = −0.04, p = 0.382) and CD8⁺ T cells (r = −0.084, p = 0.067). The CIBERSORT algorithm was used to further confirm these results. The levels of tumor‐infiltrating B cells, CD8⁺ T cells, dendritic cells, gamma delta T cells, macrophages, NK cells, and plasma cells correlated significantly with SFXN3 expression (all p < 0.05, f2). Moreover, we found that the levels of tumor-infiltrating dendritic cells, macrophages, myocytes, and plasma cells also correlated significantly with SFXN3 expression, as calculated using xCell (all p < 0.05, Figure S2)

Taken together, a high SFXN3 expression was associated with significantly higher proportions of tumor-infiltrating macrophages. Correlation analyses between SFXN3 expression and biomarkers of M1 or M2 polarized phenotypes of macrophages were further investigated. SFXN3 expression correlated significantly and positively with biomarkers of immunosuppressive M2 macrophages (CD163: r = 0.355, p < 0.001; the V-set and immunoglobulin domain-containing protein 4 (VSIG4): r = 0.318, p < 0.001; and membrane-spanning 4-domains subfamily A 4A (MS4A4A): r = 0.35, p < 0.001) and correlated negatively with one of the M1 macrophages’ biomarkers (NO synthase 2 (NOS2): r = −0.192, p < 0.001) (f2).

The data demonstrated that a high expression of SFXN3 in HNSC associated with tumor progression and poor prognosis might partly be related to enriched tumor-infiltrating macrophages.

SFXN3 expression is significantly related to the immunosuppressive microenvironment in HNSC

In addition to tumor-associated macrophages (TAMs), especially M2 macrophages, myeloid-derived suppressor cells (MDSCs) and Tregs also contribute to forming an immunosuppressive microenvironment in tumor tissues (ref. 43). Next, we explored the relationship between SFXN3 expression and biomarkers of the main immune-suppressive cells including MDSCs, TAMs, and Tregs in the tumor microenvironment. We conducted a correlation analysis using TIMER, which showed that the immune markers of MDSCs (CD33: r = 0.284, p < 0.001; integrin subunit alpha M (ITGAM): r = 0.181, p < 0.001, fucosyltransferase 4 (FUT4): r = 0.332, p < 0.001), TAMs (C–C motif chemokine ligand 2 (CCL2): r = 0.226, p < 0.001, CD68: r = 0.339, p < 0.001, interleukin (IL)-10: r = 0.354, p < 0.001), and Tregs (forkhead box P3 (FOXP3): r = 0. 258, p < 0.001, C–C motif chemokine receptor 8 (CCR8): r = 0.304, p < 0.001, signal transducer and activator of transcription 5 (STAT5B): r = 0.243, p < 0.001) were significantly and positively associated with the SFXN3 expression level ( f3) (ref. 44, ref. 45). Inhibition of T-cell activation is one of the common causes contributing to macrophage-related tumor immune evasion (ref. 46). We further performed a correlation analysis between SFXN3 expression and inhibitory factors involved in the T-cell activation (ref. 46). The analysis showed that the SFXN3 expression level correlated positively with the expression of genes related to the inhibition of T-cell activation, including CD28 (encoding CD28 molecule) (r =0.209, p < 0.001), CTLA4 [encoding cytotoxic T lymphocyte-associated antigen-4) (r = 0.151, p < 0.001), PDL1 (encoding programmed cell death 1 ligand 1 (also known as CD274)] (r = 0.267, p < 0.001), PDCD1LG2 (encoding programmed cell death 1 ligand 2 (PD-L2)) (r = 0.51, p < 0.001), VSIR (encoding V-set immunoregulatory receptor (also known as C10ORF54); r = 0.128, p = 0.004), and TIGIT (encoding T-cell immunoglobulin and ITIM domain) (r = 0.118, p = 0.007) (f3).

These data indicated that a high SFXN3 expression was associated with an immunosuppressive microenvironment, possibly via inhibition of T-cell activity.

SFXN3 expression was an independent risk factor

Univariate Cox regression analysis revealed that SFNX3 expression [hazard ratio (HR): 1.282 (1.097-1.489), p = 0.002], tumor-infiltrating macrophages [HR: 1.411 (1.077-1.847), p = 0.012], gender [HR: 0.75 (0.564-0.998), p = 0.048], and N stage [HR: 1.37 (1.030-1.822), p = 0.031] were significant predictors of OS in patients with HNSC, which are shown in f4.

All covariates with a p value less than 0.10 based on univariate analysis were included in the multivariate model. Only the following variables remained significant after multivariate adjustment: SFXN3 [HR: 1.247 (1.058-1.471), p = 0.009], gender [HR: 0.731 (0.547-0.977), p = 0.034], and N stage [HR: 1.543 (1.153-2.064), p = 0.004], which are demonstrated in f4.

The values of SFXN3 in drug resistance analysis

To investigate the potential values of SFXN3 expression to predict drug resistance in cancer, we acquired NCI-60 compound activity and corresponding NCI-60 cell line RNA-seq/composite expression from the CellMiner database. Compound activity scores and sensitivity to the corresponding compound correlated positively. We further selected 792 drugs or compounds in the database that are closely related to clinical treatment, including 574 drugs in clinical trials and 218 US FDA-approved drugs, for further analysis. Overall, SFXN3 expression was found to be associated with the increased activity of 42 drugs or compounds and the decreased activity of 20 drugs or compounds (f5).

According to the National Comprehensive Cancer Network (NCCN) guidelines, cisplatin, gemcitabine, paclitaxel, docetaxel, carboplatin, 5-FU, and methotrexate are commonly used chemotherapy drugs for head and neck cancer. Correlation analysis for these drugs showed that the activity of paclitaxel and docetaxel was negatively and significantly related to the expression of SFXN3 (r = −0.265, p = 0.041, and r = −0.371, p = 0.004, respectively; f5). We also found a trend of negative correlation between SFXN3 expression and the activity of cisplatin, carboplatin, or methotrexate, but without a significant difference (all p > 0.05). We also divided patients into high or low SFXN3 expression groups according to the median value of SFXN3 expression and then investigated whether SFXN3 expression has an impact on the activity of the drugs. Patients with a high expression of SFXN3 were found to possess a lower drug activity for paclitaxel and docetaxel (both p < 0.05, Figure S3).

In order to evaluate whether SFXN3 expression has an impact on the activity of the drugs in HNSC, we divided the expression values of SFNX3 in different HNSC cell lines into two high and low groups, which are shown in supplementary material Table S1. HNSC cell lines with a high expression of SFXN3 were confirmed to have higher IC50 values for paclitaxel (p = 0.03, f5). We also found a trend of the correlation between SFXN3 expression and IC50 of docetaxel, but without a significant difference (p = 0.06). Moreover, there was no significant difference between SFXN3 expression and the IC50 values of methotrexate, 5-FU, cisplatin, and gemcitabine (p =0.90, 0.56, 0.70, and 0.90, respectively).

These findings suggested that SFXN3 expression might function as a paclitaxel resistance predictor, which might be closely related to poor prognosis in head and neck cancer.

Prediction regulatory miRNAs of SFXN3

Non-coding RNAs can act as regulators of gene expression, which has been widely investigated. MicroRNAs (miRNAs) are short RNAs that perform their function by binding to the target mRNAs (ref. 47).

To clarify the possible mechanism by which ncRNAs might mediate the regulation of SFXN3, we first predicted miRNAs that could bind to SFXN3. Using eight prediction programs, 568 predicted miRNAs were identified, which are shown in supplementary material Table S2. According to the mechanisms of the miRNA–target interaction, the predicted upstream miRNAs should correlate negatively with SFXN3 expression. In the CancerMIRNome analysis of TCGA database, 29 positive prognostic miRNAs and 131 significantly downregulated miRNAs were found in HNSC (f6). Only hsa-miR-29c-3p met both these requirements. hsa-miR-29c-3p was significantly downregulated in HNSC, and its upregulation was associated with better prognosis of patients with HNSC (f6, p < 0.05). Pearson’s correlation analysis showed that SFXN3 expression was also significantly and negatively correlated with hsa-miR-29c-3p expression (r = -0.32, p < 0.001) (f6). Taken together, hsa-miR-29c-3p might be an upstream miRNA of SFXN3 in HNSC.

Prediction of upstream lncRNAs of hsa-miR-29c-3p

The upstream lncRNAs of hsa-miR-29c-3p were further predicted using the ENCORI database. A total of 114 possible lncRNAs were predicted. Among them, six predicted lncRNAs were negatively associated with the expression of hsa-miR-29c-3p in HNSC, including AC098828.2 (ENSG00000223734), LINC01270 (ENSG00000203999), MIR193BHG (ENSG00000262454), PVT1 (ENSG00000249859), NOP14-AS1 (ENSG00000249673), and LINC01907 (ENSG00000226125). However, only four lncRNAs were positively associated with the expression of SFXN3 (f7), which suggested that the lncRNAs were associated with a competing endogenous RNA (ceRNA) mechanism. Then, the expression levels of these four related lncRNAs in HNSC were assessed using both ENCORI and the GEPIA database. As suggested in f7, LINC01270, NOP14-AS1, and MIR193BHG were significantly upregulated in HNSC compared with that in normal controls and were consistent in the two databases. Subsequently, the prognostic values of these four lncRNAs were evaluated in HNSC. As shown in f7, patients with HNSCs with a higher expression of LINC01270 were significantly related to poorer OS (p = 0.041).

Based on the results of the correlation, expression, and survival analyses, LINC01270 was considered to be the most likely upstream lncRNA of the hsa-miR-29c-3p/SFXN3 axis in HNSC.

Discussion

In this study, the expression analysis and survival analysis for SFXN3 indicated that HNSC patients with a high expression of SFXN3 had poor prognosis. Correlation analysis and comparative analysis between groups based on SFXN3 expression indicated that SFXN3 might function as a paclitaxel and docetaxel resistance predictor and is related to enriched tumor‐infiltrating macrophages, immune-suppressive cells, and immune checkpoint expression.

SFXN1 and SFXN3 are the main mitochondrial serine transporters in human cells. Both of them are likely to be regulated by the Myc transcription factor, which is a well-known oncogene. Serine transported to mitochondria is converted into glycine and formate, while generating reactive one-carbon units. Glycine can also provide one-carbon units to the folate cycle (ref. 10). Although SFXN3 is the closest homolog of SFXN1, loss of SFXN3 causes a more severe defect in glycine synthesis defect compared with loss of SFXN1. Besides, SFXN1 has no prognostic value for HNSC (Figure S4). Moreover, SFXN3, in addition to serine, can also transport alanine or cysteine. Cysteine is another product of serine metabolism. The two products of serine metabolism, glycine and cysteine, are precursors for glutathione synthesis. Moreover, serine and glycine are essential for tumor growth and survival (ref. 48). Metabolic requirements evolve during cancer progression, and new dependencies emerged in the context of therapy resistance and metastasis (ref. 49).

Both paclitaxel and docetaxel are microtubule-stabilizing agents, which are constituents of a classical chemotherapy regimen according to the guide for the treatment of HNSC (ref. 2). Recently, resistance to docetaxel in non-small cell lung cancer cells was reported to be closely related to serine and glycine metabolism and cysteine metabolism. More serine was found to be accumulated in the docetaxel-resistant non-small cell lung cancer cells compared with the levels in the sensitive cells (ref. 50). However, the abundance of enzymes involved in glutathione biosynthesis in docetaxel-resistant cells was dramatically lower than that in sensitive cells, which ultimately reduced the glutathione levels and elevated the levels of reactive oxygen species. Increased oxidative stress stimulated the expression and function of permeability glycoproteins and thus promoted drug resistance (ref. 50). Furthermore, serine and glycine metabolism was found to be increased in the more radioresistant HNSC cancer cells (ref. 51). These observations might indicate that HNSC cancer cells that are more resistant to therapy might have altered the metabolic pathways that control redox status, DNA repair, and DNA methylation after radiation (ref. 49). In our study, we found that SFXN3 expression correlated with paclitaxel resistance by the analysis with data from head and neck tumor cell lines, and there was also a trend toward correlation with docetaxel resistance. Further investigations are needed to clarify whether resistance to paclitaxel and docetaxel in HNSC is associated with serine and glycine metabolism.

The standard clinical practices for HNSC have changed substantially in the era of immunotherapy. HNSC is a malignant tumor in which the immune surveillance mechanism is suppressed (ref. 41). Our study suggested that SFXN3 expression correlated significantly and positively with biomarkers of the main immune-suppressive cells including MDSCs, TAMs, and Tregs, in the tumor microenvironment. In oral squamous cell carcinoma, SFXN3 is not only expressed in cancer cells but also expressed in components of the tumor microenvironment, such as the stromal fibroblasts and the endothelial cells of the small arteries in the cancer nest (ref. 15). Increased levels of serine and glycine in tumor tissue were demonstrated to promote the survival of tumor stromal cells and vascular epithelial cells, which can establish a protective niche for the maintenance of the tumor (ref. 48). This may be one of the reasons why SFXN3 is related to the tumor microenvironment. Second, the inflammatory response in macrophages is supported by serine metabolism. Inflammatory cytokine IL-1β, produced by macrophages, is associated with tumorigenesis in HNSC (ref. 52). Rodriguez et al. showed that glycine, which is made from serine, could support glutathione synthesis in macrophages, leading to the production of IL-1β. Furthermore, serine also supplies glycine and one-carbon units to support T-cell proliferation (ref. 53). However, there is a competitive relationship between tumor cells and T cells, and perhaps macrophages, for cysteine and glycine, which might be related to the suppression of T-cell activation and proliferation, leading to an immunosuppressive microenvironment (ref. 54).

Extensive research has shown that ncRNAs, including miRNAs, lncRNAs, and circular RNAs (circRNAs), can act as regulators of gene expression. Many ncRNAs function as competitive endogenous RNAs (ceRNAs) to regulate the biological behaviors of cancer (ref. 47, ref. 55). To clarify whether a ceRNA regulatory mechanism existed in the regulation of SFXN3, we employed prediction programs to predict possible miRNAs that could potentially bind to SFXN3 and possible related lncRNAs. The LINC01270/hsa-miR-29c-3p/SFXN3 axis was identified as a potential regulatory pathway in HNSC. LINC01270 was upregulated and was significantly associated with worse OS in HNSC, which was consistent with the role of LINC01270 in lung adenocarcinoma, endometriosis, esophageal cancer, and breast cancer (ref. 56–ref. 59). Correspondingly, hsa-miR-29c-3p was downregulated and was significantly associated with better survival in HNSC in the present study, which was consistent with previous research in HNSC, bladder cancer, and colorectal cancer (ref. 60–ref. 63).

Our study had several limitations. First, the heterogeneity of head and neck cancer varies. The subsites of head and neck tumors and treatment modalities might have affected the identification of significant prognostic variables. Second, the potential upstream axis of SFXN3 and the estimation of infiltrated immune cell types were computationally predicted. Experimental exploration and multicenter and prospective studies are needed to further verify the clinical value of SFXN3 and its underlying mechanisms.

Conclusion

We demonstrated that SFXN3 was highly expressed and correlated with unfavorable prognosis in HNSC. We identified the LINC01270/hsa-miR-29c-3p/SFXN3 axis as a potential regulatory pathway in HNSC (f8). Furthermore, SFXN3 might exert its oncogenic roles by increasing resistance to paclitaxel and promoting the immunosuppressive microenvironment in HNSC. However, further basic research and clinical trials are needed to validate these findings.

Data availability statement

The R codes are accessible on Jianguoyun at: https://www.jianguoyun.com/p/DTLfw3QQkdi8ChiJsbkEIAA.

Author contributions

KC and YL designed the article. KC, SG and QL organized the public data and wrote the manuscript. KC, XF, JL and YL took charge for data visualization. MY, SH and YZ obtained the clinical information. NL and JM contributed to the concept and revised the article. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by grants from the Natural Science Foundation of Guangdong Province (2017A030312003).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. H Sung, J Ferlay, RL Siegel, M Laversanne, I Soerjomataram, A Jemal. Global cancer statistics 2020: Globocan estimates of incidence and mortality worldwide for 36 cancers in 185 countries.. CA Cancer J Clin (, 2021. [DOI]
2. MD Mody, JW Rocco, SS Yom, RI Haddad, NF Saba. Head and neck cancer.. Lancet (, 2021. [DOI | PubMed]
3. V Budach, I Tinhofer. Novel prognostic clinical factors and biomarkers for outcome prediction in head and neck cancer: A systematic review.. Lancet Oncol (, 2019. [DOI]
4. AA Bhat, P Yousuf, NA Wani, A Rizwan, SS Chauhan, MA Siddiqi. Tumor microenvironment: An evil nexus promoting aggressive head and neck squamous cell carcinoma and avenue for targeted therapy.. Signal Transduct Target Ther (, 2021. [DOI | PubMed]
5. JW Carlisle, CE Steuer, TK Owonikoko, NF Saba. An update on the immune landscape in lung and head and neck cancers.. CA Cancer J Clin (, 2020. [DOI]
6. X Li, D Han, R Kin Ting Kam, X Guo, M Chen, Y Yang. Developmental expression of sideroflexin family genes in xenopus embryos.. Dev Dyn (, 2010. [DOI]
7. PJ Lockhart, B Holtom, S Lincoln, J Hussey, A Zimprich, T Gasser. The human sideroflexin 5 (Sfxn5) gene: Sequence, expression analysis and exclusion as a candidate for Park3.. Gene (, 2002. [DOI]
8. MG Acoba, ESS Alpergin, S Renuse, L Fernández-Del-Río, YW Lu, O Khalimonchuk. The mitochondrial carrier Sfxn1 is critical for complex iii integrity and cellular metabolism.. Cell Rep (, 2021. [DOI | PubMed]
9. A Rivell, RS Petralia, YX Wang, MP Mattson, PJ Yao. Sideroflexin 3 is a mitochondrial protein enriched in neurons.. Neuromolecular Med (, 2019. [DOI]
10. CF Labuschagne, NJ van den Broek, GM Mackay, KH Vousden, OD Maddocks. Serine, but not glycine, supports one-carbon metabolism and proliferation of cancer cells.. Cell Rep (, 2014. [DOI]
11. B Minjarez, KG Calderón-González, ML Rustarazo, ME Herrera-Aguirre, ML Labra-Barrios, DE Rincon-Limas. Identification of proteins that are differentially expressed in brains with alzheimer’s disease using itraq labeling and tandem mass spectrometry.. J Proteomics (, 2016. [DOI]
12. IS Amorim, LC Graham, RN Carter, NM Morton, F Hammachi, T Kunath. Sideroflexin 3 is an α-Synuclein-Dependent mitochondrial protein that regulates synaptic morphology.. J Cell Sci (, 2017. [DOI]
13. MR Sullivan, KR Mattaini, EA Dennstedt, AA Nguyen, S Sivanand, MF Reilly. Increased serine synthesis provides an advantage for tumors arising in tissues where serine levels are limiting.. Cell Metab (, 2019. [DOI | PubMed]
14. M Yang, KH Vousden. Serine and one-carbon metabolism in cancer.. Nat Rev Cancer (, 2016. [DOI]
15. R Murase, Y Abe, T Takeuchi, M Nabeta, Y Imai, Y Kamei. Serum autoantibody to sideroflexin 3 as a novel tumor marker for oral squamous cell carcinoma.. Proteomics Clin Appl (, 2008. [DOI]
16. Y Li, X Ge, F Peng, W Li, JJ Li. Exaggerated false positives by popular differential expression methods when analyzing human population samples.. Genome Biol (, 2022. [DOI | PubMed]
17. W Shen, Z Song, Z Xiao, M Huang, D Shen, P Gao. Sangerbox: A comprehensive, interaction-friendly clinical bioinformatics analysis platform.. iMeta (, 2022. [DOI]
18. DS Chandrashekar, B Bashel, SAH Balasubramanya, CJ Creighton, I Ponce-Rodriguez, B Chakravarthi. Ualcan: A portal for facilitating tumor subgroup gene expression and survival analyses.. Neoplasia (, 2017. [DOI]
19. Z Tang, C Li, B Kang, G Gao, C Li, Z Zhang. Gepia: A web server for cancer and normal gene expression profiling and interactive analyses.. Nucleic Acids Res (, 2017. [DOI | PubMed]
20. T Li, J Fan, B Wang, N Traugh, Q Chen, JS Liu. Timer: A web server for comprehensive analysis of tumor-infiltrating immune cells.. Cancer Res (, 2017. [DOI]
21. B Chen, MS Khodadoust, CL Liu, AM Newman, AA Alizadeh. Profiling tumor infiltrating immune cells with cibersort.. Methods Mol Biol (, 2018. [DOI]
22. AM Newman, CL Liu, MR Green, AJ Gentles, W Feng, Y Xu. Robust enumeration of cell subsets from tissue expression profiles.. Nat Methods (, 2015. [DOI]
23. D Aran, Z Hu, AJ Butte. Xcell: Digitally portraying the tissue cellular heterogeneity landscape.. Genome Biol (, 2017. [DOI | PubMed]
24. D Aran, M Sirota, AJ Butte. Systematic pan-cancer analysis of tumour purity.. Nat Commun (, 2015. [DOI | PubMed]
25. K Yoshihara, M Shahmoradgoli, E Martínez, R Vegesna, H Kim, W Torres-Garcia. Inferring tumour purity and stromal and immune cell admixture from expression data.. Nat Commun (, 2013. [DOI | PubMed]
26. SL Carter, K Cibulskis, E Helman, A McKenna, H Shen, T Zack. Absolute quantification of somatic DNA alterations in human cancer.. Nat Biotechnol (, 2012. [DOI]
27. WC Reinhold, M Sunshine, H Liu, S Varma, KW Kohn, J Morris. Cellminer: A web-based suite of genomic and pharmacologic tools to explore transcript and drug patterns in the nci-60 cell line set.. Cancer Res (, 2012. [DOI]
28. T Hastie, R Tibshirani, B Narasimhan, G Chu. Impute: Imputation for microarray data. r package version 1.58.0. 2019. (, 2019
29. ME Ritchie, B Phipson, D Wu, Y Hu, CW Law, W Shi. Limma powers differential expression analyses for rna-sequencing and microarray studies.. Nucleic Acids Res (, 2015. [DOI | PubMed]
30. W Yang, J Soares, P Greninger, EJ Edelman, H Lightfoot, S Forbes. Genomics of drug sensitivity in cancer (Gdsc): A resource for therapeutic biomarker discovery in cancer cells.. Nucleic Acids Res (, 2013. [DOI]
31. IS Vlachos, MD Paraskevopoulou, D Karagkouni, G Georgakilas, T Vergoulis, I Kanellos. Diana-Tarbase V7.0: Indexing more than half a million experimentally supported Mirna:Mrna interactions.. Nucleic Acids Res (, 2015. [DOI]
32. Y Chen, X Wang. Mirdb: An online database for prediction of functional microrna targets.. Nucleic Acids Res (, 2020. [DOI]
33. CE Vejnar, M Blum, EM Zdobnov. Mirmap web: Comprehensive microrna target prediction online.. Nucleic Acids Res (, 2013. [DOI]
34. H Dweep, N Gretz. Mirwalk2.0: A comprehensive atlas of microrna-target interactions.. Nat Methods (, 2015. [DOI | PubMed]
35. KC Miranda, T Huynh, Y Tay, YS Ang, WL Tam, AM Thomson. A pattern-based method for the identification of microrna binding sites and their corresponding heteroduplexes.. Cell (, 2006. [DOI]
36. B Fromm, T Billipp, LE Peck, M Johansen, JE Tarver, BL King. A uniform system for the annotation of vertebrate microrna genes and the evolution of the human micrornaome.. Annu Rev Genet (, 2015. [DOI]
37. MD Paraskevopoulou, G Georgakilas, N Kostoulas, IS Vlachos, T Vergoulis, M Reczko. Diana-Microt web server V5.0: Service integration into mirna functional analysis workflows.. Nucleic Acids Res (, 2013. [DOI]
38. HY Huang, YC Lin, J Li, KY Huang, S Shrestha, HC Hong. Mirtarbase 2020: Updates to the experimentally validated microrna-target interaction database.. Nucleic Acids Res (, 2020. [DOI]
39. R Li, H Qu, S Wang, JM Chater, X Wang, Y Cui. Cancermirnome: An interactive analysis and visualization database for mirnome profiles of human cancer.. Nucleic Acids Res (, 2021. [DOI]
40. JH Li, S Liu, H Zhou, LH Qu, JH Yang. Starbase V2.0: Decoding mirna-cerna, mirna-ncrna and protein-rna interaction networks from Large-scale clip-seq data.. Nucleic Acids Res (, 2014. [DOI]
41. RL Ferris. Immunology and immunotherapy of head and neck cancer.. J Clin Oncol (, 2015. [DOI]
42. R Mandal, Y Şenbabaoğlu, A Desrichard, JJ Havel, MG Dalin, N Riaz. The head and neck cancer immune landscape and its immunotherapeutic implications.. JCI Insight (, 2016. [DOI | PubMed]
43. AJ Petty, DH Owen, Y Yang, X Huang. Targeting tumor-associated macrophages in cancer immunotherapy.. Cancers (Basel) (, 2021. [DOI | PubMed]
44. BZ Qian, JW Pollard. Macrophage diversity enhances tumor progression and metastasis.. Cell (, 2010. [DOI | PubMed]
45. SK Biswas, A Mantovani. Macrophage plasticity and interaction with lymphocyte subsets: Cancer as a paradigm.. Nat Immunol (, 2010. [DOI]
46. G Morad, BA Helmink, P Sharma, JA Wargo. Hallmarks of response, resistance, and toxicity to immune checkpoint blockade.. Cell (, 2021. [DOI]
47. GJ Goodall, VO Wickramasinghe. Rna in cancer.. Nat Rev Cancer (, 2021. [DOI | PubMed]
48. SL Geeraerts, E Heylen, K De Keersmaecker, KR Kampen. The ins and outs of serine and glycine metabolism in cancer.. Nat Metab (, 2021. [DOI]
49. B Faubert, A Solmonson, RJ DeBerardinis. Metabolic reprogramming and cancer progression.. Science (, 2020. [DOI | PubMed]
50. SJ Li, B Cao, ZY Lu, RB Sun, SH Guo, Y Xie. Cystine supplementation rebalances the redox homeostasis of microenvironment in non-small cell lung cancer cells and reverses their resistance to docetaxel.. Acta Pharmacol Sin (, 2021. [DOI]
51. E Lindell Jonsson, I Erngren, M Engskog, J Haglöf, T Arvidsson, M Hedeland. Exploring radiation response in two head and neck squamous carcinoma cell lines through metabolic profiling.. Front Oncol (, 2019. [DOI | PubMed]
52. L Chen, CF Huang, YC Li, WW Deng, L Mao, L Wu. Blockage of the Nlrp3 inflammasome by Mcc950 improves anti-tumor immune responses in head and neck squamous cell carcinoma.. Cell Mol Life Sci (, 2018. [DOI]
53. EH Ma, G Bantug, T Griss, S Condotta, RM Johnson, B Samborska. Serine is an essential metabolite for effector T cell expansion.. Cell Metab (, 2017. [DOI | PubMed]
54. T Wang, G Liu, R Wang. The intercellular metabolic interplay between tumor and immune cells.. Front Immunol (, 2014. [DOI | PubMed]
55. N Yang, K Liu, M Yang, X Gao. Cernas in cancer: Mechanism and functions in a comprehensive regulatory network.. J Oncol (, 2021. [DOI | PubMed]
56. Y Wang, J Fu, Z Wang, Z Lv, Z Fan, T Lei. Screening key lncrnas for human lung adenocarcinoma based on machine learning and weighted gene Co-expression network analysis.. Cancer biomark (, 2019. [DOI]
57. L Liu, J Cheng, F Wei, L Pang, Z Zhi, W Yang. Lnc-rna Linc01279 induces endometriosis Via targeting of Hoxa10.. J Obstet Gynaecol Res (, 2021. [DOI]
58. N Li, Z Zhao, F Miao, S Cai, P Liu, Y Yu. Silencing of long non-coding rna Linc01270 inhibits esophageal cancer progression and enhances chemosensitivity to 5-fluorouracil by mediating Gstp1methylation.. Cancer Gene Ther (, 2021. [DOI]
59. J Ping, S Huang, J Wu, P Bao, T Su, K Gu. Association between lincrna expression and overall survival for patients with triple-negative breast cancer.. Breast Cancer Res Treat (, 2021. [DOI]
60. R Fang, Y Huang, J Xie, J Zhang, X Ji. Downregulation of mir-29c-3p is associated with a poor prognosis in patients with laryngeal squamous cell carcinoma.. Diagn Pathol (, 2019. [DOI | PubMed]
61. T Inamoto, H Uehara, Y Akao, N Ibuki, K Komura, K Takahara. A panel of microrna signature as a tool for predicting survival of patients with urothelial carcinoma of the bladder.. Dis Markers (, 2018. [DOI | PubMed]
62. S Zhang, J Jin, X Tian, L Wu. Hsa-Mir-29c-3p regulates biological function of colorectal cancer by targeting sparc.. Oncotarget (, 2017. [DOI]
63. K Hudcova, M Raudenska, J Gumulec, H Binkova, Z Horakova, R Kostrica. Expression profiles of mir-29c, mir-200b and mir-375 in tumour and tumour-adjacent tissues of head and neck cancers.. Tumour Biol (, 2016. [DOI]