Identification of a metabolic signature genes predicting overall survival for Head and neck squamous cell carcinoma

Abstract

Background: Head and Neck Squamous Cell Carcinoma (HNSCC) is a common head and neck malignancy that originates from lips, mouth, paranasal sinuses, oropharynx, larynx, nasopharynx, and other pharyngeal cancers. Since high-throughput expression data is now available, it is feasible to use global gene expression data to analyze the relationship between metabolic-related gene expression and clinical outcomes in HNSCC patients.

Method: In this study, we used RNA sequencing (RNA-seq) data from cancer genomic maps (TCGA) and validated in the GEO dataset to shine the dysfunctional metabolic microenvironment and given potential biomarkers for metabolic therapy.

Results: TCGA database contained 529 patients and 327 metabolic genes (198 upregulated and the largest logFC is CA9; 129 downregulated and the largest logFC is CA6. Cox regression model found 51 prognosis‐related genes and we found the most high-risk gene was LCLAT1 (HR=1.144, 95% CI=1.044~1.251); the most low-risk gene was CHDH (HR=0.580, 95% CI=0.400~0.839). We next study whether metabolic-related genes patterns could serve as an early predictor of incidence of HNSCC by ROC curve and the model demonstrated an AUC of 0.745 in TCGA and 0.618 in GEO. Meanwhile, the predictor of clinicopathological in HNSCC was also analyzed and we found the AUC for age, gender, grade, stage, T, M and N were 0.520, 0.495, 0.568, 0.606, 0.577, 0.476 and 0.673, respectively in TCGA and the AUC for age, gender, stage, T, M, N, smoking and HPV16-pos were 0.599, 0.531, 0.622, 0.606, 0.616, 0.550, 0.614, 0.519 and 0.397 respectively in GEO.

Conclusion: Our research found a novel metabolic signature gene for HNSCC prognosis prediction based on the TCGA dataset. Our signatures may reflect metabolic microenvironment disorders and provide useful biomarkers for metabolic therapy and response treatment prediction.

Background

Head and Neck Squamous Cell Carcinoma (HNSCC) is a head and neck malignancy that originates from lips, mouth, paranasal sinuses, larynx, nasopharynx, and other pharyngeal cancers¹. As the sixth most common type of malignant tumors, with more than 655,000 new cases and 90,000 deaths every year². Currently, smoking, drinking, and human papillomavirus (HPV) infection are considered risk factors for HNSCC occurrence and prognosis³. Unfortunately, the 5-year survival rate is still below 50% due to the usual lack of early symptoms when HNSCC is detected early, while the survival rate is reduced to 35% due to local recurrence and metastasis⁴. The occurrence and development of HNSCC is a complex process including multiple molecules. LASP1, a chaperone of heat shock protein family A member 1A (HSPA1A) can lead to colony formation, cell proliferation, invasion and cell cycle G2/M phase transition in HNSCC⁵. Kim et al reported that mechanistic and functional roles of CXCR7 as a key regulator of oncogenic TGF-β1/Smad2/3 signaling in HNSCC. In addition, Hsu et al⁶ notes the driver-oncogene role of atypical cadherin 1 (FAT1) in mediation of proliferation, cell-death evasion, chemoresistance and oncogenicity in oral squamous cell carcinoma (OSCC). As for the histological type of HNSCC and multiple anatomical sites, tumor markers vary widely. Therefore, it is possible to identify more valuable HNSCC drug targets by screening for changes in gene function networks linked to tumor formation and progression.

Loss of tumors suppressors and activation of oncogenes contribute to metabolic reprogramming in cancer, leading to improve nutrient uptake to provide energetic and biosynthetic pathways⁸. In the 1920s, Otto Warburg first reported that tumors took up distinctly more levels of glucose when compared with normal tissues indicating that these cells were shuttling glucose via the glycolytic fermentation pathway⁹. Recent studies have shown that immune cells have unique metabolic characteristics that affect their immune function. For example, macrophage polarization is associated with unique metabolic characteristics related to iron metabolism, energy metabolism, and lipid metabolism^10,11. More and more studies have investigated the prognostic role of these metabolic genes in cancer, but the role and mechanism of metabolism are still unclear. Hu et al¹² found that mutationally activated KRAS strikingly increased the glutathione biosynthesis and intracellular cystine level in lung adenocarcinoma. Yoo et al found that the SLC1A5 variant is a mitochondrial glutamine transporter used for metabolic reprogramming of pancreatic cancer and the knockout and overexpression of the SLC1A5 variant alters the growth of cancer cells and tumors, thereby supporting carcinogenesis¹³. However, there is no systematic assessment of metabolic-related genes and predicting the overall survival (OS) of HNSCC patients or the characteristics of response to immunotherapy. Since high-throughput expression data is now accessible, it is realistic to use global gene expression data to analyze the relationship between metabolic-related gene expression and clinical outcomes in HNSCC patients. In this study, we used RNA sequencing (RNA-seq) data from cancer genomic maps (TCGA) and validated in the Gene Expression Omnibus (GEO) dataset to show the dysfunctional metabolic microenvironment and given useful biomarkers for metabolic therapy.

Methods

Results

Built and verification of prognostic metabolic gene signatures

TCGA database contained 529 patients and 327 metabolic genes (198 upregulated and the largest logFC is CA9; 129 downregulated and the largest logFC is CA6; Table S1) were involved in TCGA‐LIHC to train the prognostic model. Meanwhile, the GSE65858 dataset contained 270 HNSCC tissue samples. The univariate Cox regression shown 51 survival‐related genes and we found the most high-risk gene was LCLAT1 (HR=1.144, 95% CI=1.044~1.251); the most low-risk gene was CHDH (HR=0.580, 95% CI=0.400~0.839) Table 1. Then Lasso‐penalized Cox analysis found 30 genes to build the prognostic model. The cohorts were then divided into high- and low-risk groups using the median risk score value as a cut-off.

Survival results and multivariate examination

OS analysis demonstrated that HNSCC with high risk group had a more terrible prognosis than that with low risk group (P<0.01) Fig. 1. Detailed prognostic signature information of HNSCC groups is visualized in Fig. 2, which indicated that the mortality rate was higher for HNSCC patients in the high-risk group and related to the lower OS. The result of univariate and multivariate Cox regression analysis showed that our prognostic model is an independent prognostic factor for OS Fig. 3. We next study whether metabolic-related genes patterns could serve as an early predictor of incidence of HNSCC by ROC curve and the model demonstrated an AUC of 0.745 in TCGA and 0.618 in GEO. Taken together, these results indicate that the prognostic model has moderate sensitivity and specificity. Meanwhile, the predictor of clinicopathological in HNSCC was also analyzed and we found the AUC for age, gender, grade, stage, T, M and N were 0.520, 0.495, 0.568, 0.606, 0.577, 0.476 and 0.673, respectively in TCGA and the AUC for age, gender, stage, T, M, N, smoking and HPV16-pos were 0.599, 0.531, 0.622, 0.606, 0.616, 0.550, 0.614, 0.519 and 0.397 respectively in GEO Fig 4.

Constructing and validating a predictive nomogram in TCGA and in GEO cohort

Nomogram was construct by involving clinicopathological and the prognostic model. Through the LASSO logistic regression algorithm, the most important prediction markers are selected through the training data set, and a strong contribution is made to the final prediction model. The model ultimately included 7 features in TCGA: age, gender, grade, stage, T, M and N and 8 features in TCGA: age, gender, stage, T, M, N, smoking and HPV16-pos Fig 5. Integrating our prognostic model with clinicopathological raise the forecasting sensitivity and specificity for forecasting 1‐, 3‐, and 5‐year OS and bring some net benefit which might useful for clinical management.

Gene set enrichment analyses

Then samples were divided into high‐ and low‐risk groups as training set to distinguish the potential function and elucidate the significant survival difference utilizing GSEA. Annotated gene sets c2.cp.kegg.v6.0.symbols.gmt was selected as the reference gene sets, which includes terms with NOM<0.05. Gene set permutations were executed multiple times for every examination. Great majority the metabolism‐related pathways were enriched in the high‐risk group (such as the metabolism of galactose, nicotinate/nicotinamide and pantothenate/COA biosynthesis or metabolic disease related such as the maturity‐onset diabetes of the young), whereas most of the non-metabolism related pathways were enriched in the low‐risk group (base excision repair, spliceosome, homologous recombination, nucleotide excision, DNA replication) Fig 6 and Table 2.

Online database analysis

We used multidimensional survey ways to explore CA6, CA9, LCLAT1 and CHDH based on thousands of variations in copy numbers or gene expressions in patients whom with HNSCC in detail to provide new insights into the potential functions, expression patterns, molecular mechanisms and distinct prognostic. Just like our results, CA6 were significantly under-expressed, while CA9 were significantly overexpressed in HNSCC cancer both in the TIMER database (Figure 7) Oncomine (Figure 8). In spite of lack in the Oncomine datas, the mRNA expression of LCLAT1 was also obvious overexpressed and CHDH under-expression in HNSCC in the TIMER database. Representative protein expressions of CA9, LCLAT1, and CHDH were explored in the HPA database, as shown in Figure 9. However, CA6 cannot be found on the website. CA9 has the most frequent genetic changes (10%), and amplified mutations are the most common changes (Figure 10). In summary, the abnormal expression of these genes has been further verified in HNSCC, and genetic changes may help explain the abnormal expression of these genes to a certain extent.

Discussion

Proliferating cancer cells must maintain sufficient energy and a library of metabolic intermediates to build the macromolecules required for proliferation, including DNA, proteins, and lipids¹⁶. Because the metabolic profile of tumor cells distinguishes them from normal cells and is essential for their growth and survival, metabolic signaling pathways have become ideal targets for therapeutic intervention in cancer patients. In this study, we identified a novel and effective metabolic-related prognostic signature gene based on the TCGA dataset, and verified its validity in the GSE65858 dataset, indicating that the signature has strong prognostic value. Our signature may represent the metabolic status of patients with HNSCC and provide potential biomarkers and targets for metabolic signaling pathways therapeutic intervention.

In our study, we downloaded transcriptome data from the TCGA database and verification in GEO dataset, metabolic-related genes extract from the GSEA metabolic signaling pathways. We first study the connections between differentially expressed of RNA sequencing, immune-related genes and transcription factors in LUAD patients. Univariate Cox regression model found 51 survival-related genes and we found the most high-risk gene was LCLAT1, the most low-risk gene was CHDH. Lysocardiolipin acyltransferase 1 (LCLAT1, also known as LYCAT; AGPAT8; ALCAT1) whichis an Acyl-CoA:lyso-CL acyltransferase and an isoform of MLCLAT1¹⁷. Cardiolipin (CL) types of polyunsaturated fatty acids, especially DHA (C22: 6n3), increased in ALCAT1-expressing cells, while C16-C18 fatty acids decreased significantly¹⁸. A recent study found that ALCAT1 is critical for coupling mitochondrial respiration and metabolic plasticity. Wang et al reported that forced expression of ALCAT1 in primary hepatocytes led to multiple defects including steatosis, defective autophagy, and mitochondrial dysfunction²⁰. Until now, there was no study about the role of ALCAT1 in tumors, thus we suppose ALCAT1 plays a regulatory role in cardiolipin remodeling in response to oxidative stress and stimulate mitochondrial activity in HNSCC cancers and needs more relate studies to conduct. Choline dehydrogenase (CHDH) encoded a choline dehydrogenase that localizes to the mitochondrion and variations in this gene can affect susceptibility to choline deficiency. CHDH strongly predicted clinical outcome in breast cancer patients receiving tamoxifen monotherapy²¹. Choline is an essential nutrient required for methyl group metabolism and CHDH associated with an increased risk of breast cancer²². While there was no study focus on the role of CHDH in HNSCC and more researches are needed.

We next study whether metabolic-related genes patterns could serve as an early predictor of incidence of HNSCC by ROC curve and the model demonstrated an AUC of 0.745 in TCGA and 0.618 in GEO. Integrating our prognostic model with clinicopathological raise the predicting sensitivity and specificity for predicting 1-, 3‐, and 5‐year OS and bring some net benefit which might help clinical management. Then we try to distinguish the potential function and elucidate the significant survival difference utilizing GSEA. Great majority of the metabolism‐related pathways were enriched in the high‐risk group (such as the metabolism of galactose, nicotinate/nicotinamide and pantothenate/COA biosynthesis or metabolic disease related such as the maturity‐onset diabetes of the young), whereas most of the non-metabolism related pathways were enriched in the low‐risk group (base excision repair, spliceosome, homologous recombination, nucleotide excision, DNA replication). Galactose is essential for human metabolism and plays a definitive role in energy transfer and galactosylation of complex molecules. Nicotinamide adenine dinucleotide (NAD) plays a central role in energy metabolism and integrates cell metabolism with signaling and gene expression²³. NAD biosynthesis is dependent on nicotinamide/nicotinate single nucleotide adenylate transferase²⁴. From this perspective, high-risk patients may benefit from metabolic therapy, while low-risk patients may benefit from non-metabolic therapy. However, much effort is needed to study the relationship between signatures and metabolic microenvironments and metabolic therapies. On the other hand, the results provide a promising direction for elucidating the underlying molecular mechanisms of signatures. In conclusion, our signatures may reflect metabolic microenvironment disorders and provide potential biomarkers for metabolic therapy and treatment response prediction.

TCGA database contained 529 patients and 327 metabolic genes (198 upregulated and the largest logFC is CA9; 129 downregulated and the largest logFC is CA6. Carbonic anhydrases (CAs) are a large class of zinc metal enzymes that catalyze the reversible hydration of carbon dioxide. They are involved in various biological processes, including bone resorption, respiration, calcification, and acid-base balance. Study found that pancreatic ductal adenocarcinomas (PDACs) cells that express activated KRAS increase expression of CA9, via stabilization of hypoxia-inducible factor 1 subunit alpha (HIF1A) and HIF2A, to regulate pH and glycolysis²⁵. Similarly, CA9 is an independent prognostic factor for OSCC patients and therefore a potential therapeutic target²⁶. Carbonic anhydrase 6 (CA6, also known as CA-VI; GUSTIN) encodes several isoenzymes of carbonic anhydrase, which are only found in salivary glands. Saliva and proteins may play a role in the reversible hydration of carbon dioxide. CA6 is a specific marker of salivary gland serous acinar cells and acinar cell carcinoma (AciCC). CA6 has the same sensitivity and specificity as DOG1 in the differential diagnosis of AciCC and breast analogs (MASC)²⁷. However, the role of CA9 and CA6 in human HNSCC prognostic remains to be elucidated. Thus, carbonic anhydrases may help us understand the protein disease pathogenesis and progression in future study. Our study identified a novel metabolic signature gene for HNSCC prognosis prediction base on TCGA data set. Our signatures may reflect metabolic microenvironment disorders and provide potential biomarkers for metabolic therapy and treatment response prediction. However, signatures must be verified in more independent cohorts and functional experiments of predictive metabolic genes. At the same time, our study has limitations. On the one hand, our results have not been verified in clinical samples, and on the other hand, because of the relatively small number of patients, our results cannot provide accurate clinical data. We will carry out a lot of related research in the follow-up work, and finally translate it into the development of new strategies for precision cancer medicine.

Conclusion

Our research found a novel metabolic signature gene for HNSCC prognosis prediction based on the TCGA dataset. Our signatures may reflect metabolic microenvironment disorders and provide useful biomarkers for metabolic therapy and response treatment prediction.

Abbreviations

HNSCC:head and neck squamous cell carcinoma; TCGA:cancer genome atlas; HPV:human papillomavirus; ROC:receiver operating characteristic; GSEA:Gene set enrichment analysis; GEO:Gene Expression Omnibus; OS:over survival.

Declarations

Ethics approval and consent to participate

Not applicable.

Consent to public

Not applicable.

Availability of data and materials

Data sharing is not applicable to this article as no datasets were generated or analysed during the current study.

Competing interests

The authors declare that they have no competing interests.

Funding

This work is not supported by grants.

Authors’ contributions

W.Z.H. and Z.W. designed and analyzed the research study; W.Z.H. wrote and revised the manuscript, W.Z.H. collected the data and all authors have read and approved the manuscript.

Acknowledgements

Not applicable

References

1. Xiao-Nan F, Miao Y, Hua L, et al. Comprehensive analysis of competitive endogenous RNAs network associated with head and neck squamous cell carcinoma. Scientific Reports, 2018, 8(1):10544-. https://doi.org/10.1038/s41598-018-28957-y.
2. Ferlay, J. et al. Cancer incidence and mortality worldwide: sources, methods and major patterns in GLOBOCAN 2012. Int. J. Cancer 136, E359–E386.
3. Castellsague, X. et al. HPV involvement in head and neck cancers: comprehensive assessment of biomarkers in 3680 patients. J. Natl. Cancer Inst. 108, djv403.
4. Posner, M. & Vermorken, J. B. Induction therapy in the modern era of combined-modality therapy for locally advanced head and neck cancer. Seminars in oncology 35, 221–228, https://doi.org/10.1053/j.seminoncol.2008.03.007 (2008).
5. Chen Q, Wu K, Qin X et al. LASP1 promotes proliferation, metastasis, invasion in head and neck squamous cell carcinoma and through direct interaction with HSPA1A. J Cell Mol Med. 2019 Dec 2. doi: 10.1111/jcmm.14854.
6. Kim N, Ryu H, Kim S et al. CXCR7 promotes migration and invasion in head and neck squamous cell carcinoma by upregulating TGF-β1/Smad2/3 signaling. Sci Rep. 2019 Dec 2;9(1):18100. doi: 10.1038/s41598-019-54705-x.
7. Hsu TN, Huang CM, Huang CS et al. Targeting FAT1 Inhibits Carcinogenesis, Induces Oxidative Stress and Enhances Cisplatin Sensitivity through Deregulation of LRP5/WNT2/GSS Signaling Axis in Oral Squamous Cell Carcinoma. Cancers (Basel). 2019 Nov 27;11(12). pii: E1883. doi: 10.3390/cancers11121883.
8. Boroughs, L., DeBerardinis, R. Metabolic pathways promoting cancer cell survival and growth. Nat Cell Biol 17, 351–359 (2015). doi:10.1038/ncb3124.
9. Hanahan, D., & Weinberg, R. A. (2011). Hallmarks of cancer: The next generation. Cell, 144(5), 646–674.
10. Jha, A.K., Huang, S.C., Sergushichev, A., Lampropoulou, V., Ivanova, Y., Loginicheva, E., Chmielewski, K., Stewart, K.M., Ashall, J., Everts, B., et al. (2015). Network Integration of Parallel Metabolic and Transcriptional Data Reveals Metabolic Modules that Regulate Macrophage Polarization. Immunity 42, 419–430.
11. Biswas, S.K., and Mantovani, A. (2010). Macrophage plasticity and interaction with lymphocyte subsets: cancer as a paradigm. Nat. Immunol. 11, 889–896.
12. Hu K, Li K, Lv J et al. Suppression of the SLC7A11/glutathione axis causes synthetic lethality in KRAS-mutant lung adenocarcinoma. J Clin Invest. 2019 Dec 24. pii: 124049. doi: 10.1172/JCI124049.
13. Yoo HC, Park SJ, Nam M et al. A Variant of SLC1A5 Is a Mitochondrial Glutamine Transporter for Metabolic Reprogramming in Cancer Cells. Cell Metab. 2019 Dec 13. pii: S1550-4131(19)30664-3. doi: 10.1016/j.cmet.2019.11.020.
14. Tibshirani, R. (1997). The lasso method for variable selection in the Cox model. Statistics in Medicine, 16(4), 385–395。
15. Iasonos A, Schrag D, Raj G V, et al. How To Build and Interpret a Nomogram for Cancer Prognosis[J]. Journal of Clinical Oncology, 2008, 26(8):1364-1370.
16. Weyandt J D, Thompson C B, Giaccia A J, et al. Metabolic Alterations in Cancer and Their Potential as Therapeutic Targets[J]. American Society of Clinical Oncology - Educational Book, 2017, 37:825-832.
17. Cao, Y. Liu, J. Lockwood, P. Burn, and Y. Shi, “A novel cardiolipin-remodeling pathway revealed by a gene encoding an endoplasmic reticulum-associated acyl-CoA:lysocardiolipin acyltransferase (ALCAT1) in mouse,” Journal of Biological Chemistry, vol. 279, no. 30, pp. 31727–31734, 2004. View at Publisher · View at Google Scholar · View at Scopus
18. Li, C. Romestaing, X. Han et al., “Cardiolipin remodeling by ALCAT1 links oxidative stress and mitochondrial dysfunction to obesity,” Cell Metabolism, vol. 12, no. 2, pp. 154–165, 2010. View at Publisher · View at Google Scholar · View at Scopus
19. Rieger B, Krajčová A, Duwe P et al. ALCAT1 Overexpression Affects Supercomplex Formation and Increases ROS in Respiring Mitochondria. Oxid Med Cell Longev. 2019 Dec 6;2019:9186469. doi: 10.1155/2019/9186469. eCollection 2019.
20. Wang L, Liu X, Nie J et al. ALCAT1 controls mitochondrial etiology of fatty liver diseases, linking defective mitophagy to steatosis[J]. Hepatology, 2015, 61(2).
21. Xiao-Jun M, Hilsenbeck S G, Wilson W, et al. The HOXB13:IL17BR expression index is a prognostic factor in early-stage breast cancer[J]. Journal of Clinical Oncology Official Journal of the American Society of Clinical Oncology, 2006, 24(28):4611-9.
22. Xu X, Gammon M D, Zeisel S H, et al. Choline metabolism and risk of breast cancer in a population-based study[J]. The FASEB Journal, 2008, 22(6):2045-2052.
23. Coelho A I, Berry G T, Rubio-Gozalbo M E. Galactose metabolism and health[J]. Current Opinion in Clinical Nutrition & Metabolic Care, 2015, 18(4):422-7.
24. Ortiz-Joya L, Contreras-Rodríguez LE, Ramírez-Hernández MH. Et al. Protein-protein interactions of the nicotinamide/nicotinate mononucleotide adenylyltransferase of Leishmania braziliensis. Mem Inst Oswaldo Cruz. 2019 Mar 21;114:e180506. doi: 10.1590/0074-02760180506.
25. McDonald PC, Chafe SC, Brown WS et al. Regulation of pH by Carbonic Anhydrase 9 Mediates Survival of Pancreatic Cancer Cells With Activated KRAS in Response to Hypoxia. Gastroenterology. 2019 Sep;157(3):823-837.
26. Eckert AW, Horter S, Bethmann D et al. Investigation of the Prognostic Role of Carbonic Anhydrase 9 (CAIX) of the Cellular mRNA/Protein Level or Soluble CAIX Protein in Patients with Oral Squamous Cell Carcinoma. Int J Mol Sci. 2019 Jan 16;20(2).
27. Hsieh M S, Jeng Y M, Jhuang Y L, et al. Carbonic anhydrase VI: a novel marker for salivary serous acinar differentiation and its application to discriminate acinic cell carcinoma from mammary analogue secretory carcinoma of the salivary gland[J]. Histopathology, 2016, 68(5):641-647.

Tables

Table 1. Prognostic related metabolic genes, HR> 1 is a high-risk gene, and HR <1 is a low-risk gene.

Table 2. Gene sets enriched in phenotype high and low.