Identification of genetic mechanisms underlying lipid metabolism-mediated tumor immunity in head and neck squamous cell carcinoma
BMC Medical Genomics volume 16, Article number: 110 (2023)
To identify the genetic mechanisms underlying lipid metabolism-mediated tumor immunity in head and neck squamous carcinoma (HNSC).
Materials and methods
RNA sequencing data and clinical characteristics of HNSC patients were procured from The Cancer Genome Atlas (TCGA) database. Lipid metabolism-related genes were collected from KEGG and MSigDB databases. Immune cells and immune-related genes were obtained from the TISIDB database. The differentially expressed genes (DEGs) in HNSC were identified and weighted correlation network analysis (WGCNA) was performed to identify the significant gene modules. Lasso regression analysis was performed to identify hub genes. The differential gene expression pattern, diagnostic values, relationships with clinical features, prognostic values, relationships with tumor mutation burden (TMB), and signaling pathways involved, were each investigated.
One thousand six hundred sixty-eight DEGs were identified as dysregulated between HNSC tumor samples and healthy control head and neck samples. WGCNA analysis and Lasso regression analysis identified 8 hub genes, including 3 immune-related genes (PLA2G2D, TNFAIP8L2 and CYP27A1) and 5 lipid metabolism-related genes (FOXP3, IL21R, ITGAL, TRAF1 and WIPF1). Except CYP27A1, the other hub genes were upregulated in HNSC as compared with healthy control samples, and a low expression of these hub genes indicated a higher risk of death in HNSC. Except PLA2G2D, all other hub genes were significantly and negatively related with TMB in HNSC. The hub genes were implicated in several immune-related signaling pathways including T cell receptor signaling, Th17 cell differentiation, and natural killer (NK) cell mediated cytotoxicity.
Three immune genes (PLA2G2D, TNFAIP8L2, and CYP27A1) and immune-related pathways (T cell receptor signaling, Th17 cell differentiation, and natural killer (NK) cell mediated cytotoxicity) were predicted to play significant roles in the lipid metabolism-mediated tumor immunity in HNSC.
Imbalance in lipid metabolism homeostasis is an important metabolic hallmark of cancer [1, 2]. Water-insoluble molecules, e.g., triacylglycerides, phosphoglycerides, sterols and sphingolipids, formed by lipids are principal components of biological membranes, and critical building blocks used for energy storage. In cancer cells, lipid metabolism and related processes are critical for energy generation along with the maintenance of membrane components and signaling molecules necessary for various cancer cell processes such as cell proliferation, survival, invasion, metastasis [1, 2]. The impairment of lipid metabolic homeostasis not only influences the cellular processes of cancer cells, but also play critical roles in regulating anti-tumor immunity . By secreting lipid-derived signaling molecules in the tumor microenvironment (TME), tumor cells can alter the metabolic programming of immune cells toward fatty acid oxidation, which is associated with pro-tumorigenic and immune-suppressive phenotypes of immune cells .
Head and neck squamous cell carcinoma (HNSCC) accounts for more than 90% of all head and neck cancers worldwide, ranking sixth in the world . A recent reported that 23 genes associated with lipid metabolism showed significant prognostic values in HNSCC, where 11 genes (ARSI, CYP27B1, CYP2D6, DGKG, DHCR7, LPIN1, PHYH, PIP5K1B, PLA2G2D, RDH16, and TRIB3) were associated with clinicopathological features (stage, gender, and pathological stage) . However, this study  did not investigate the underlying genetic and immunological mechanisms involved in lipid metabolism-mediated regulation of tumor immunity in HNSC. Recent work has demonstrated that HNSC tumor cell-derived fatty acids can mediate the maturation of macrophages into a M2 phenotype which plays anti-inflammatory and pro-tumor role . Another recent study found that lipid nanoparticle-based targeted delivery led to a decrease in the adenosine A2A receptor, which in turn increased chemotaxis and T cell infiltration . Considering the emerging evidence regarding a pivotal role of lipid metabolism in regulating tumor immunity in HNSC, an improved understanding of the molecular mechanisms involved could provide basis for novel therapeutic strategies to overcome tumor-induced immunosuppression in HNSC.
The advances in computational biology approaches have resulted in several analytical approaches that could generate insights into molecular processes involved in lipid metabolism-mediated regulation of tumor immunity in HNSC. In the present study, we leveraged weighted gene co-expression network analysis (WGCNA), Lasso regression analysis, survival analysis, nomogram plot analysis, tumor mutation burden (TMB) analysis, and gene-pathway interaction analysis for this purpose. By comprehensively applying these analyses, the current study aimed to investigate the most significant mechanisms involved in lipid metabolism-mediated tumor immunity in HNSC.
Materials and methods
Procurement of HNSC datasets
We downloaded RNA-seq data of HNSC samples from TCGA-GDC (https://portal.gdc.cancer.gov/) , and the extracted gene expression values as a TPM (Transcripts Per Million) dataset. We selected the samples with sample numbers beginning with 01 and 11 for analysis, which was based on the sample type codes in the TCGA code table (URL: https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables/sample-type-codes). The samples with sample numbers starting with 01 comprised the Tumor group (such as TCGA-BA-5151-01A), and those with sample numbers starting with 11 comprised the Normal group (such as TCGA-CV-7406-11A). We also downloaded the relevant clinical information data and the SNV (simple nucleotide variation) dataset , where the data type of the SNV dataset was ‘Masked Somatic Mutation’.
Procurement of lipid metabolism-associated genes and immune-related gene set
The lipid metabolism-related feature gene set was gathered from KEGG (Kyoto Encyclopedia of Genes and Genomes, http://www.kegg.jp/blastkoala/)  and MSigDB (Molecular Signatures Database v7.5.1, https://www.gsea-msigdb.org/gsea/msigdb/index.jsp)  using the following keywords; fatty acyls, glycerolipids, glycerophospholipids, sphingolipids, sterol lipids, prenol lipids, saccharolipids, and polyketides. We thus obtained 21 pathways related to lipid metabolism from the KEGG database, and 6 gene sets related to lipid metabolism from MSigDB database. After integrating the data from the two databases, a total of 1079 lipid metabolism-related genes and 27 pathways were finally obtained. We downloaded immune cells and related genes from the TISDB database (http://cis.hku.hk/TISIDB/download.php) , and obtained 28 immune cells and 782 immune-related genes.
We selected tumor samples that contained clinical information, and then used the selected tumor samples and normal samples to form an expression matrix for subsequent analysis. We de-duplicated the duplicate genes in the expression matrix according to the mean, deleted the genes whose expression values were 0 in more than 50% of the samples, and performed log2 scaling on the gene expression values.
Gene set variation analysis
Gene set variation analysis (GSVA) is a non-parametric and unsupervised method for identifying patterns in gene expression profiles that correlate with biological pathways . We first downloaded datasets of pathways and genes from the MSigDB database . We then merged the lipid metabolism-related genes and immune-related genes, extracted their expression values in HNSC samples. And performed GSVA enrichment analysis on the gene expression matrix to obtain the enrichment scores of pathways in each sample. Next, we applied the “limma” package  in R to perform differential expression analysis  for Tumor vs Normal sample groups. The thresholds used included |logFC|> 0.25 and P.adjust < 0.05, where log2FC > 0.25 indicated differentially up-regulated pathways, and log2FC < 0.25 indicated differentially down-regulated pathways.
Univariate analysis to screen survival-related genes
We merged lipid metabolism-related genes and immune-related genes, and then extracted the gene expression values for the combined gene set in the Tumor group. The ‘survival’ package  in R was used to build a Cox proportional hazards regression model (Cox-PH)  for each gene. The genes identified by the Cox-PH univariate analysis with P-value < 0.05 were regarded as survival-related genes and included in the subsequent analysis.
WGCNA for the significant module
The expression values of the survival-related genes in the tumor case group were obtained. The WGCNA package  in R was applied to establish a scale-free weight network for these genes and samples, and the consistency module was extracted. An appropriate adjacency matrix weight parameter β value was selected by setting the β value from 1 to 30, and then calculating the corresponding correlation coefficient and gene adjacency function mean of the dataset. The higher the correlation coefficient (R2) (maximum = 1), the closer the network is to a no network scale distribution. At the same time it is also necessary to ensure a certain degree of gene connectivity.
Based on the correlation coefficient and gene connectivity, we selected the β value, and then applied the TOMsimilarity method  to establish a given adjacency matrix based on the expression matrix. Next, we used the ‘dynamic cut tree’ algorithm to cut the given adjacency matrix and perform module mining and thereafter used the ‘mergeCloseModules’ method to merge modules with correlation coefficients greater than 0.8. We then obtained the P value of the genes in the univariate analysis and looked at the significance of each module. Among the modules obtained by WGCNA, the grey module comprised an unassigned genes module and was not included in the subsequent analysis. We counted the number of genes in other significant modules, and selected modules with both lipid metabolism-related genes and immune-related genes for the subsequent analysis.
Hub gene screening in significant modules
Through WGCNA analysis, we obtained significant modules related to lipid metabolism genes and immune genes. We then extracted the expression values of these genes in the module in tumor and normal samples and performed differential expression analysis using limma in R  with the comparison method as Tumor vs Normal. The resultant genes with |logFC|> 0.5 and P.adjust < 0.01 were considered as the significant differently expressed genes (DEGs). Next, we used LASSO (Least absolute shrinkage and selection operator) logistic regression analysis  to further screen the DEGs. Pearson’s correlation coefficient analysis was applied to investigate the correlation between the lipid metabolism-related genes and immune-related genes. If a lipid metabolism-related gene was highly correlated with all immune-related genes, then we extracted this gene and recorded it as a ‘Hub’ gene. We extracted the expression levels of the Hub genes in tumor and normal samples and applied the Wilcoxon test to the Hub genes in different groups. Next, we applied ROC analysis  to test the prognostic values of the hub genes.
DNA methylation analysis of immune-related hub genes
DNA methylation is currently the most studied epigenetic modification and is essential for promoting important biological processes such as embryonic development, genomic imprinting and X chromosome inactivation. MethSurv (URL: https://biit.cs.ut.ee/methsurv/) is an intuitive web-based tool for multivariable survival analysis based on CpG methylation patterns . The MethSurv tool was used to analyze methylation level differences in the immune-related hub genes between groups of samples with different clinical characteristics. Using the common regions of the immune-related hub genes, the relationship between these genes’ methylation level and sample clinical characteristics was analyzed. In this analysis, DNA methylation values were represented using beta values (ranging from 0 to 1) . The parameters used in the methylation analysis are listed in Table 1. In addition, MethSurv was also used to analyze the influence of immune-related hub genes’ DNA methylation on the survival prognosis of TCGA-HNSC .
Immune cell infiltration analysis for immune-related hub genes
The expression values of the immune-related hub genes in HNSC tumor samples were collected and immune infiltration analysis was performed using Tumor Immune Estimation Resource (TIMER, URL: http://timer.cistrome.org/) database . TIMER analyzes the enrichment scores of immune cells in different samples using 6 computational methods (e.g., CIBERSORT, CIBERSORT-ABS, EPIC, MCPCOUNTER, QUANTISEQ, and TIMER) [26,27,28]. The scores of the immune cells using the different methods were obtained, and then their relationship with immune cells was analyzed using Pearson’s correlation coefficient analysis . An absolute r value between 0 and 0.3 indicated a weak linear relationship, that between 0.3 and 0.7 indicated a moderate relationship, and that between 0.7 and 1.0 indicated a strong relationship .
Prediction of targeted small molecule drugs
The significant module associated with survival identified by the WGCNA analysis was used and a connectivity Map (CMap, URL: https://clue.io/)  was used to discover small molecule targets and functionally annotate genetic variants of disease genes . Both upregulated genes (n ≥ 10) and downregulated genes (n ≥ 10) were uploaded to the CMap database to predict potential small molecular drug targets for the genes in the significant module. The top 10 drugs/molecules with positive normalized WTCS (weighted connectivity score) value and the top 10 drugs/molecules with negative normalized WTCS value were obtained using CMap .
Multivariate analysis of hub genes
In order to investigate the relationship between the hub genes and survival outcomes of HNSC, we extracted the expression values of hub genes in tumor samples, and a Cox-PH model was built to predict OS and OS_Event using multivariate analysis. Thus, the risk scores for the hub genes were obtained. The samples were divided into high-risk and low-risk groups on based on the median risk scores and we examined whether low-risk survival differed at different time periods (3 years, 5 years, and 10 years). Next, we extracted the expression values of the hub genes in different risk groups and performed Wilcoxon’s test for different risk groups.
To analyze the relationship between clinical characteristics and high- and low- risk groups, the clinical characteristics and sample risk scores were integrated. We then constructed a nomogram plot by utilizing the nomogram method in the “rms” package  in R. The Wilcoxon’s test was applied to analyze whether there were significant differences in risk scores between different groups of clinical characteristics.
Somatic mutation analysis of hub gene
Based on SNV data, we analyzed the mutation information regarding hub genes in the tumor samples and calculated the TMB score  of the tumor samples. We extracted the expression values of hub genes in HNSC samples and used the TMB scores for Pearson correlation coefficient analysis  to analyze the relationship between hub genes and TMB.
Pathway analysis of hub genes
We obtained the genes related to all pathways from the KEGG database , and then extracted the pathways where the hub genes were represented. Next, we extracted all the genes under these pathways, and selected survival-related lipid metabolism genes and immune genes. Thereafter, we used Cytoscape (version 3.9.1) software  to construct a hub genes-pathways-survival genes interaction network.
We obtained the HNSC dataset from TCGA and obtained 501 tumor samples and 44 normal samples. Table 2 summarizes the clinical features of the 501 HNSC tumor samples.
Gene set variation analysis
Fourty three common genes (Fig. 1A) were found overlapped between 1079 lipid metabolism-related genes, 782 immune-related genes, and 1688 HNSC-DEGs. These 1688 HNSC-DEGs regulated the T cell receptor signaling pathway, cytosolic dna sensing pathway, ecm receptor interaction, nod like receptor signaling pathway, toll like receptor signaling pathway, and metabolism of xenobiotics by cytochrome p450 among others (Fig. 1B).
Screening of functional genes associated with survival by univariate analysis
Three hundred seventy-nine genes related to survival, including 192 immune-related genes and 181 lipid metabolism-related genes and 6 common genes (Fig. 2).
Mining modules highly related to lipid metabolism- and immune-related genes
In the process of WGCNA, we set the β value from 1 to 30, and then calculated the corresponding correlation coefficient and gene adjacency function mean. At a β of 24, the established network was found closest to the scale-free network (Fig. 3A–B). By using this β value as the network construction parameter, we built the WCGNA model. Using the dynamic cut tree algorithm to mine modules, we set at least 5 genes in each module, and the maximum connection height was 0.99 (minModuleSize = 5, cutHeight = 0.95). As a result, 5 modules were obtained including the grey, yellow, brown, blue, and turquoise modules (Fig. 3C).
After obtaining the modules, we merged the modules with correlation coefficients greater than 0.8, and obtained 4 modules, among which the brown module had the highest significance (Fig. 3D). We extracted the genes in the modules and the P values for univariate analysis of these genes, based on the 'mean (P value)' statistical module significance (Fig. 3E). Then, we counted the gene groupings in brown, blue, and yellow modules (Fig. 3F), and found that the blue module contained genes related to lipid metabolism and immunity, while other modules only contained lipid metabolism-related genes. Therefore, the blue module was regarded as the module of interest and used for the subsequent analyses.
Hub gene screening in significant modules
Based on the thresholds of P.adjust < 0.01 and |logFC|> 0.5 for significantly differentially expressed genes, a total of 34 DEGs were obtained. The LASSO regression analysis was used to remove redundant genes from the 34 DEGs (Fig. 4A–B) and obtained 19 differentially expressed genes (Table 3). The 19 genes included 16 lipid metabolism-related genes, 2 immune-related genes (PLA2G2D and TNFAIP8L2) and one common gene (CYP27A1). Figure 4C shows that 3 immune-related genes (PLA2G2D, TNFAIP8L2 and CYP27A1) were highly correlated with FOXP3, IL21R, ITGAL, TRAF1 and WIPF1 (cor > 0.5) (Fig. 4C). Therefore, these 3 immune-related genes (PLA2G2D, TNFAIP8L2 and CYP27A1) and 5 highly related lipid metabolism-related genes (FOXP3, IL21R, ITGAL, TRAF1 and WIPF1) were considered as hub genes.
The relationship between the methylation levels of the 3 immune-related genes and clinical characteristics of HNSC
The methylation levels of CYP27A1 were higher than that of TNFAIP8L2 and CYP27A1, irrespective of clinical feature analyzed. Figure 5A shows that the methylation level of CYP27A1 in patients younger than 61 years old was higher compared to that in patients older than 61 years and the lowest methylation level of CYP27A1 was observed in patients older than 69 years. The methylation level of PLA2G2D in patients older than 69 years was higher compared with that among patients with the age ranges 65–61 years and 61–69 years. The methylation expression level of TNFAIP8L2 in patients older than 61 years was lower compared to that of patients younger than 61 years old and the lowest methylation expression level of TNFAIP8L2 was observed in patients older than 69 years (Fig. 5).
The methylation level of CYP27A1 in African American patients was lower compared to that in Asian patients and White patients. The methylation level of PLA2G2D in White patients was higher than that in Asian and African American patients. The methylation level of TNFAIP8L2 in Asian patients was higher compared to that for African American patients and White patients (Fig. 5B).
No significant differences in the methylation level of CYP27A1 and TNFAIP8L2 were noted between female and male patients, while the methylation level of PLA2G2D in male patients was lower than that in female patients (Fig. 5C).
The highest methylation level of CYP27A1 was seen in patients with clinical stage I; while the lowest methylation expression of CYP27A1 was seen in patients with clinical stage III. The highest methylation level of PLA2G2D was seen in patients with clinical stage I; while the lowest methylation level of PLA2G2D was seen in patients of clinical stage IV. The highest methylation level of TNFAIP8L2 was seen in patients of clinical stage III, while the lowest methylation expression of TNFAIP8L2 was seen in patients of clinical stage II (Fig. 5D).
Kaplan–Meier (Fig. 5E) revealed that methylation levels of CYP27A1 (P = 0.43), TNFAIP8L2 (P = 0.71), and CYP27A1 (P = 0.53) were not significantly associated with the prognosis of HNSC patients.
Correlation between 3 immune-related hub genes and TIICs
A negative correlation with the highest |r| value was observed between TNFAIP8L2 and Macrophage M1 (r = − 0.51, P = 1.56E−34), identified by the CIBERSORT method. A positive correlation with the highest |r| value was observed between TNFAIP8L2 and Myeloid dendritic cell_TIMER (r = 0.85, P = 2.361–138). The lowest P value for correlation was also observed between TNFAIP8L2 and Myeloid dendritic cell_TIMER (r = 0.85, P = 2.361–138) (Fig. 6, Table 4).
The targeting relationship between immune-related hub genes and small molecule drugs
Table 5 shows the top 10 positive and negative scoring small molecules with normalized WTCS (Weighted Connectivity Score) values selected as the predicted target small molecules. The small molecule drug targeted by PLA2G2D was mepacrine. No small molecule drug was predicted for the other two immune-related hub genes (TNFAIP8L2 and CYP27A1).
The expression pattern of hub genes
A heatmap was used to show the expression values of the eight hub genes in the HNSC tumor and normal samples. CYP27A1 was downregulated in HNSC tumor samples, while the other genes were upregulated in HNSC tumor samples (Fig. 7A) A heatmap (Fig. 7B) depicted the expression levels related to different clinical features.
The hub genes were significantly dysregulated between tumor and healthy control samples (Fig. 8A). The overall expression level of CYP27A1 was higher than that of the other hub genes. The ROC analysis showed five lipid metabolism-related hub genes (FOXP3, IL21R, ITGAL, TRAF1 and WIPF1) had higher predictive values than three immune-related hub genes (PLA2G2D, TNFAIP8L2 and CYP27A1) (Fig. 8B).
Grouping tumor samples of HNSC based on multivariate analysis
We extracted the expression values of the hub genes in the tumor case samples, and established a Cox-PH model based on the clinical characteristics of OS and OS_Event for multivariate analysis. The results showed that the expression patterns of 8 hub genes were significantly related with the overall survival outcomes of HNSC patients (Fig. 9A). The samples were divided into high-risk and low-risk groups based on the median of their risk scores. The number of non-surviving samples in the high-risk group was higher than that in the low-risk group (Fig. 9B). The survival probability of high-risk and low-risk groups within 3 years was not significantly different (P = 0.37) (Fig. 10A). However, there were significant differences between the high-risk group and the low-risk group in the other three time periods (i.e., 5 year (P = 0.017), 10 year (P = 0.0078), and 20 year (P = 0.0035), showing that the survival rate of the high-risk group was significantly lower than that of the low-risk group (Fig. 10B–D). These findings indicated that as the survival time of HNSC patients increases, lipid metabolism-related genes and immune-related genes have a greater impact on survival. The 8 hub genes were expressed differentially between high- and low-risk groups, with the high-risk group expressing each hub gene at a lower level compared with the low-risk group (Fig. 11).
Association of clinical characteristics and risk groups with survival
A nomogram of 467 samples was plotted to demonstrate the effect of different clinical characteristics and risk scores on survival (Fig. 12). Figure 13 shows the differences in risk scores among subjects grouped by clinical characteristics. Figure 14 shows that the patients with clinical characteristics age < 60, gender = male, stage = (IV, T3, N0, N2)) showed significant differences in terms of overall survival outcomes between high- and low- risk groups, and the survival probability of patients in the high risk group was lower than that in the low risk group.
The somatic mutation of hub genes
ITGAL and WIPF1 were the most mutated in tumor samples (Fig. 15A). TNFAIP8L2, CYP27A1, FOXP3, IL21R, ITGAL, TRAF1, WIPF1 and TMB were significantly correlated; however, PLA2G2D was not significantly correlated with TMB (Figs. 15B, 13I).
Pathway network analysis of hub genes
Figure 16 shows that CYP27A1 and CYP27B1 were involved in metabolic pathways. Both CYP27A1 and CYP27B1 are genes related to lipid metabolism and immunity. The lipid metabolism-related genes CYP27A1, ACOX3, ADCY1 and ACAA1, are also involved in the regulatory pathway PPAR signaling pathway. The hub genes co-regulate cancer-related pathways by directly or indirectly interacting with other lipid metabolism genes or immune-related genes.
The main findings of the current research identified 3 immune-related genes (PLA2G2D, TNFAIP8L2, and CYP27A1) which are likely to play highly significant roles in lipid metabolism-mediated tumor immunity regulation in HNSC. PLA2G2D (Phospholipase A2, Group IID) encodes the secreted phospholipase A2 enzyme that lyses glycerophospholipids to yield lysophospholipids and free fatty acids . PLA2G2A is induced and secreted by immune cells including dendritic cells and macrophages when stimulated by proinflammatory mediators, and also results in augmented production of proinflammatory lipid and cytokines by acting on phospholipids . Dendritic cells in lymphoid tissues abundantly express PLA2G2D, which regulates the steady-state levels of anti-inflammatory lipids in order to resolve Th1 immune response . Chronic inflammation induced by PLA2G2A could contribute to the promotion of carcinogenesis. Miki et al. found that the overexpression of PLA2G2D suppressed anti-tumor immunity by increasing tumor-promoting M2-like macrophages and decreasing tumor-suppressing M1-like macrophages and cytotoxic T cells in skin carcinoma . Contrary to the findings of the current study, another study reported that PLA2G2D expression was positively correlated with immune infiltration and indicated better prognosis in HNSC patients . The knockdown of PLA2G6 (Phospholipase A2 Group VI) was found to suppress the progression of melanoma by affecting the phenotypes of melanoma cells, inhibiting cell proliferation, migration, and invasion, promoting tumor cell apoptosis . To our knowledge, experimental research investigating the influence of PLA2G2D on the phenotypes of HNSC cells is lacking. Of note, the small molecule drug targeted by PLA2G2D mepacrine (also named as quinacrine), which was originally used as an anti-malarial drug has been recently repurposed as an anticancer agent in treating ovarian cancer, gynecologic, breast cancer, and colon cancer [39, 40]. Mepacrine was also found to inhibit cell viability and clonogenic survival, as well as promote apoptosis of five types of HNSC cell lines (CAL27, SCC040, FaDu, SCC47 and VU147) . Mepacrine induced TP53 mRNA and protein expression, increased TP53 reporter activity and p21 protein expression, and induced growth inhibition in the wild-type TP53 HNSC cell lines .
TNFAIP8L2 (TNF Alpha Induced Protein 8 Like 2; also named as TIPE2) is a lipid transfer protein and can inhibit lipid biosynthesis pathways by negatively regulating lipid biosynthesis-related gene signatures . TNFAIP8L2 acts as a negative regulator of innate and adaptive immunity by inhibiting the function of Toll-like receptor and T-cell receptor . TNFAIP8L2 was found to positively regulate and enhance the anti-tumor immune response in head and neck cancer . A previous research comprehensively assessed the expression and prognosis of TNFAIPs family members in head and neck cancer and found a positive correlation between TNFAIP8 and tumor infiltrating immune cells macrophage, neutrophil, CD8+ T cell, CD4+ T cell, and dendritic cell) . In addition, the overexpression of TNFAIP8L2 predicted improved overall survival in head and neck cancer patients . TNFAIP8L2 was also found significantly associated with cancer stem cell index, indicating its potential role as a novel immune checkpoint gene for the immunotherapy of cancers . The transfection of TNFAIP8L2 in hepatocellular carcinoma (HCC) cell lines markedly inhibited tumor cell growth, migration and invasion in vitro , but research investigating the effects of TNFAIP8L2 upregulation on the phenotypes of HNSC cells is lacking.
CYP27A1 (Cytochrome P450 Family 27 Subfamily A Member 1) is a key enzyme involved in the process of bile acid synthesis and involved in regulating cellular cholesterol homeostasis. Mutations in CYP27A1 result in reduced bile acid synthesis, increased production of cholestanol, and subsequently cholestanol accumulation . The dysregulation of CYP27A1 expression has been found to be a prognostic biomarker in breast cancer , prostate cancer , and ovarian cancer . CYP27A1 was found to be highly expressed in myeloid immune cells and macrophages and played a pro-tumorigenic role in breast cancer by impairing T cell expansion . However, the current research indicated that a high expression of CYP27A1 indicated the low risk of death in HNSC. In agreement, CYP27A1 was found to inhibit proliferation and migration of clear cell renal cell carcinoma by activating the LXRs/ABCA1 pathway . Experimental research investigating the involvement and regulatory roles of CYP27A1 in the phenotypes of HNSC tumor cells is warranted.
The current study showed several immune cell-related signaling pathways were involved in lipid metabolism-mediated tumor immunity in HNSC, including T cell receptor signaling, Th17 cell differentiation, and natural killer (NK) cell mediated cytotoxicity. Several studies have highlighted these relationships. Modulating fatty acid metabolism can enhance CD8 T-cell memory generation, thereby further amplifying the anti-tumor immunity . A high content of cholesterol in tumors inhibits anti-tumor immunity by upregulating immune checkpoint genes and further inducing CD8+ T Cell exhaustion . The cholesterol esterification enzyme acetyl-CoA acetyltransferase (ACAT1) knockout in CD8+ T cells was found to downregulate membrane cholesterol and improve T cell receptor clustering and signaling, enhancing CD8+ T cell function and anti-tumor immunity in mouse tumor models . Th17 lymphocytes express a proinflammatory phenotype by secreting cytokines (e.g., IL-10, IL-17, IFN-γ, and IL-22)  and the Th17/Treg balance relies heavily on fatty acid metabolism . Short-chain fatty acids (SCFAs) activate and promote the differentiation of Th17 lymphocytes into T regulatory (Treg) cells, which play immunosuppressive functions in the tumor microenvironment . The innate lymphoid cells-NK cells play a crucial role in preventing tumor metastasis . NK cells in the lipid-rich environment are found to be immature and defective in the ability to lyse target tumor cells . The cytotoxicity of NK cells is impaired by lipid accumulation in cancer patients, which contributes to an immunosuppressive tumor microenvironment . Although these immune-related pathways have been implicated in lipid metabolism-mediated cancer biology, evidence in the context of head and neck cancer is lacking and merits future investigation.
This study utilized a computational biology approach to identify candidate immune-related genes and pathways involved in imbalanced lipid metabolic homeostasis-mediated tumor immunity in HNSC. Using a multi-level approach, the most significant candidates were identified, and the findings provide a theoretical foundation for future discovery of lipid metabolism-targeting therapies in HNSC. The candidate genes and mechanisms identified in this study need to be validated by designing relevant experiments. We suggest that co-culture environment models of lipid laden-HNSC cells and immune cells be established, and the genetic reciprocal interactions between HNSC cells and immune cells be examined. In vivo experiments could be also designed to observe the effects of lipid metabolism alterations on the immune function of tumor-bearing animal models. The current study was mainly focused on the effects of abnormal lipid metabolism on individual tumor infiltrating immune cells, which is far from sufficient as a particular lipid metabolic pathway might produce contradictory results for different types of immune cells . Therefore, specific mechanisms underpinning lipid metabolic homeostasis imbalance in HNSC tumor immune regulation merit comprehensive research.
The current research identified immune genes PLA2G2D, TNFAIP8L2, and CYP27A1 and immune cells-related pathways T cell receptor signaling, Th17 cell differentiation, and NK cell mediated cytotoxicity as the primary candidates involved in lipid metabolism-mediated tumor immunity in HNSC. These may comprise promising therapeutic biomarkers and targets for modulation in the field of cancer immunotherapy for head and neck cancer.
Availability of data and materials
The data analyzed during the current study are available in the TCGA database with the accession numbers TCGA-HNSC. The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.
Bian X, Liu R, Meng Y, Xing D, Xu D, Lu Z. Lipid metabolism and cancer. J Exp Med. 2021;218:1.
Santos CR, Schulze A. Lipid metabolism in cancer. FEBS J. 2012;279:2610–23.
Liu X, Zhang P, Xu J, Lv G, Li Y. Lipid metabolism in tumor microenvironment: novel therapeutic targets. Cancer Cell Int. 2022;22:1–13.
Broadfield LA, Pane AA, Talebi A, Swinnen JV, Fendt S-M. Lipid metabolism in cancer: new perspectives and emerging mechanisms. Dev Cell. 2021;56:1363–93.
Klein JD, Grandis JR. The molecular pathogenesis of head and neck cancer. Cancer Biol Ther. 2010;9:1–7.
Xiong Y, Si Y, Feng Y, Zhuo S, Cui B, Zhang Z. Prognostic value of lipid metabolism-related genes in head and neck squamous cell carcinoma. Immun Inflamm Dis. 2021;9:196–209.
Albakri MM, Huang SC-C, Tashkandi HN, Sieg SF. Fatty acids secreted from head and neck cancer induce M2-like macrophages. J Leukoc Biol. 2022;112:617.
Newton HS, Chimote AA, Arnold MJ, Wise-Draper TM, Conforti L. Targeted knockdown of the adenosine A2A receptor by lipid NPs rescues the chemotaxis of head and neck cancer memory T cells. Mol Ther-Methods Clin Dev. 2021;21:133–43.
Jensen MA, Ferretti V, Grossman RL, Staudt LM. The NCI genomic data commons as an engine for precision medicine. Blood J Am Soc Hematol. 2017;130:453–9.
Katsonis P, Koire A, Wilson SJ, Hsu T-K, Lua RC, Wilkins AD, et al. Single nucleotide variations: biological impact and theoretical interpretation. Protein Sci. 2014;23:1650–66.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28:27–30.
Liberzon A, Subramanian A, Pinchback R, Thorvaldsdóttir H, Tamayo P, Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics. 2011;27:1739–40.
Wan J, Qian S-B. TISdb: a database for alternative translation initiation in mammalian cells. Nucleic Acids Res. 2014;42:D845–50.
Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf. 2013;14:1–15.
Smyth GK. Limma: linear models for microarray data. In: Bioinformatics and computational biology solutions using R and Bioconductor. Springer; 2005. p. 397–420.
Anders S, Huber W. Differential expression analysis for sequence count data. Nat Preced. 2010;1–1.
Therneau TM, Lumley T. Package ‘survival.’ R Top Doc. 2015;128:28–33.
Harrell FE. Cox proportional hazards regression model. In: Regression modeling strategies. Springer; 2015. p. 475–519.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinf. 2008;9:1–13.
Shuai M, Chen X. Algorithm optimization for weighted gene co-expression network analysis: accelerating the calculation of Topology Overlap Matrices with OpenMP and SQLite. bioRxiv. 2021.
Weppler S, Schinkel C, Kirkby C, Smith W. Lasso logistic regression to derive workflow-specific algorithm performance requirements as demonstrated for head and neck cancer deformable image registration in adaptive radiation therapy. Phys Med Biol. 2020;65: 195013.
Fawcett T. An introduction to ROC analysis. Pattern Recognit Lett. 2006;27:861–74.
Modhukur V, Iljasenko T, Metsalu T, Lokk K, Laisk-Podar T, Vilo J. MethSurv: a web tool to perform multivariable survival analysis using DNA methylation data. Epigenomics. 2018;10:277–88.
Anuraga G, Wang W-J, Phan NN, An Ton NT, Ta HDK, Berenice Prayugo F, et al. Potential prognostic biomarkers of NIMA (Never in Mitosis, Gene A)-Related Kinase (NEK) family members in breast cancer. J Pers Med. 2021;11:1089.
Song Y, Ma R. Identifying the potential roles of PBX4 in human cancers based on integrative analysis. Biomolecules. 2022;12:822.
Li T, Fu J, Zeng Z, Cohen D, Li J, Chen Q, et al. TIMER2. 0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res. 2020;48:W509–14.
Kao T-J, Wu C-C, Phan NN, Liu Y-H, Ta HDK, Anuraga G, et al. Prognoses and genomic analyses of proteasome 26S subunit, ATPase (PSMC) family genes in clinical breast cancer. Aging. 2021;13:17970.
Laham AJ, El-Awady R, Lebrun J-J, Ayad MS. A bioinformatics evaluation of the role of dual-specificity tyrosine-regulated kinases in colorectal cancer. Cancers. 2022;14:2034.
Benesty J, Chen J, Huang Y, Cohen I. Pearson correlation coefficient. In: Noise reduction in speech processing. Springer; 2009. p. 1–4.
Musa A, Ghoraie LS, Zhang S-D, Glazko G, Yli-Harja O, Dehmer M, et al. A review of connectivity map and computational approaches in pharmacogenomics. Brief Bioinform. 2018;19:506–23.
Subramanian A, Narayan R, Corsello SM, Peck DD, Natoli TE, Lu X, et al. A next generation connectivity map: L1000 platform and the first 1,000,000 profiles. Cell. 2017;171:1437–52.
Wang C-Y, Chiao C-C, Phan NN, Li C-Y, Sun Z-D, Jiang J-Z, et al. Gene signatures and potential therapeutic targets of amino acid metabolism in estrogen receptor-positive breast cancer. Am J Cancer Res. 2020;10:95.
Harrell FE Jr, Harrell MFE Jr, Hmisc D. Package ‘rms.’ Vanderbilt Univ. 2017;229:Q8.
Kristensen K, Nielsen A, Berg CW, Skaug H, Bell B. TMB: automatic differentiation and Laplace approximation. ArXiv Prepr ArXiv150900660. 2015.
Smoot ME, Ono K, Ruscheinski J, Wang P-L, Ideker T. Cytoscape 2.8: new features for data integration and network visualization. Bioinformatics. 2011;27:431–2.
Miki Y, Kidoguchi Y, Sato M, Taketomi Y, Taya C, Muramatsu K, et al. Dual roles of group IID phospholipase A2 in inflammation and cancer. J Biol Chem. 2016;291:15588–601.
Murakami M, Yamamoto K, Miki Y, Murase R, Sato H, Taketomi Y. The roles of the secreted phospholipase A2 gene family in immunology. Adv Immunol. 2016;132:91–134.
Wang Y, Song H, Miao Q, Wang Y, Qi J, Xu X, et al. PLA2G6 silencing suppresses melanoma progression and affects ferroptosis revealed by quantitative proteomics. Front Oncol. 2022;12.
Oien DB, Pathoulas CL, Ray U, Thirusangu P, Kalogera E, Shridhar V. Repurposing quinacrine for treatment-refractory cancer. In: Seminars in Cancer Biology. Elsevier; 2021. p. 21–30.
Kumar M, Sarkar A. Repurposing of anti-malarial drug quinacrine for cancer treatment: a review. Sci Pharm. 2022;90:12.
Bryant J, Batis N, Franke AC, Clancey G, Hartley M, Ryan G, et al. Repurposed quinacrine synergizes with cisplatin, reducing the effective dose required for treatment of head and neck squamous cell carcinoma. Oncotarget. 2019;10:5229.
Friedman J, Nottingham L, Duggal P, Pernas FG, Yan B, Yang XP, et al. Deficient TP53 expression, function, and cisplatin sensitivity are restored by quinacrine in head and neck cancer. Clin Cancer Res. 2007;13:6568–78.
Li T, Wang W, Gong S, Sun H, Zhang H, Yang A-G, et al. Genome-wide analysis reveals TNFAIP8L2 as an immune checkpoint regulator of inflammation and metabolism. Mol Immunol. 2018;99:154–62.
Sun H, Gong S, Carmody RJ, Hilliard A, Li L, Sun J, et al. TIPE2, a negative regulator of innate and adaptive immunity that maintains immune homeostasis. Cell. 2008;133:415–26.
Lan G, Yu X, Sun X, Li W, Zhao Y, Lan J, et al. Comprehensive analysis of the expression and prognosis for TNFAIPs in head and neck cancer. Sci Rep. 2021;11:1–12.
Bai K-H, Zhang Y-Y, Li X-P, Tian X-P, Pan M-M, Wang D-W, et al. Comprehensive analysis of tumor necrosis factor-α-inducible protein 8-like 2 (TIPE2): A potential novel pan-cancer immune checkpoint. Comput Struct Biotechnol J. 2022;20:5226.
Cao X, Zhang L, Shi Y, Sun Y, Dai S, Guo C, et al. Human tumor necrosis factor (TNF)-alpha-induced protein 8-like 2 suppresses hepatocellular carcinoma metastasis through inhibiting Rac1. Mol Cancer. 2013;12:1–10.
Lorbek G, Lewinska M, Rozman D. Cytochrome P450s in the synthesis of cholesterol and bile acids–from mouse models to human diseases. FEBS J. 2012;279:1516–33.
Kimbung S, Inasu M, Stålhammar T, Nodin B, Elebro K, Tryggvadottir H, et al. CYP27A1 expression is associated with risk of late lethal estrogen receptor-positive breast cancer in postmenopausal patients. Breast Cancer Res. 2020;22:1–13.
Alfaqih MA, Nelson ER, Liu W, Safi R, Jasper JS, Macias E, et al. CYP27A1 loss dysregulates cholesterol homeostasis in prostate cancer CYP27A1 loss is involved in prostate cancer progression. Cancer Res. 2017;77:1662–73.
He S, Ma L, Baek AE, Vardanyan A, Vembar V, Chen JJ, et al. Host CYP27A1 expression is essential for ovarian cancer progression. Endocr Relat Cancer. 2019;26:659–75.
Ma L, Wang L, Nelson AT, Han C, He S, Henn MA, et al. 27-Hydroxycholesterol acts on myeloid immune cells to induce T cell dysfunction, promoting breast cancer progression. Cancer Lett. 2020;493:266–83.
Liang Z, Jiao W, Wang L, Chen Y, Li D, Zhang Z, et al. CYP27A1 inhibits proliferation and migration of clear cell renal cell carcinoma via activation of LXRs/ABCA1. Exp Cell Res. 2022;419: 113279.
Pearce EL, Walsh MC, Cejas PJ, Harms GM, Shen H, Wang L-S, et al. Enhancing CD8 T-cell memory by modulating fatty acid metabolism. Nature. 2009;460:103–7.
Ma X, Bi E, Lu Y, Su P, Huang C, Liu L, et al. Cholesterol induces CD8+ T cell exhaustion in the tumor microenvironment. Cell Metab. 2019;30:143–56.
Yang W, Bai Y, Xiong Y, Zhang J, Chen S, Zheng X, et al. Potentiating the antitumour response of CD8+ T cells by modulating cholesterol metabolism. Nature. 2016;531:651–5.
Hubler MJ, Kennedy AJ. Role of lipids in the metabolism and activation of immune cells. J Nutr Biochem. 2016;34:1–7.
Endo Y, Kanno T, Nakajima T. Fatty acid metabolism in T-cell function and differentiation. Int Immunol. 2022;34:579.
Ichise H, Tsukamoto S, Hirashima T, Konishi Y, Oki C, Tsukiji S, et al. Functional visualization of NK cell-mediated killing of metastatic single tumor cells. Elife. 2022;11:e76269.
Niavarani SR, Lawson C, Bakos O, Boudaud M, Batenchuk C, Rouleau S, et al. Lipid accumulation impairs natural killer cell cytotoxicity and tumor control in the postoperative period. BMC Cancer. 2019;19:1–14.
Yu W, Lei Q, Yang L, Qin G, Liu S, Wang D, et al. Contradictory roles of lipid metabolism in immune response within the tumor microenvironment. J Hematol OncolJ Hematol Oncol. 2021;14:1–19.
Nothing to declare.
This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Liu, S., Wang, S. & Wang, Z. Identification of genetic mechanisms underlying lipid metabolism-mediated tumor immunity in head and neck squamous cell carcinoma. BMC Med Genomics 16, 110 (2023). https://doi.org/10.1186/s12920-023-01543-6
- Lipid metabolism
- Tumor immunity
- Head and neck squamous carcinoma