A Novel Lipid Prognostic Signature of ADCY2, LIPE, and OLR1 in Head and Neck Squamous Cell Carcinoma

Simple Summary Clinically, aberrant lipid metabolism is responsible for overweight and/or obesity. Overweight is considered as an independent factor of cancer risk in 2019. Therefore, lipid metabolic reprogramming is an emerging hallmark of malignancy. It is an urgent need to comprehensively understand the relationship among lipid metabolism and HNSCC and identify a valuable biomarker for predicting prognosis of HNSCC patients. Three new findings were found in this study. Firstly, we identified the lipid-related differentially expressed genes (DEGs) by using the GEO microarrays and TCGA dataset. A novel lipid-related mRNA prognostic signature (LRPS, consisting of ADCY2, LIPE and OLR1) was developed, which could predict the survival and prognosis of HNSCC patients as an independent effective prognostic factor. Secondly, we found that the LRPS could indicate the type of infiltrated immune cells in HNSCC tumor microenvironment. Thirdly, we verified that the LPPS score could interpret the TP53 status of HNSCC. Our new findings indicated that LRPS has a potential to be a promising indicator of overall survival, TP53 status, and immune characteristics in HNSCC, and perhaps can monitor and guide the treatment efficacy and prognosis of HNSCC in the future. Background Head and neck squamous cell carcinoma (HNSCC) is characterized by a high frequency of lymph node metastasis and a high mortality. Lipid metabolic reprogramming is an emerging carcinogen as its role in fulfilling cancer growth and spread. However, little is known about the correlation between lipid metabolism and HNSCC. Materials and Methods Expressions of lipid-related genes were obtained from the Cancer Genome Atlas (TCGA) and Gene expression Omnibus (GEO) databases for differential and functional analyses. A total number of 498 patients from TCGA with complete information were included to identify a lipid-related prognostic signature (LRPS), based on ADCY2, LIPE, and OLR1, by using univariate and multivariate Cox regression analyses. LRPS-high and LRPS-low groups were accordingly divided to pathway and cell enrichment analyses. Results LRS-low patients had a better overall survival and relapse - free survival than LRS-high ones in HNSCC. The LRPS-high group was significantly related to perineural invasion of cancer, cancer-related pathways, high TP53 mutation rate, high proportion of natural killer T cells (NKT), dendritic cells, monocytes, Treg, and M1 and M2 macrophage infiltration in HNSCC tumor tissues. Conversely, the LRPS-low group correlated with DNA damage-related and T-cell-regulated pathways, low frequency of mutated TP53, and high infiltration of B cells and CD4+ effector cells including Th1 and Th2. Conclusion LRPS has a potential to be a promising indicator of overall survival, prognosis, TP53 status, and immune characteristics in HNSCC.


INTRODUCTION
Head and neck squamous cell carcinoma (HNSCC) is the most common type of the head and neck cancers, with a high risk for recurrence and poor survival under the advanced treatment approaches. The incidence of HNSCC was increased by 36.3% during the past 10 years, from~482,000 HNSCC patients in 2008 to~657,000 cases in 2018 (1,2). Smoking, alcohol assumption, and virus infection are recognized as important carcinogenic factors (3). Recent studies implicate that abnormal lipid metabolism may be related with HNSCC development and progression (4,5).
Lipid metabolic reprogramming is an emerging hallmark of malignancy (6). Overwhelming lipid anabolic and catabolic processes are essential for the uncontrolled cell proliferation and rapid cancer growth. Simultaneously, lipids constitute most of the cell membranes and serve as signaling molecules. Theoretically, fatty acids and cholesterol synthesis provide carcinogenesis and metastasis with a range of metabolic fuels and substrates, as well as pro-tumor signaling cytokines (7)(8)(9)(10)(11). Furthermore, the roles of lipid metabolites in protecting cancer cells from harmful conditions (like endoplasmic reticulum stress, reactive oxygen species, and drug toxicity) have been substantiated in various cancers (12,13). Some oncolipidactivated signaling pathways, such as sterol regulatory element-binding proteins and stearoyl-CoA desaturases, have been identified to be the potential targets for cancer treatment in the future (6,14,15).
Clinically, aberrant lipid metabolism is responsible for overweight and obesity. Overweight is considered as an independent factor of cancer risk by the American Cancer Society, which released a report entitled Cancer Facts & Figures in 2019 (16). Nowadays, it is estimated that 5% of cancers in men and 11% in women are attributed to overweight (17). Experimental evidence indicates that high-fat diet-induced obesity not only promotes carcinogenesis, but also induces lymphangiogenesis and lymphatic metastasis in vivo (18)(19)(20). Conversely, diet-caused weight loss was shown to reduce cancer risk (21). Furthermore, a deliberate weight loss has been proved to reverse the effects of obesity-induced oxidative stress, inflammatory activities, and oncogenesis (22). Reduction of DNA damage responses in overweight mice was also observed after an administration of energy restriction (23).
In this study, we firstly identified the lipid-related differentially expressed genes (DEGs) by using the GEO microarrays and the TCGA dataset. A novel lipid-related mRNA prognostic signature (LRPS, consisting of ADCY2, LIPE, and OLR1) was developed for predicting survival of HNSCC patients. Accordingly, HNSCC patients were divided into high-risk and low-risk groups according to their LRPS signature, and gene-set enrichment analysis (GSEA) and cell enrichment analysis were used to elucidate the potential mechanisms.

Ethics Approval
The original datasets in our study were downloaded from the TCGA database and GEO dataset. We downloaded and analyzed the study data in accordance with the relevant data policies of TCGA database and GEO datasets, and therefore, no additional ethics approval was needed.

Data Source
The original datasets comparing the mRNA expression profiles between tumors and adjacent normal tissues were obtained from the three GEO databases [GSE30784 (containing 167 oral squamous cell carcinoma, 17 dysplasia, and 45 normal oral tissues), GSE37991 (containing 40 male oral squamous cell carcinoma biopsies), and GSE65858 (containing 290 HNSCC biopsies)] and a TCGA dataset (containing 498 HNSCC biopsies). The clinical samples from the TCGA database with complete clinical information of patients were selected. The microarray data of GSE30784, GSE37991, and GSE65858 were based on GPL570 (Affymetrix Human Genome U133 Plus 2.0 Array), GPL6883 (Illumina HumanRef-8 v3.0 expression beadchip), and GPL10558 (Illumina HumanHT-12 V4.0 expression beadchip), respectively. The corresponding clinical information of patients with HNSCC was also acquired from the TCGA database (up to July 19, 2019). A total of 498 HNSCC patients with detailed follow-up time were included for the following analyses.

Functional Enrichment Analyses
Gene ontology (GO) and KEGG pathway enrichment analyses were performed using the DAVID online database (the Database for Annotation, Visualization, and Integration Discovery) (24,25). The enriched biological processes (BP), cellular component (CC), and molecular function (MF) were obtained to analyze the DEGs.

Protein-Protein Interaction Analysis
The STRING online database (http://string-db.org) was performed for PPI analysis. Cytoscape software was employed to construct the PPI network (26). MCODE tool of Cytoscape was performed to identify gene cluster of the PPI network. Degree cutoff ≥ 2, node score cutoff ≥ 0.2, K-core ≥ 2, and max depth = 100 was set as the threshold value.

Prognostic Signature Generation and Validation
The TCGA original dataset was performed as a training cohort. Univariate and multivariate Cox proportional hazards regression analyses were carried out to identify potential genetic predictors for HNSCC survival. Kaplan-Meier survival analysis with logrank test was performed in R package. An internal dataset derived from the original TCGA served as a validation cohort using the bootstrap resampling method (27). Multivariate survival analysis was then performed to assess the association between the signature and clinical pathological index, namely, age, gender, lymphovascular invasion, margin status, recurrence, lymphatic metastasis, perineural invasion, cancer status, and nodal extracapsular spread.

Pathway Enrichment and Immune Enrichment Analyses
Gene-set enrichment analysis (GSEA) was performed using GSEA software with the criteria p < 0.05 and FDR < 0.25 (28,29) and visualized using clusterProfiler packages of R (30). mRNA expression profiles were uploaded to xCell online software to evaluate the immunocyte heterogeneity of LRPShigh and LRPS-low groups (31).

Statistical Analysis
Statistical analyses were performed using GraphPad Prism 8.4 and R software (version 3.6.3); p < 0.05 was considered statistically significant. A nonparametric t-test was performed to compare continuous variables and c 2 test was used to compare categorical variables between two groups. ANOVA test was utilized to compare more than two groups.

Lipid-Associated DEGs in HNSCC
To explore the lipid-related genes in HNSCC, 37 lipid-metabolic channels and 4 lipid-related signaling pathways (Supplementary Files S1) were selected based on KEGG pathway databases, and then analyzed in the TCGA database and GEO datasets using R packages. The result showed that a total of 65 genes significantly abnormally expressed in all three independent cohorts including TCGA, GSE30784, and GSE37991. The 26 upregulated and 39 downregulated lipid-related DEGs in total are listed in Table 1. The top 20 DEGs from the TCGA database are listed in Figure 1, and all lipid-related DEGs in GEO datasets are shown in Supplementary Figure S1 (p < 0.01, |log2FC| > 2).

Comprehensive Analysis of Molecular Characteristics in DEGs
Potential functions were then investigated in biological processes of the DEGs in HNSCC. GO analysis was performed and visualized in Figure 2A. The module of BP showed that the oxidation-reduction process, cholesterol homeostasis, sphingolipid biosynthetic process, and lipid metabolic and catabolic processes were commonly enriched. CC showed that the DEGs were significantly enriched in extracellular exosome, endoplasmic reticulum (membrane), and lipid particle. With regard to the module of molecular function (MF), the DEGs were mainly involved in iron ion binding and oxidoreductase activity.

Identification of a Lipid-Related Prognostic Signature of HNSCC
The LRPS was validated in the GSE65858 database. A total of 290 patients with HNSCC were subdivided into LRPS-high and -low groups according to the risk score, and survival analysis showed that the LRPS-high group had poorer 5-year overall survival than the LRPS-low group with a high effectivity ( Figures  S2A, B). We also established an internal TCGA dataset by bootstrap resampling method to validate the effectiveness of the LRPS (Supplementary Files S3, Validation dataset). The clinical characteristics within the two datasets had no significant differences using t-test analysis (Supplementary Table S1). The validation database was calculated and divided in the same way as the original group ( Figure S3). Five-year overall survival analysis demonstrated that the high-risk group (37.92%, 95% CI: 30.71%-46.8%) was significantly poorer than its counterpart (59%, 95% CI: 51.1%-68.2%, p < 0.001, Supplementary Files S5). Taken together, the lipid-based signature of ADCY2, LIPE, and OLR1 could effectively predict the HNSCC patients' survival.

LRPS Was an Independent Indicator of Prognosis and Closely Correlated to HNSCC Recurrence
Univariate Cox regression analysis showed that age, gender, lymphovascular invasion and metastasis, nodal extracapsular spread, perineural invasion, margin status, recurrence, cancer status, and LRPS score were significantly correlated with HNSCC prognosis ( Figure 5A). Multivariate Cox analysis determined LRPS as an independent predictor after adjustment by other pathologic characteristics ( Figure 5B). We evaluated the clinicopathologic factors in HNSCC among LRPS-high and LRPS-low groups ( Table 2). Meanwhile, the LRPS could also independently predict the overall survival of HNSCC patients from the GEO database ( Figures S2C, D). Statistically, LRPS-high patients were more prone to suffering relapse than the LRPS-low counterparts (52.73% vs. 22.92%, p = 0.0024). A high LRPS score had affinity relation with perineural invasion, compared with a low LRPS score (54.59% and 38.41% respectively, p = 0.0027). There were no significant differences within age, sex, alcohol and smoking history, tumor size and stage, lymphovascular invasion, and metastasis between the LRPS-high and LRPS-low groups. Furthermore, based on the primary tumor sites, the data were classified into four subgroups: oral cavity, tongue, larynx, and pharynx. The results showed that the proportion of oral cavity and larynx samples were almost equally distributed between the high risk and low risk, but there were more cases of tongue and fewer pharynx in the high-risk than in the low-risk group (p < 0.0001, Fisher test; Figure 6A). Meanwhile, the samples with HPV test results were subdivided into positive and negative groups, and the data were performed to analyze the association between HPV status and the lipid signature. Surprisingly, we found that almost all HPV-positive samples showed a low risk for LRPS, while HPV-negative samples had a high risk for LRPS (p < 0.0001, Fisher test; Figure 6B).

Molecular and Immune Characteristics in Different LRPS Subgroups
Since the lipid signature could increase the risk for recurrence, we sought to illuminate potential mechanisms regulating cancer relapse. Different LRPS groups were performed to GSEA and xCell analyses. We observed that the LRPS-high group significantly positively correlated to focal adhesion, MAPK signaling pathway, neuroactive ligand-receptor interaction, cancer-related pathway, and TGFb signaling pathway. The LRPS-low group was mainly enriched and negatively related to apoptosis, cell cycle, oxidative phosphorylation, p53 signaling pathway, and T-cell receptor signaling pathway ( Figure 7A, p < 0.05, FDR < 0.25). xCell analyses revealed that compared with adjacent normal tissues, the tumors that had a high LRPS score were more infiltrated in NKT, dendritic cells, monocytes, Treg, and M1 and M2 macrophages, which is in line with the inflammatory niche of the HNSCC tumor microenvironment. In addition, a proportion of B cells and CD4 + T effector cells including Th1 and Th2 significantly decreased in the LRPS-high group compared with the LRPS-low group, implicating that there is a suppressive immunity in the LRPS-high group ( Figure 7B). Notably, hematopoietic stem cells (HSCs) and mesenchymal stem cells (MSCs) were observed enriched in the LRPS-high group (p < 0.001). The above results indicated a possible changed immune milieu of primary tumor sites with an increased risk for HNSCC progression.
Finally, gene mutations were further analyzed to explore the molecular nature of the LRPS subgroups, and the top 10 genes with the highest mutation rates were identified ( Figure 8A). Genomic analysis showed that HNSCC endowed a high frequency of TP53 mutations, as high as about 71%, suggesting a vital role on tumor bioactivities. Our results showed a higher mutation rate of TP53 in LRPS-high patients than those with low LRPS score (82% vs. 60%), underlying a potential crosstalk between altered lipid metabolism and TP53 status. Missense variations were the most popular in both LRPS-high and -low groups. Importantly, TP53 showed a significantly higher mutation rate in the LRPS-high than in the LRPS-low group ( Figure 8B, p < 0.0001, Fisher's exact test). In addition, TTN, CDKN2A, FAT1, FRGB1, MUC16, CSMD3, PIK3CA, and SYNE1 were higher than 16% in both groups. The correlation between LRPS score and total mutation burden (TMB) was further explored, suggesting that the LRPS score was slightly correlated with total mutation burden (r = −0.11, p = 0.015, Figure S4).

DISCUSSION
The first new finding of the manuscript is that we identify the novel lipid prognostic signature of ADCY2, LIPE, and OLR1, which can predict the survival and prognosis of HNSCC patients as an independent effective prognostic factor. Meanwhile, our data may explain how lipidomics affects the prognosis and survival of patients with HNSCC through affecting tumor microenvironment via immunosuppression. In recent years, lipid metabolism has come into a sharp focus on cancer initiation and progression owing to its essential role in HNSCC and striking contribution to cancer development. In this study, we, for the first time, identified a novel LRPS of HNSCC through univariate and multivariate Cox analyses, which were performed to analyze the lipid-related DEGs as predictors for survival in TCGA patients with HNSCC. ADCY2, OLR1, and LIPE significantly predicted the overall survival of HNSCC  among the lipid DEGs. When the three genes were combined to indicate the prognosis of HNSCC patients, it showed that the LRPS-high group was highly related to poor prognosis. OLR1 is a stimulator of epithelial-mesenchymal transition (EMT) and involved in PPAR pathway, regulated by the secondary messenger cyclin adenosine monophosphate (cAMP). OLR1 promotes migration and metastatic spread in different pathways, such as TBC1D3/OLR1/TNFa/NF-kB, OLR1/c-Myc/HMGA2, oxLDL/OLR1/VEGF-C, and PI3K/Akt/ GSK3b (32)(33)(34)(35)(36). Recently, LOX-1D4, an alternative OLR1 isoform, has been shown to directly drive non-tumorigenic breast epithelial cells into fast proliferation status (37). More importantly, OLR1 is also reported to be positively correlated with the occurrence of lymphatic metastases in pancreatic cancers (38).
Adenylate cyclase 2 (ADCY2) encodes the adenylate cyclase that catalyzes ATP to transform to the second messenger cyclic adenosine monophosphate (cAMP). The latter is a crucial signal in cell fate, inflammation, and many other bioactivities, and is also greatly involved in the growth and differentiation of MSCs.
Zhao et al. have reported E2-induced ADCY2 as a positive regulator in MSCs (39). In colorectal cancer, ADCY2 could also be an important prognostic marker (40).
Lipase E, hormone-sensitive type (LIPE) increases both the levels of free cholesterol and free fatty acids, and plays an important role in adipocyte function and lipid and glucose homeostasis (41). More importantly, LIPE encodes the ratelimiting enzyme of lipolysis, and homozygous null mutation of LIPE results in marked inhibition of lipolysis, leading to multiple symmetric lipomatosis (42). Our studies showed that LIPE played a central role in the protein-protein interactions of the DEGs, significantly related to the survival rate of patients with HNSCC.
The second new finding of the study is that the LRPS can indicate the type of the infiltrated immune cells in the HNSCC tumor microenvironment. Comprehensive analyses indicated a diverse characteristic of LRPS subgroups. Lumps of the LRPShigh group showed a higher infiltration of inflammationassociated cells, including dendritic cells, M1 and M2 macrophages, and monocytes, yet a lower proportion of

immunocytes (B cells and pro B cells, CD8+
Tcm) compared with the LRPS-low group. The finding is also supported by the new concept that the infectious, chronic irritated, and inflammatory infiltration induces cancer and promotes neoplastic risk.
Macrophages are the main source of tissue repairmentrelated growth factors and cytokines after activation, such as TGFb1, TNFa, TGFa, and IL1 (43). These factors partly contributed to carcinogenesis via different signaling pathways. Numerous studies further showed a strong correlation between macrophage abundance and poor cancer prognosis, including thyroid cancer, lung cancer, and hepatocellular cancer (44)(45)(46)(47). Compromised immunity was also observed in our results, consistent with the research that shows that high-fat dietinduced obesity accelerates tumor growth by impairing CD8+ T-cell function (48). These observations could partly elucidate that the patients with LRPS high score were more subject to a poor survival because of the inflammatory-rich and immunodeficient conditions. Intriguingly, we found that LRPS-low harbors more Th1 and Th2 cells and fewer Treg cells in HNSCC. By contrast, LRPS-high has more Treg cells, consistent with the results from Whiteside group, which found a large number of Treg cells in the peripheral circulation of patients with HNSCC (49). Treg cells serve as one of the culprits that suppress anti-tumor immune response. Tumor within a niche of Treg cells is recognized as an unfavorable factor of cancer prognosis (50).
We further found that the patients with high LRPS score were more susceptible to recurrence because of increased infiltration of MSCs and HSCs in the tumor microenvironment. Recent data have proposed lipid metabolic rewiring as a new hallmark of cancer stem cells (CSCs) owing to its modification on stem-like cells' properties (51).
Taken together, these results affirmed that abnormal lipid metabolism exerts a great impact on immune cells' function in the tumor microenvironment, just influencing the progression and prognosis of HNSCC.
The third new finding of the study is that LPPS score can interpret the TP53 status of HNSCC. Our results showed that there were fewer LRPS-high samples in HPV-positive HNSCCs, which was in accordance with the negative relationship between HPV status and p53 mutation frequency (52). We also found a significant higher mutation rate of TP53 in LRPS-high patients than those with low LRPS score, underlying a potential crosstalk between altered lipid metabolism and TP53 status. Wild-type p53 supervises the cell damage response to various stimuli, and recent findings increasingly link p53 to lipid metabolism. P53 suppresses lipid biosynthesis via inhibiting lipogenesis, yet induces fatty acid oxidation as an alternative energy source to glycolysis in the condition of nutritional deficiency (53,54), implicating p53 as a positive regulator of catabolism (increase fatty acid levels) and an inhibitor of anabolism (decrease fatty acid levels) in the process of fatty acid metabolism. Otherwise, loss of p53 can lead to cell malignant transformation. As expected, mutated p53 exert a great impact on carcinogenesis through regulating gene transcription related to cell cycle, DNA repair, immunity and energetic activities, and so on. This gain of function of mutated p53 has been validated in various human cancers including breast, prostate, colon, pancreas, and head and neck cancers (55)(56)(57)(58)(59)(60). Our data further supported that p53 mutations did cooperate with abnormal lipid metabolism to promote cancer progression in HNSCC, though more laboratory investigations are needed in the future.
Though the LRPS has great potential for predicting HNSCC survival and p53 status, there are some limitations. The training and validation cohorts were retrospective, and more findings need to be validated prospectively. Moreover, the value of LRPS is not validated by in vitro and in vivo assays. Therefore, more studies are needed in the future.

CONCLUSION
Our data confirmed that the three lipid-related genes play a pivotal role in tumorigenesis and recurrence of HNSCC,  potentially by suppressing anti-tumor immunity and reflecting TP53 mutations status. LRPS has a potential to be a promising indicator of overall survival, prognosis, TP53 status, and immune characteristics in HNSCC, and perhaps could monitor and guide the treatment efficacy and prognosis of HNSCC in the future.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. These data can be found here:. The original datasets comparing the mRNA expression profiles between tumors and adjacent normal tissues were obtained from the GEO databases (GSE30784, GSE37991  File S1 | List of full lipid-related pathways and the corresponding genes from KEGG database. File S2 | List of protein-protein interaction (PPI) results from 65 lipid DEGs and two key modules within the PPI.
File S3 | The LRPS expression and the corresponding clinical information of highrisk and low-risk groups in the model dataset and the validation dataset, respectively.
File S4 | Results of all DEGs in HNSCC using univariate Cox analysis.
File S5 | Kaplan-Meier results of LRPS-high and -low subgroups in the model dataset and validation dataset, respectively.