Fatty acid metabolism-related signature predicts survival in patients with clear cell renal carcinoma

Objective: To explore fatty acid metabolism-related genes and signature, which could predict survival outcomes of clear cell renal carcinoma (ccRCC) patients. Materials and Methods: Transcriptional and survival data of fatty acid genes in ccRCC patients were retrieved from UCSC Xena and Geo DataSets. We first performed Lasso Cox regression analysis to identify survival-related genes. These genes were then used to construct metabolic-related gene signature and risk score. Enrichment analysis and immune component and chemotherapy response prediction were also performed. Results: In total, five survival-related genes were identified: AGR2, HAO2, IGF2BP1, MCCD1 and OLFM4 (p < 0.05). A series of survival value analyses revealed survival-related signature and risk score, including KM analysis (training set: p < 0.001; test set: p = 0.008). Four clinical indexes (T stage, N stage, M stage, and pathology) were positively correlated with risk score. Time-dependent ROC analysis yielded AUC value of 0.813. Immune landscape analysis revealed that risk score was strongly correlated with TAM score and cytotoxic score. Patients with high risk score and TAM score or cytotoxic score had the shortest survival time. Finally, inhibition of fatty acid metabolism in human ccRCC cell line produced corresponding changes in five genes, consistent with our preliminary results. Conclusion: We identified five survival-related genes (AGR2, HAO2, IGF2BP1, MCCD1 and OLFM4) in ccRCC patients. Our results also indicated that survival-related signature based on these genes is a potential robust prognostic biomarker for ccRCC in patients.


INTRODUCTION
Self-sufficiency in growth, tissue invasion and metastasis are typical tumor characteristics [1]. To support these proliferative activities, tumor tissue must undergo metabolic reprogramming. A classic example of metabolic reprogramming was the Warburg's effect in which glucose undergoes anaerobic glycolysis in tumor cells [2]. In recent years, several studies have shown that a few metabolic substrates, including glucose, lactic acid, fatty acid, and amino acid, are preferentially utilized for cancer growth, invasion and metastasis. Lipid metabolism has also been reported to perform a role in cancer metastasis. Previously understood as simple metabolic substrates, fatty acids are now known to influence lipid level in tissues, becoming a subject of much research interest [3].
Clear cell renal carcinoma is a urological malignant tumor with high incidence in the population. However, effective treatment for advanced renal clear carcinoma in patients is still lacking. Developing new molecular drug targets for treating advanced renal clear carcinoma would potentially improve survival outcomes of patients. Based on previous literature, key fatty acid proteins were reported to be strongly correlated with renal clear cell carcinoma's growth, invasion and metastasis [4,5]. Inhibition of fatty acid synthase was shown to activate P53 and STAT pathways, decreasing growth of renal cancer cells [6]. Moreover, fatty acid binding proteins (FABP) have been implicated in renal cell carcinogenesis [7].
However, fatty acid metabolism-related genes set associated with clear cell renal carcinoma in patients have not been systematically studied. In this study, we investigated key fatty acid-related protein and potential signature predicting survival of renal clear cell carcinoma patients to provide a basis for future research directions.

Data source
Two clear cell renal carcinoma patients' datasets with survival results were used in this study: TCGA-KIRC dataset, including RNA sequences and clinical phenotypes, downloaded from the UCSC Xena [8] and GSE29609 dataset, including expression matrix and survival results, downloaded from the GEO DataSets [9]. These two datasets were publicly available from the databases. We used TCGA-KIRC dataset as a training set to perform a series of survival analyses and identify a survival-related signature. In contrast, GSE29609 was designed as a test set for testing the signature's accuracy and stability.

Survival-related fatty acid genes
Based on a review of previous literature, 90 fatty acidrelated genes were selected [10]. To find the survival value of each of the 90 fatty acid genes, we performed Lasso cox regression analysis on the TCGA-KIRC dataset. Subsequently, univariate and multivariate Cox regressions were performed on genes with survival value to find survival-related fatty acid genes.

Signature of survival-related fatty acid genes
We predicted a new signature containing survivalrelated fatty acid genes. The following formula was used to calculate each sample risk score: i Risk Score = (Coefi × ExpGenei)  "Coef" represents non-zero regression coefficients calculated with univariate Cox regression analysis. "ExpGene" represents expression values of genes from the prognostic risk score model. We calculated risk score level for TCGA-KIRC dataset' patients and verified the survival value of the signature in both the training set and test set using KM analysis.

Relationship between risk score and clinical features
First, using test set, we performed a time-dependent ROC analysis to predict risk score on 1-, 3-and 5-year survival results. Second, important clinical prognostic features, such as age and TMN stage, were compared in the training set using Cox regression analysis. R package "survival", and "glmnt" were used for ROC and Cox regression analyses. We constructed a forest plot showing these results. In addition, prognostic characteristics with significantly different risk scores selected for further analysis. A nomogram of overall survival for clear cell renal carcinoma patients in the training set with risk scores and other important clinical prognostic indexes was drawn.

Enrichment analysis for low-and high-risk score patient groups
Patients included in the training set were divided into two groups: high risk-score group and low risk-score group. R package 'limma' was used to identify differentially expressed genes between two patient groups. Genes with FDR <0.05 were considered statistically significant. Three functional enrichment analyses were performed on differentially expressed genes for the two groups: GO analysis [11], KEGG analysis [12], and Gene set variation analysis (GSVA) [13]. R packages "cluster Profile" and "GSVA" were used for enrichment analysis. Based on enrichment analysis results, we identified candidate different functions between the high-score group and low-score group patients.

Immune-related characteristics in low-and high-risk score groups
We performed multi-dimension immune-related analysis on the two groups of patients. First, we used cibersort website (https://cibersort.stanford.edu/) to analyze the immune component of the two patients groups and divided patients into immune subtypes [14]. Second, we used R package "GSVA" to analyze immune-related score. Third, we explored immune status of the two groups using R package "estimate". Fourth, the TAM score and Cytotoxic score, important prognostic indexes generated using R package "GSVA", were compared between the two patient groups. Finally, we compared the expression levels of two important immune checkpoint genes: CTLA4 and PDL1.

Chemotherapy drug sensitivity analysis
Ridge regression model was constructed with R package "pRRophetic" and "oncoPredict" using the remaining data. Drug sensitivity (IC50) was chosen as an outcome index. Chemotherapeutic response was predicted using tumor genes expression and drug sensitivity data of cell lines from the cancer genome project (CGP) and drug sensitivity in cancer (GDSC) [15].

Verification of survival-related fatty acid genes in cell lines
Increasing the expression level of fatty acid synthase (FASN) was positively correlated with aggressive cell proliferation, migration, apoptosis, and lipid droplet formation. It was also shown to regulate metabolic disorders associated with clear cell renal carcinoma.
Researchers also proved that pharmacological inhibitor of FANS could suppress the growth and invasiveness of clear cell renal carcinoma. Therefore, we used two human clear cell renal carcinoma cell line (caki-1 and 786-0 cell lines) from authenticated cell cultures of the Chinese national collection. We down-regulated the expression of FANS in caik-1 cells by the way of siRNA virus and exposed FASN inhibitorC75 (HY-12364, MedChemExpress, China) with two cell line (caki-1 and 786-0). We verified the expression level of survival-related fatty acid genes with quantitative realtime polymerase chain reaction (Q-PCR) [16].

Data availability statement
The original data presented in this study are included in the article. All data were retrieved from public databases.

Identity of survival-related fatty acid genes
We downloaded 90 fatty acid genes, 525 patient clinical records and RNA sequences from the TCGA-KIRC dataset (the training set) and survival results of 39 renal cancer patients and the matrix of genes' expression from the GSE29609 (the test set). We performed Lasso regression analysis on the training set for all fatty acid genes ( Figure 1A, 1B). The results showed that only 14 genes had survival value for renal cancer patients. After univariate and multivariate Cox regression analyses ( Figure 1D), only five survival-related genes (AGR2, HAO2, IGF2BP1, MCCD1 and OLFM4) were significantly different to have the statistical value (p < 0.05). The correlation among the five genes was also established in the Figure 1C.

Risk score construction and survival value verification
The risk score was calculated using the following formula: The patients in the training set were divided into lowand high-risk score groups based on their risk score level. We also presented the distribution of risk score in the training set ( Figure 1E). To explore the survival value of each risk score, we tested it on both training and test sets using KM analysis. The results showed that renal cancer patients with high risk score had shorter survival time for both the training set (p < 0.001, Figure 1F) and test set (p = 0.008, Figure 1G). We also used time-dependent ROC curve to analyze the risk score in the test set. The AUC value of risk score was 0.813, with satisfactory results for ROC curves (Figure 2A, 2B).

Relationship between risk score and clinical features
We compared clinical features (like TMN stages) in the training set using Cox regression analysis. The results showed that six clinical features (age, T stage, N stage, M stage, gender and pathology stage) had significant survival value (p < 0.05, Figure 2C) and were presented in a forest plot. A further separate comparison of these six indexes between the two groups showed that patients with high-risk score had worse T stage, N stage, M stage, and pathology stage (P < 0.001) than low-risk score patients ( Figure 2D-2I). A nomogram containing several clinical features and risk score was drawn predicting 1-, 3-and 5-year overall survival ( Figure 3A-3D). Using R package "limma", we identified differentially expressed genes, which were represent in a volcano plot ( Figure 4A). Three forms of functional enrichment analysis were used to determine potential function of risk score in two groups of renal cancer patients. First, GO analysis showed that functions of extracellular matrix organization, extracellular structure organization, and urogenital system development would be activated in high-risk score patients ( Figure 4C). In contrast, the functions of negative thymic T cell selection and actin cytoskeleton reorganization were inhibited in high-risk score patients. Second, KEGG analysis showed that pathways of PIK3, papillomavirus infection, and MAPK signaling would be activated in high-risk score patients ( Figure 4D). Third, GSVA analysis ultimately found that a high-risk score was strongly correlated with inflammatory response and epithelial mesenchymal response in patients ( Figure 4B).

Immune-related characteristics in low-and high-risk score groups and chemotherapy drug sensitivity analysis
Using cibersort website, we identified immune components of two patient groups. However, differences between main immune cells could not be identified by scanning whole immune landscapes of the two groups ( Figure 5A). However, high-risk score patient group were more likely to exhibit immune insensitivity subtypes, such as immune C2 and C5, than low-risk score patients. Subsequently, we used R package "estimate" to analyze the estimate score and immune score ( Figure 5B, 5C). The results showed that immune scores were significantly different between the two patient groups (p = 0.043). Scanning the immune scores with R package "GSVA" revealed a strong correlation of TAM score and cytotoxic score were with risk score. TAM score and cytotoxic score have been reported to have a prognostic value for many cancers.
Here, we analyzed the two scores for renal cancer patients. The results showed a strong positive survival value of TAM score and cytotoxic score in renal cancer patients ( Figure 5D). KM analysis further revealed that patients with high-risk score and TAM score or cytotoxic score had poorer survival results than low-risk score patients (p < 0.001, Figure 5E, 5F). We further compared the correlation of risk score with two important immune therapy genes: PDL1 and CTLA4.
The results showed a strong positive correlation between risk score and expression level of PDL1 or CTLA4 ( Figure 6B, 6C). Finally, we performed chemotherapy drug sensitivity analysis on CGP and GDSC platforms. The landscape of drug sensitivity analysis was presented in Figure 6A.

Verification of survival-related fatty acid genes in caki-1 cell line
We verified survival related fatty acid genes in caki-1 cell and 786-0 cell (two human renal clear cell carcinoma cell lines) in two ways: down-regulating FASN expression and exposure to pharmacological inhibitor of FASN. We measured expression levels of five survival-related fatty acid genes. The results showed that three genes (AGR2, IGF2BP1 and OLFM4) had been inhibited whereas two genes (HAO2 and MCCD1) had been activated in both two ways ( Figure 6D, 6E). These changes suggested that the five genes regulated fatty acid synthase function in renal cancer patients. Therefore, the five-gene signature could be used to predict survival outcomes of renal cancer patients.

DISCUSSION
Lipid metabolism has attracted significant attention in cancer research in recent years because it provides a novel connection between blockers of lipid metabolism and inhibitors of tumor growth and invasion [17]. FASN regulates fatty acid metabolism and plays a critical role in the growth and tissue invasion by many cancers [18,19]. We found a few studies providing deep analysis of the role of FASN in clear cell renal carcinoma development. The study by Xu and team found that FASN expression was positively correlated with aggressive cell proliferation, migration, apoptosis and lipid droplet formation and regulated metabolic disorders in clear cell renal carcinoma microenvironment [6]. Further, using a pharmacological inhibitor of fatty acid synthase has been shown to suppress growth and invasiveness of renal cancer cells [20].
With increasing research evidence for the role of fatty acids in clear cell renal carcinoma development in patients, analysis of fatty acid genes was warranted in the present study. With this purpose, we conducted a series of analyses to identify these genes. In total, we identified five survival-related fatty acid genes implicated in ccRCC. Consistent with these results, HAO2 was reported to inhibit malignancy of clear cell renal cell carcinoma by promoting lipid catabolic process [21]. In addition, AGR2 and IGF2BP1 promoted tumorigenesis by accelerating hypoxia state in renal cell carcinoma cell line [22][23].
However, no direct metabolic function could be identified for the other three genes (MCCD1, and OLFM4) in research. We used two ways to inhibit fatty acid metabolism in clear cell renal carcinoma cell line. The five survival-related genes exhibited corresponding expression variation upon Q-PCR test. Consistent with these results, previous studies showed that inhibition of fatty acid metabolism in a cell line using the two methods could suppress cancer cell line growth or invasion. Our results illustrated that the expression levels of the five genes not only reflected the fatty acid   [24][25].
Metabolic reprogramming usually transforms immune ability and immune components in cancer tissues [2]. Interdisciplinary research examining metabolic reprogramming and immune function is increasingly pursued in cancer studies. We used three methods to identify the risk score, immune ability and immune component: cibersort website, estimate and GSVA. We found no significant differences in the immune component between low-risk and high-risk score groups using cibersort. In contrast, using estimate, we found stronger immune ability in the high-risk score group than in the low-risk score group. In addition, we analyzed different immune-related scores to determine TAM score and cytotoxic score. Our results showed that patients with high-risk score and TAM score or cytotoxic score had lower survival time than those with low-risk score. Recent studies on immune cells have found a correlation between TAM cells were and cancer progression and anti-immunotherapy [26]. So far, few researchers have studied the connection between the fatty acid metabolism and percentage or function of TAM cells in the cancer patients. Our results suggested a potential connection, providing a novel direction for future studies. Cytotoxic score, which has not been clearly defined in literature, was positively correlated with risk score or TAM score in clear cell renal carcinoma patients. Therefore, cytotoxic score provided a new approach to predicting survival outcomes of clear cell renal carcinoma patients.
Studies examining fatty acid metabolism, a major metabolism reprogramming mechanism in cancer patients, have elicited much research interest in the topic in past several years. Whereas FASN has been a subject of intense research, the other fatty acid-related genes have not been studied for their role in cancer. The multifarious functions of fatty acid proteins in cancer provide another research direction that needs further exploration.

CONCLUSION
We identified five survival-related genes (AGR2, HAO2, IGF2BP1, MCCD1 and OLFM4) in clear cell renal carcinoma patients. The survival signature based on these genes was proved to have a significant prognostic value for clear cell renal carcinoma in patients.

AUTHOR CONTRIBUTIONS
WRJ wrote the draft manuscript, and SJW did the experiment and edited the manuscript; CY, and GJG analyzed the data; YJX developed the images. All authors revised the manuscript and approved the submitted version.

CONFLICTS OF INTEREST
The authors declare no conflicts of interest related to this study.

ETHICAL STATEMENT AND CONSENT
Ethical review and approval were waived in accordance with the local legislation and institutional requirements. Written informed consent for participation was not required for this study in accordance with the national legislation and institutional requirements.

FUNDING
This study was financially supported by the medicine and health program of Zhejiang province (2021KY345); Huzhou Science and Technology Bureau Project (2021GY34), and Zhejiang public welfare technology project (LGF22H160022).