Construction of the prognostic model
A total of 309 genes related with fatty acids metabolism were included in this study. After TCGA datasets filtering, quality assessment and data processing, 96 DEGs were finally extracted through the “limma” R package and the results indicated that 62 DEGs were downregulated and 34 DEGs were upregulated. The heatmap and volcano map was used to visualize DEGs (Figures 1b, 1c). The prognostic model was developed based on univariate and multivariate Cox regression analysis as well as LASSO analysis (Figures 1d, 1e). Ultimately, the prognostic model was visualized by “forest” (Figure 1f) consisting of 20 genes: HACD1, HPGD, ALOX15B, ABCD1, HMGCS2, CPT1B, TDO2, SCD5, PCCA, DPEP1, ALAD, ACADM, ACADSB, ACAT1, PLA2G4A, ALOX12B, IL4I1, ACAD11, HIBCH, LTC4S. Risk score= (0.319240 * HACD1) + (-0.084391 * HPGD) + (-0.180228 * ALOX15B) + (0.459381 * ABCD1) + (-0.063391 * HMGCS2) + (0.325758 * CPT1B) + (0.140669 * TDO2) + (-0.029671 * SCD5) +(0.127550 * PCCA) + (-0.067707 * DPEP1) + (-0.342020 * ALAD) + (-0.040588 * ACADM) + (0.131207 * ACADSB) + (-0.110134 * ACAT1) + (-0.161595 * PLA2G4A) + (-0.407284 * ALOX12B) + (0.280250 * IL4I1) + (-0.028668 * ACAD11) + (-0.158444 * HIBCH) + (-0.766512 * LTC4S).
Prognostic model validation
The ROC curves and the KM curves were utilized to explore the predictive power of the prognostic model on the internal and external test sets. The area under the ROC curve (AUC) of the 3-, 5-, and 7-year in TCGA-KIRC train and test cohort as well as E-MTAB-1980 cohort presented good prediction performance (Figures 2a-2c). Using the same prognostic model, we classified the remaining patients in all the test sets into different risk groups based on the median of all risk scores. The risk score of the prognostic model in TCGA-KIRC train and test as well as E-MTAB-1980 cohort was an independent protective factor of overall survival (OS). Patients with high-risk scores exhibited significantly lower OS compared to those in the low-risk group (P<0.05, Figures 2d, 2e). K-M analysis further revealed that ccRCC patients with high-risk score had lower progression-free-survival (PFS) in TCGA-KIRC cohort (P<0.05, Figures S1a-S1d).
The clinicopathological characteristics in TCGA-KIRC cohort
As shown in Figure 3a, a clinical correlation analysis was performed, demonstrating that grade and stage levels increased with the risk score of the prognostic model (P<0.05, Figure 3a). Moreover, the risk score exhibited superior predictive power, with the largest area under the ROC curve compared to age, gender, grade, and stage (Figure 3b). Furthermore, the AUC values of the 3-, 5-, and 7-year in TCGA-KIRC total cohort were 0.745, 0.765, and 0.742 respectively (Figure 3c). Subsequently, based on the clinical features and risk score, a nomogram for OS prediction was conducted and composed of age, stage as well as risk score as the independent prognostic factors (P<0.05) (see Figure 3d). In the calibration diagram (Figure 3e), the 1-, 3-, and 5-year OS for ccRCC individuals had a good predictive performance of this personalized nomogram model (C-index = 0.779).
The immune landscape of the prognostic model and response to anti-PD-1/ PD-L1 therapy
The response to anti-PD-1/ PD-L1 therapy was examined based on the immune-related cohorts (CheckMate-025, IMmotion150 and IMmotion151), and results indicated that different risk groups between SD/PD and CR/PR possessed a statistically significant difference. Not surprisingly, results from the immune-related databases revealed that the high-risk group processes a significantly lower response to anti-PD-1/ PD-L1 therapy (Figure 3f). Further, through the TIDE algorithm, high-risk group in the TCGA-KIRC database had a high potential of immune escape with a worse effect of immunotherapy (Figure 3g). Applying the ESTIMATE algorithm, we calculated the overall levels of the immune score. The immune score of the high-risk group was higher than low-risk group (Figures 3h).Moreover, as expected, in IMmotion150 and IMmotion151 cohorts, the low-risk group for survival was linked to higher PFS (P<0.05, Figure 3i). We then performed a hierarchically K-M analysis of the CheckMate-025 study. The results indicated there was a significant difference in the survival rate between different risk groups in the ccRCC patients with Nivolumab monotherapy from the CheckMate-025 cohort while there was no significant difference in those treated with Everolimus (a mammalian target of mTOR inhibitor) (Figures S1d,1e), suggesting that our model is more effective and suitable for anti-PD-1/ PD-L1 therapy in the ccRCC patients. Subsequently, in the comparison of several immune signatures, inflammation-promoting, T cell co-stimulation, checkpoint, antigen-presenting cell (APC) co-stimulation, chemokine receptors (CCR), and Type I IFN response were found to be higher in the high-risk group than the low-risk group, while type II IFN response was obviously downregulated (see Figure 3j). Additionally, TME cell composition and fraction of individual immune cell types in three Immune-related cohorts were computed and generated with CIBERSORT (http://cibersort.stanford.edu/). We observed that the high-risk groups presented a higher fraction of M0 macrophages than the low-risk group (P<0.05, see Figures S2a-2c), thus highlighting the risk score of our model could predict clinical response to ICI-based immunotherapy.
Risk scores and drug sensitivity analysis
We conducted drug sensitivity analysis based on the different risk groups through the “pRRophetic” package and the IC50 value was utilized to measure the sensitivity to drugs (Figure 4a). The high-risk group possessed significantly higher IC50 values but was less sensitive to the drugs than the group with low-risk (P<0.05). For quantifying the individual patients, the correlations between risk score of prognostic model and targeted drugs was assessed through the “ggplot” package and identified that Sorafenib, Erlotinib, Saracatinib and Crizotinib had high degree of correlation with risk score of prognostic model, suggesting that the constructed model could effectively predict efficacy and sensitivity to chemotherapy (Figures 4b). We further evaluated the risk score of ccRCC cell lines through CCLE and drug sensitivity data from GDSC. The results indicated that A498 had the highest risk score while BFTC-909 had the lowest risk score, which implied that higher malignancy existed in A498 compared with other cell lines (Figure S3a). Pearson correlation analysis suggested that the risk score had a positive correlation with IC50 of C-75 (an inhibitor of fatty acids synthase), while was negatively correlated with IC50 of PD0325901 (an oral potent ERK inhibitor) (Figure S3b). These revealed that the higher risk score was more sensitive to PD0325901 and less to C-75, which might be an effective targeted therapy of ccRCC.
Identifying related pathways and immune checkpoints
We utilized ssGSEA to examine the correlations between the risk score and the enrichment scores of related pathways or immune checkpoints to explore the immune-related functional processes. Our findings revealed a positive correlation between the risk score and JAK-STAT3 signaling, as well as major immune checkpoints (Figures 4c, 4d). These implied that fatty acids metabolism in ccRCC was associated with JAK-STAT3 signaling and the risk score levels could reflect the therapeutic effect of ccRCC patients treated with anti-PD-1/ PD-L1.
Enrichment analysis and GSEA hallmark visualization
To confirm the underlying mechanism of JAK-STAT3 signaling and fatty acids metabolism in ccRCC progression, we performed GSEA on TCGA-KIRC cohort. The results based on GSEA analysis indicated that differentially expressed target genes are enriched in immune and metabolism-related functional pathways. The upregulated DEGs groups, in particular, were enriched in JAK-STAT3 signaling. Meanwhile, fatty acids metabolism was observed in the downregulated DEGs groups (see Figure 4e). KEGG enrichment analysis was developed for tumorigenic pathways enrichment analysis in TCGA-KIRC cohort which demonstrated an association with cytokine-cytokine receptor interaction, fatty acids metabolism and complement and coagulation cascades (Figure 4f). GO analysis was performed to further explore the potential biological processes of DEGs and revealed that most immune responses and associated activities were significantly enriched in these genes (Figure 4g). Thus, these results indicated that the fatty acids metabolism might contribute to ccRCC development mainly concentrated on immune response.
Evaluation of joint indicators and multi model comparison
To establish the superiority of our model, we conducted a comparison of the accuracy of several immune indicators and previous studies in the TCGA-KIRC cohort. [20]. The ROC curve indicated that the AUC of the risk score is significantly higher than other indicators and other fatty acids metabolism-related model, suggesting that our model was more accurate(Figure 4h). Overall, the constructed model demonstrated greater representativeness in fatty acids metabolism compared to other network models or indicators.Validation of mRNA expression
To analyze the mRNA expression profiles, we explored the expression levels of 6 fatty acids metabolism-related genes in our risk model in tumor and normal tissues as well as ccRCC cell lines. The qPCR results indicated that the expression levels of ABCD1, ALOX12B, ALOX15B, HACD1, IL4I1 significantly upregulated in tumor samples, while the expression of CPT1B was low in RCC tissues (Figures 5a, 5b). Similarly, the expression patterns of these six genes were also observed in ccRCC cell lines (Figure 5c). Additionally, oligo sequences in the qPCR were displayed in Supplementary Table S2.
Silencing IL4I1 suppresses the growth and invasion of ccRCC cells
To elucidate the role of IL4I1 in ccRCC migration and invasion in vitro, three siRNA specifically targeting IL4I1 (si-IL4I1-1, si-IL4I1-2 and si-IL4I1-3) were constructed respectively (Figure 5d). Transwell assays revealed that the migration and invasion capabilities of 786-O and 769-P were inhibited after the silence of IL4I1 (Figure 6a). Furthermore, in accordance with the above results, colony formation assays and wound healing demonstrated that silencing IL4I1 significantly suppresses the proliferation of 786-O and 769-P (Figures 6b, 6c). Altogether these results collectively elucidated that IL4I1 could promote the growth and invasion of ccRCC cells.
Silencing IL4I1 impacts JAK1/STAT3 signaling pathway
Based on the above analysis about our model, GSEA suggested that our model strongly associated with JAK1/STAT3 signaling pathway (Figure 5e). To further clarify the effects of IL4I1 on JAK1/STAT3 signaling pathway, western blot was assessed to test IL4I1’s contribution to the functional outcomes of JAK1/STAT3 signaling pathway. The result showed that silencing IL4I1 downgrades the expression of phosphorylated JAK1 and phosphorylated STAT3 (Figure 5f), suggesting that IL4I1 could modulate JAK1/STAT3 signaling pathway and lead to JAK1/STAT3 phosphorylation.
Silencing IL4I1 inhibits M2-like macrophages polarization
To investigate whether IL4I1 signaling to Chemokine ligand 2 (CCL2) could mediate macrophages polarization, an indirect co-culture condition between si-IL4I1-transfected 786-O and 769-P as well as M0-like macrophages was conducted. Flow cytometry results revealed that silencing IL4I1 inhibited M2-like macrophages polarization (Figures 7b,7c). Additionally, western blot also implied that silencing IL4I1 could suppress the level of CD206 (M2 macrophage surface marker) (see Figure 7d) and even CCL2 expression of si-IL4I1-transfected ccRCC cells (786-O and 769-P) (Figure 5f). Moreover, IHC staining validated that the high expression of IL4I1 and CD206 in tumor samples were simultaneously higher than that in adjacent normal tissues (Figure 7a). In summary, these results proved that the knockdown of IL4I1 in si-IL4I1-transfected 786-O and 769-P inhibited M2-like macrophage polarization engaged by the regulation of CCL2.