Construction of a prognostic model based on FA acid metabolism-related lncRNAs
To construct FA metabolism-related lncRNA signatures in HCC, 158 genes were extracted from the Hallmark FA metabolism pathway and a single sample GESA were performed to generate an FA metabolism score for each patient. Totally, 2981 and 182 lncRNAs from TCGA LIHC and GSE76427 datasets, respectively, were identified to be significantly related to FA metabolism scores (|Cor|> 0.25 and P < 0.05). However, only 70 lncRNAs from the two datasets overlapped and exhibited a co-positve or co-negative correlation with FA metabolism scores (Fig. 1A). Then, “ConsensusClusterPlus”, an unsupervised machine learning algorithm, was used to explore optimal molecular clusters. Two clusters were observed to be the most stable using the Cumulative distribution function Delta area curve (Fig. 1B, C). The overall survival rate of cluster 1 (C1) was significantly longer than that of cluster 2 (C2) (Fig. 1D). Additionally, the same phenomenon was observed in the GSE76427 cohort (Fig. 1E). Furthermore, the FA metabolism scores in C1 were significantly higher than those in C2 in the two cohorts (Fig. 1F, G).
Relationship between FA metabolism related-lncRNA signatures and clinical characteristics
The different clinical features between the two molecular subtypes, which were classified using the FA metabolism related-lncRNA signatures, were investigated. No significant difference was observed in age between the two clusters in the TCGA cohort (Fig. 2A). In line with the conventional speculations that males are more likely to develop liver cancer, the patients in C1 were observed to have a higher percentage of males (Fig. 2B). However, other clinical features, including T stage (tumor stage), N stage (lymph node stage), M stage (tumor metastasis stage), and Grade stage (pathological grade), were showed significant difference between the molecular C1 and C2 (Fig. 2C-G). The patients in C2 were associated with high-grade clinicopathological features, which partly explains the poor prognosis of patients in C2.
Genomic DNA damage and mutational characterization of high- and low-risk clusters
Malignancy and immune infiltration are associated with DNA damage assessment, including aneuploidy, homologous recombination deficiency (HRD), fraction of altered and number of segments. Aneuploidy is ubiquitously observed in various cancer, reflecting the degree of immune evasion and lower response to immunotherapy [22]. The aneuploidy scores of C1 were significantly lower than that of C2. Moreover, a strong negative correlation between aneuploidy scores and FA metabolism scores were observed in the TCGA cohort (Fig. 3A, B). Additionally, the degree of HRD was negatively correlated with the scores of FA metabolism. The patients in C1 had lesser HRD compare with those in C2. Previous studies reported that apoptosis-resistant tumors are correlated with a higher degree of HRDs and fraction altered but not with mutation rates [23]. Similarly, the level of fraction altered was significantly negatively correlated with the FA metabolism score, and the patients in C1 had lesser fraction altered than those in C2 (Fig. 3C-F). The patients in C1 had a lower number of segments than those in C2 (Fig. 3G). The number of segments were negatively correlated with the FA metabolism score without statistical significance (Fig. 3H). Moreover, the difference in the nonsilent mutation rate between C1 and C2 was not statistically significant (Fig. 3I, J). The mutation burden of C1 was smaller than that of C2, especially, in TP53, BAP1, SPEG, and ADRGL3 (Fig. 3K). The patients in C2 had a higher mutation rate of TP53 (42% vs. 24%), BAP1 (11% vs. 4%), SPEG (8% vs. 3%), and ADRGL3 (9% vs. 3%), the lower mutation rate of ATP10D, MYO1B, and RNF17 (0% vs. 4%) than those in C1 (Supplementary excel 4). Thus, C2 exhibited a higher degree of malignancy.
Distinct pathway patterns between FA metabolism-related lncRNA subgroups
To investigate the activated cellular signaling pathways underlying FA metabolism-associated lncRNAs, a different expressed genes analysis and GSEA between the two clusters were performed. As shown in Fig. 4A, compared with C1, 18 pathways were activated in C2, while 20 pathways were activated and 10 pathways were inhibited in the GSE76427 dataset. In total, 14 oncogenic pathways, including G2M checkpoint, E2F targets, MYC targets, mitotic spindle, allograft rejection, epithelial mesenchymal transition, WNT β-Catenin signaling, unfolded protein response, inflammatory response, mTORC1 signaling, apical junction, PI3K AKT mTOR signaling, IL2 STAT5 signaling, and TNFA signaling via NFKB overlapped in the two datasets (Fig. 4A). Moreover, Radar plots were used to describe consistently up-regulated pathways between C1 and C2 in the two cohorts. Furthermore, GSEA revealed that patients with subtype C2 showed up regulation in EMT, indicating that those lncRNAs could play an important role in connecting the hair surface of EMT, WNT, and tumor necrosis factor (Fig. 4B, C). Thus, the patients in C2 were speculated to possess high level EMT features and immune regulatory-related pathways, indicating that these lncRNAs could play important roles in FA metabolism, EMT, and the immune microenvironment.
Calculation of immune cell infiltration scores
The infiltration of immune cells, according to RNA-sequencing, was determined based on a previous study [24]. As shown in Fig. 5A, activated B cell, activated CD4 T cell, immature B cell, regulatory T cell, T follicular helper cell, type 2 T helper cell, activated dendritic cell, MDSC, and plasmacytoid dendritic cell in C2 were more significant than those in C1 in the TCGA LIHC and GSE76427 cohorts (Fig. 5). Additionally, CIBERSORT was used to assess the degree of immune cell infiltration in the TCGA LIHC and GSE76427 cohorts. The differences between C1 and C2 were mainly concentrated in the number of macrophages. The M0 macrophages were significantly upregulated in C2 in TCGA LIHC and GSE76427, while M1 and M2 macrophage were significantly downregulated in C2. Besides macrophages, mast cells resting were significantly decreased in C2 compared to C1. The proportion of activated NK cells was decreased in the TCGA LIHC cohort. However, activated NK cells in C2 were higher than that in C1 in the GSE76427 cohort (Fig. 5B).
Predictive value of the immunotherapy response of FA metabolism-related lncRNA signatures
Given the importance of immunotherapy in HCC treatment, the immune checkpoint genes in the HisgAtlas dataset were analysed [25], revealing that many immune checkpoint genes were significantly highly increased in C2 than that in C1. As shown in Fig. 6A, B and Supplementary Fig. 2, the high-risk subgroups of C2 had overexpression of several critical immune checkpoint genes, such as BTLA, CD160, CD27, CTLA4, HAVCR2(TIM-3), ICOS, IDO1, LAG3, TNFSF4, PDCD1 (PD-1), and ADORA2A, compared to those of C1. Additionally, tumor Immune dysfunction and exclusion (TIDE) software was used to assess the potential response to immunotherapy in the two clusters. The higher the TIDE predictive score, the higher the chance of immune escape, indicating that patients are less likely to benefit from immunotherapy (Fig. 6C, D). TIDE scores in C2 were higher than that of C1 in the TCGA cohort, suggesting that patients in C2 have a higher possibility of immune escape and lower benefits from immunotherapy (Fig. 6C). However, no significant difference was observed in the predictive value of immunotherapy response in the GSE76427 cohort (Fig. 6D).
Characterization of FA metabolism-related lncRNAs
Unlike protein-coding genes (PCG), most FA metabolism-related lncRNAs were negatively correlated with FA scores (Fig. 7A). Previous studies have shown that the function of lncRNA partly depends on its cellular localization. Usually, lncRNAs in the nucleus can activate or inhibit the transcriptional activities of target genes by directly binding to them, and regulating gene expression by participating in histone modification or recruitment of TFs [26]. lncRNAs in the cytoplasm usually interact with miRNAs as competitive endogenous RNA and mediate target gene expression [27]. “Relative concentration index (RCI)”, as defined by LncATLAS database [28] was used to describe the proportion of lncRNAs in the cytoplasm and nucleus. Additionally, 64.94% of the RCI was negative in the TCGA cohort, while 68.57% of the RCI was negative in the GSE76427 cohort, suggesting that approximately 70% of FA metabolism-related lncRNAs are localized in the nucleus (Fig. 7B).
Furthermore, the relationship between FA metabolism-related lncRNAs and dysregulated TFs was analyzed. To compare the TF activity between the two clusters, the TF activity score of each sample from TCGA and GSE76427 was calculated using the method developed by Garcia-Alonso [20]. A total of 123 and 70 significantly activated TFs were identified in C2 in the TCGA and GSE76427 cohorts, respectively. Additionally, analyzing correlation between nuclear lncRNAs and differential expressed TFs revealed a significantly negative correlation between a group of FA metabolism-related lncRNAs and a set of TFs (Fig. 7C). A total of 18 lncRNAs were positively correlated with FA scores, while the remaining 52 lncRNAs were negatively correlated with FA metabolism scores. Based on the correlation between lncRNAs and TF activity, as shown in Fig. 7D, most of the 12 lncRNAs in the signature were located in the nucleus. Compare with C1, nine TFs were significantly decreased and 17 TFs were significantly increased in C2 in the TCGA dataset, while four TFs were significantly down-regulated and three TFs were significantly up-regulated in C2 compared with C1 in the GSE76427 dataset (Fig. 7E). These 26 TFs were speculated to contribute to poor prognosis in C2. To validate this hypothesis, a cellular signaling pathway enrichment analysis of the 26 TFs was performed. The results showed that the downstream targets by the TFs were significantly enriched in a few critical cancer-related pathways, such as PI3K-Akt signaling pathway, cell cycle, FoxO signaling pathway, small cell lung cancer, p53 signaling pathway (Fig. 7F, Supplementary excel 5). Furthermore, 17 TFs expression were significantly increased and 9 TFs expression were significantly decreased in C2 compared with C1 in the TCGA dataset (Fig. 7G). Therefore, these FA metabolism-related lncRNAs and TFs could together promote HCC development.
Identification of key lncRNAs related to FA metabolism
To investigate the most crucial lncRNA in FA metabolism regulation, a first-order partial correlation analysis between FA metabolism scores and lncRNAs was performed. When the contribution of three key lncRNAs (SNHG1, LINC00261, and SNHG7) was removed, the correlation between FA metabolism scores and FA metabolism-related protein-coding genes (PCGs) decreased significantly (Fig. 8A), indicting their critical roles in linking FA metabolism and related signaling pathways. Then, using the FA metabolism-related PCGs and three key lncRNAs, consistent FA metabolism-related genes were identified, which were significantly enriched in FA degradation, drug metabolism-cytochrome P450, retinol metabolism, PPAR signaling pathway, and other signaling pathways (Fig. 8B). Base on the expression of the three key lncRNAs, the samples were divided into high- and low-risk groups using the best segmentation point. The overall survival of the high-risk groups was significantly shorter than the low-risk group in the two HCC datasets (Fig. 8C, D). Thus, the three lncRNAs occupy the core position in FA metabolism and are associated with the prognosis of patients with HCC.
The biological function of SNHG1 and SNHG7 in HCC cells
To further investigate the molecular function of the three FA metabolism-related lncRNAs, cox analysis was performed, which revealed that the expressions of SNHG1 and SNHG7 were significantly correlated with prognosis (P < 0.05), while LINC00261 was correlated with the prognosis without statistically significance (P = 0.092) (Fig. 9A). Thus, SNHG1 and SNHG7 were used for subsequent analysis. The results of qRT-PCR showed that the expression of SNHG1 was up-regulated in the hepG2 and Huh7 cell lines, while that of SNHG7 was upregulated in the most HCC cell lines (Fig. 9B). Additionally, lncRNA silencer was used to knock down SNHG1 expression, because SNHG1 is in the nucleus, while siRNA was used to inhibit SNHG7 expression, because SNHG7 is located in the cytoplasm. Smart silencer-mediated knockdown of SNHG1 and siRNA-mediated knockdown of SNHG7 were confirmed using qRT-PCR (Fig. 9C). FA Metabolism PCR Array was applied to explore the potential targets of SNHG1 and SNHG7, wherein, many FA metabolism-related genes were decreased on SNHG1 and SNHG7 knockdown. Among these potential downstream genes, five genes (ECHS1, MCEE, ACOT12, CPT1B, and BDH2) were up-regulated and 17 genes (ACSBG1, FABP6, ACADVL, ACSM3, ACOX2, CPT2, ECI2, ECHS1, DECR, SLC27A6, MUT, SLC27A4, ACAD10, FASN, ACSL4, ACADSB, and GK2) were decreased after SNHG1 knockdown in the HepG2 cell line (I Fold change I > 2) (Fig. 9D). Moreover, GK was up-regulated, and 14 genes (ACAA1, SLC27A4, ACSBG1, PRKAG3, FABP5, ACSL1, CPT1C, MUT, ACAD11, GPD2, ECHS1, GK2, CPT2, and CPT1B) were down-regulated after SNHG7 silencing (I Fold change I > 2) (Fig. 9E). Additionally, the signature negatively regulated fatty acid oxidation and fatty acid beta-oxidation using CoA oxidase among FA metabolism through GSEA analysis (Fig. 9F). ACSBG1 possesses long-chain Acyl-CoA synthetase, which activates fatty acids to their CoA derivatives, plays a central role in fatty acid metabolism. Previous literature revealed that ACSBG1 promote tumors to be more sensitive to ferroptosis. ACSBG1 was identified as a key pro-ferroptotic factor in this heat-induced ferroptosis process [29]. Thus, we chose ACSGB1 to perform western blot analysis. The results of western blot showed that the expression of ACSBG1 were significantly increased after SNHG1 inhibition, while the level of ACSBG1 is slightly increased after SNHG7 inhibition (Fig. 9G). Moreover, we performed the Oil Red O staining experiments and found that the amounts of intracellular lipid droplets are significantly decreased after knockdown of SNHG1 or SNHG7 in HepG2 and Huh 7 cell lines (Fig. 9H). Ferroptosis is recognized as a non-traditional form of programmed cell death, characterized by iron overload and lipid peroxidation. The cross-talk between ferroptosis and FA metabolism plays an important role in cancer progression. Hence, several key proteins were investigated by knocking down SNHG1 and SNHG7 in the HepG2 and Huh 7 cell lines. Moreover, GPX4, NRF2, NCOA4, and CD98 were downregulated whereas KEAP1 was upregulated after the knockdown of SNHG1 in the HepG2 and Huh7 cell lines (Fig. 9I). The levels of GPX4, NRF2, NCOA4, and CD98 were downregulated while that of KEAP1 was upregulated after the knockdown of SNHG7 in the HepG2 and Huh7 cell lines (Fig. 9J). Thus, SNHG1/7 not only play important role in FA metabolism, but also in ferroptosis.
The validation of SNHG1 and SNHG7 in HCC tissues
To validate our findings, we performed qRT-PCR analysis using 36 liver cancer tissue specimens from our center. According to the expression levels of SNHG1 or SNHG7, we divided 36 HCC tissues into high and low expression groups. We observed that the immune-infiltration related genes (HAVCR2, ICOS, LAG3, and PDCD1) were significantly increased in the SNHG1 high expression group (Fig. 10A), however, the level of HAVCR2, ICOS, LAG3, and PDCD1 were slightly increased without significance in the SNHG7 high expression group (Fig. 10B). On the other hand, the transcription factors (FOXO1, HNF1A, PPARA, PPARG, CEBPA, ONECUT1, CREB1, FOXA1, ZBTB7A, CEBPB, E2F4, and SMAD3) were remarkably overexpressed in the SNGH1 high expression group (Fig. 10C), the expression of FOXO1, HNF1A, PPARA, PPARG, CEBPA, ONECUT1, CREB1, FOXA1, ZBTB7A, and SMAD3 were significantly overexpressed in the SNHG7 high expression group (Fig. 10D).