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.

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.

The GEO data were processed and analyzed using GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/). TCGA mRNA counts were normalized and analyzed by R packages (DESeq2 package) (p < 0.01, |log2FC| > 2).

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.

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.

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 log-rank 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.

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 LRPS-high and LRPS-low groups (31).

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 χ ² test was used to compare categorical variables between two groups. ANOVA test was utilized to compare more than two groups.

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).

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.

Next, KEGG pathway enrichment analysis was used to figure out functions of the proteins encoded by the DEGs. As shown in Figure 2B , the DEGs were closely associated with arachidonic acid metabolism, PPAR signaling pathway, regulation of lipolysis in adipocytes, and metabolic pathways. p < 10⁻⁵ was recognized as significantly enrichment categories.

To figure out the relationship between the DEGs in HNSCC, the PPI network was constructed by STRING online database ( Figure 2C ). The central node genes (more than 10 connections or interactions) and the top 10 highly connected genes were identified, namely, PPARG, LIPE, SLC27A6, CYP2E1, ADIPOQ, PLA2G16, PLIN1, PLA2G2A, CYP2J2, and SLC2A4 ( Supplementary Files S2 ). MCODE plugin from Cytoscape was used for the key module within the PPI network. The two most significant modules were identified for further pathway enrichment analysis. Module 1 consisted of 13 hub genes, namely, LIPE, PPARG, ADIPOQ, SLC2A1, SLC2A4, SLC27A6, MGLL, PLIN1, PLIN4, FABP3, CYP2E1, ALOX12, and CYP2J2 ( Figure 2D ). Module 2 included 5 hub genes, namely, SQLE, DHCR7, HMGCS2, TM7SF2, and CH25H ( Figure 2E ). Pathway enrichment analysis revealed that the hub genes in module 1 were closely correlated with PPAR signaling pathway (p = 2.58 × 10⁻⁶) and AMPK signaling pathway (p = 8.29 × 10⁻⁴). Module 2 was mainly enriched in steroid biosynthesis (p = 4.80 × 10⁻⁵) ( Supplementary Files S2 ).

To verify whether the lipid-related DEGs could be potential prognostic markers for HNSCC, the univariate and multivariate Cox analyses were performed to analyze the lipid-related DEGs as predictors for survival in TCGA patients with HNSCC ( Supplementary Files S3 , Model dataset). Univariate Cox analysis showed that ADCY2, OLR1, and LIPE significantly affected the overall survival of patients with HNSCC among the DEGs ( Figures 3A, B and Supplementary Files S4 ). Next, a lipid-related prognostic signature (LRPS), containing LIPE, ADCY2, and OLR1, was constructed based on the coefficient of multivariate Cox analysis and mRNA expression of the three genes, the risk score = (−0.15) × LIPE + (0.08) × ADCY2 + (0.09) × OLR1 ( Figure 3C ). Subsequently, the LRPS containing LIPE, ADCY2, and OLR1 was selected to predict the prognosis of HNSCC patients through the TCGA and GEO databases.

The patients were accordingly divided into the high-risk group (n = 249) and low-risk group (n = 249) based on the risk score ( Figures 4A, B ), finding that the HNSCC patients in the high-risk group had a poorer 5-year overall survival (36.9%, HR = 0.377, 95% CI = 29.8%-45.7%) than those in the low-risk group (55.9%, HR = 0.566, 95% CI = 47.07%–66.5%) (p = 4.889 × 10⁻⁶) ( Figure 4C and Supplementary Files S5 ). The concordance indices (C-index) of the lipid signature showed a higher specificity and sensitivity for predicting 3-, 5-, and 10-year overall survival (C-index = 0.645, 0.592, and 0.66, respectively, Figure 4D ).

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.

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 ).

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 TGFβ 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 ).