Identification of Fatty acid metabolism-Related Genes in the Progression from Non-Alcoholic Fatty Liver to Nonalcoholic fatty liver disease

doi:10.21203/rs.3.rs-5030401/v1

Aim

To elucidate the regulatory mechanisms of fatty acid metabolism-related genes (FAMRGs) and the gene expression clustering in nonalcoholic fatty liver disease (NAFLD).

Methods

The NAFLD dataset was sourced from the Gene Expression Omnibus database. Differentially expressed genes (DEGs) were identified and their specific functions analyzed. Biomarkers were identified using machine learning and Weighted Gene Co-Expression Network Analysis. CIBERSORT evaluated immune cell infiltration and the relationship between biomarkers and immune cells. Fatty acid metabolism-related DEGs (FAMRDEGs) were identified, and consensus clustering differentiated NAFLD patients into two clusters. Clinical differences between subtypes were compared. Principal component analysis calculated cluster-specific gene scores, and single sample gene set enrichment analysis assessed the proportion of immune cells between clusters.

Results

A total of 2124 DEGs were identified, primarily associated with immune-related pathways. Among 44 FAMRDEGs, FMO1 was identified as a biomarker for NAFLD and validated using an independent dataset, qRT-PCR, and WB. Immune cell infiltration analysis suggested that NAFLD may be co-regulated by immune cells and FMO1. Clustering of NAFLD individuals based on the 44 FAMRDEGs revealed that genes in cluster A were predominantly related to immune pathways, while those in cluster B were associated with metabolic pathways. Disease severity was higher in cluster A, which also had a larger proportion of differing immune cells compared to cluster B.

Conclusion

FMO1 was identified as a biomarker for NAFLD. High expression of PPT1 and PTGS2 correlated with disease severity. The identification of NAFLD subgroups based on has enhanced our knowledge of NAFLD etiology.

non-alcoholic fatty liver disease

Fatty acid metabolism

machine learning

immune infiltration

molecular subgroups

Nonalcoholic fatty liver disease (NAFLD) is characterized by hepatic steatosis, marked by triglyceride accumulation within hepatocytes. It is closely associated with obesity, hypertension, diabetes, and dyslipidemia(1). NAFLD progression may lead to nonalcoholic steatohepatitis (NASH), posing a risk for cirrhosis and hepatocellular carcinoma (HCC), the leading cause of liver-related mortality in many Western nations and the second leading cause of liver transplantation. Approximately 25% of the global population is believed to have NAFLD(2). Besides its hepatic manifestations, NAFLD elevates the risk of various metabolic disorders, including diabetes, cardiovascular disease, and chronic kidney disease(3). Understanding the underlying mechanisms of NAFLD development is crucial for its prevention and treatment.

The pathogenesis of NAFLD remains elusive, with lipid metabolic disorders playing a pivotal role. Following the interplay of various pathogenic factors resulting in liver injury, dysregulation of lipid metabolism, endoplasmic reticulum stress, apoptosis, mitochondrial dysfunction, insulin resistance, and immune response contribute to hepatic steatosis and disease progression(2, 4, 5). Gene expression studies have investigated fatty acid metabolism in NAFLD patients(6, 7). However, the relationship between certain fatty acid metabolism-related genes and NAFLD remains unclear. Identification of key genes regulating NAFLD can provide insights into its pathogenesis and reveal potential diagnostic markers and therapeutic targets.

Gene chip technology offers a novel and efficient approach to investigating the molecular mechanisms of various diseases by discerning gene expression disparities between diseased and healthy tissues(8). Despite its potential, limited studies have utilized machine learning techniques such as least absolute shrinkage and selection operator (LASSO), support vector machine-recursive feature elimination (SVM-RFE), random forests (RF), and weighted gene co-expression network analysis (WGCNA) to identify individual-specific NAFLD biomarkers. Employing four machine learning algorithms enabled the accurate identification of biomarkers and predictive models of NAFLD, but most investigations have solely compared NAFLD cases with normal controls, overlooking nuanced case-by-case distinctions. In the study, we profiled NAFLD cases and juxtaposed the results with clinical characteristics which delineated varied expression patterns and severity within subgroups, suggesting an association between fatty acid metabolism-related genes and NAFLD.

Data collection and processing data

Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo/) provided GSE107231(9), GSE89632(10), GSE66676(11), GSE48452(12)and GSE49541(13). GSE107231 contains normal (n = 5) and NAFLD (n = 5) liver tissues; GSE89632 contains normal (n = 24) and NAFLD (n = 19) liver tissues; GSE66676 contains normal (n = 34) and NAFLD (n = 33) liver tissues. GSE49541 contains 15 mild (fibrosis stage 0–1) and 57 advanced (fibrosis stage 3–4) NAFLD liver tissues. GSE107231, GSE89632, and GSE66676 were based on GPL6244 platform. GSE49541 was based on GPL570 platform. Expression profiling of GSE48452 based on GPL11532 platform validates our results (including 12 normal samples and 17 NAFLD samples). Table 1 has dataset details. Figure S1 illustrates a work flow of this study.

Table 1

Characteristics of the GSE107231, GSE89632, GSE66676, GSE48452 and GSE49541 datasets.
The characteristic of the five datasets
Datasets	NAFLD		Normal	Platform
GSE107231	5		5	GPL6244
GSE89632	19		24	GPL6244
GSE66676	33		34	GPL6244
GSE48452	17		12	GPL11532
GSE49541	15(fibrosis stage 0–1)	57(fibrosis stage 3–4)	NA	GPL570
NAFLD, Nonalcoholic fatty liver disease

Research utilized annotation files to convert probes to gene symbols, excluding unmatched probes per dataset's annotation profile. For genes with multiple probes, the average expression level was used. All datasets were normalized using limma's normalize Between Arrays. ComBat was employed to adjust for batch effects across different platforms. The "sva" package was used to remove batch effects from GSE107231, GSE89632, and GSE666761, which were then merged due to their similar platforms. Principal component analysis (PCA) plots of training matrices before and after batch effect removal demonstrated the impact of inter-sample rectification.

Identification of Fatty acid metabolism-related Differentially Expressed Genes Between Normal and NAFLD samples

The Molecular Signature Database was queried for fatty acid metabolism-related genes (n = 309) (Table S1). The "limma" package identified DEGs, and the "ggplot2" package visualized their differential expression. Genes with an adjusted P < 0.05 were considered statistically significant. Overlapping DEGs with fatty acid metabolism-related genes (FAMRGs) yielded 44 fatty acid metabolism-related differentially expressed genes (FAMRDEGs).

Functional Enrichment Analysis of DEGs

The "clusterProfiler" package was used to analyze Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment of DEGs. Gene Set Enrichment Analysis (GSEA) was performed on the gene expression matrix using the clusterProfiler package with reference gene sets "c5.go.v7.4.symbols.gmt" and "c2.cp.kegg.v7.0.symbols.gmt." Significance was defined as q < 0.05 and adjusted P < 0.05.

Screening and Verifification of Diagnostic Markers

Four algorithms—RF, LASSO, SVM-RFE and WGCNA—were employed to identify novel and key biomarkers for NAFLD. The "randomForest" package implemented the RF method, while the "glmnet" package was used for LASSO, with the optimal lambda set to the minimum. The caret package's RFE function was used to select key genes, and the SVM classifier was developed using the "svmRadial" method from the e1071 package. WGCNA was performed using the "WGCNA" package. Overlapping genes from these four models were selected for further analysis. The GSE48452 dataset was used for validation, and the area under the curve (AUC) was computed to assess predictive efficacy.

Cell culture and treatment

HepG2 cells were purchased from xxxx with STR authentication and cultured in DMEM medium containing 10% fetal bovine serum (BI, USA) at 37°C in 5% CO₂. The cells were then treated with free fatty acids at a final concentration of 0.5 mM (oleic acid/palmitic acid ratio 2:1) for 24 hours. Cell steatosis was assessed by Oil red O and Lipi-Red staining, and total triglyceride content was determined by triglyceride enzyme method. Finally, qRT-PCR and western blot (WB) were used to detect the relative expression of FMO1.

The specific reagents and cell experimental methods are detailed in the Supplementary Information.

Evaluation and Correlation Analysis of Infifiltration-Related Immune Cells

The CIBERSORT website was used to filter 22 types of immune cell matrices, obtaining the immune cell infiltration matrix with P < 0.05. For PCA cluster analysis of this matrix, the "ggplot2" package was used. The "corrplot" package visualized the correlation among the 22 types of infiltrating immune cells with a correlation heatmap. The "ggstatsplot" and "ggplot2" packages analyzed and displayed the Spearman correlation between diagnostic markers and immune infiltrating cells.

Construction of subgroups based on consensus clustering

Using consensus clustering from the 'ConsensusClusterPlus' package, we classified NAFLD patients in GSE49541 into subgroups based on the expression matrix of FAMRDEGs from normal and NAFLD samples. Clustering was performed using the 'pam' method and 'Euclidean' distance, with a maximum of nine clusters. The consensus matrix determined the final number of clusters and consensus scores (> 0.8).

Identification of DEGs Between Two Subgroups and Gene Ontology Functional Enrichment Analysis

The "limma" package was employed to identify DEGs between two subgroups, with criteria of |logFC|>0.5 and adjusted P < 0.05. Subsequently, 102 DEGs were identified through differential analysis (Table S2). GO and KEGG functional enrichment analyses were conducted to elucidate the potential mechanisms of up- and down-regulated DEGs in the two subgroups, and the findings were visualized using an enrichment circle diagram.

Comparing the clinical characteristics of two subgroups

Clinical data, including disease severity, were obtained from the series matrix file downloaded from the GSE49541 dataset in the GEO. The pairwise proportion test was used to compare disease stage proportions between two subgroups. Disease stage proportions were analyzed as categorical variables, and their ratios were displayed in a histogram.

Subgroup-specific genes screening between subtypes

By comparing the two clusters, the subgroup-specific genes (PTGS2, GPD1, CYP2C19, IDH1, THRSP, AMACR, TP53INP2, ALAD, ACADL, ENO3, ERP29, PPT1, CYP4F2) were identified. We used boxplots and heatmaps to show the differences in 44 FAMRDEGs across subgroups.

Estimation of the Subgroup-specific genes Signature

To quantify gene patterns across different subgroups, we employed PCA to compute subgroup-specific gene scores for each sample. Initially, PCA was conducted to identify subgroups. Subsequently, subgroup-specific gene scores were calculated as follows: score for subgroup-specific genes = PC1i, where PC1 represents principal component 1 and i denotes the expression of subgroup-specific genes.

Estimation of Immune Cell Infiltration

Single sample GSEA (ssGSEA) was utilized to assess the abundance of immune cells in asthmatic samples. Initially, ssGSEA sequenced gene expression levels in the samples to establish their ranks. Subsequently, the expression levels of these genes were summed based on their presence in the input dataset. The quantity of immune cells in each sample was then calculated from the preceding analysis.

Statistics analysis

The moderate t-test was used to identify DEGs, whereas Fisher's exact test was used to evaluate GO and KEGG annotation enrichments, respectively. The Wilcoxon test was utilized to evaluate immune cells. All statistical analysis and image visualization was performed using the R programming language (version 4.1.1)

Data Preprocessing and DEG Screening

The expression matrices from GSE107231, GSE89632, and GSE66676 were merged, comprising 63 normal samples and 57 NAFLD samples. Subsequently, batch effects were removed, and PCA plots illustrated sample distributions pre- and post-normalization for batch effects (Figure S2). 2124 DEGs were identified from the normalized data, visualized through heat maps and volcano plots (Fig. 1A-B; Table S3). Furthermore, 44 FAMRGs overlapping with DEGs were presented via Venn diagrams and heat maps, consisting of 32 upregulated and 12 downregulated genes (Fig. 1C-D). The list of up- and down-regulated genes is provided in Table S4.

Functional Correlation Analysis

In the GO analysis (Fig. 2A; Table S2), enriched biological process (BP) terms included response to bacterial molecules, lipopolysaccharide response, regulation of autophagy, and chemokine production regulation. Cellular component terms were mainly associated with membrane rafts and nuclear chromosomes. Significant molecular function terms comprised DNA binding transcription activator activity, RNA polymerase II-specific activity, and DNA binding transcription activator activity. KEGG analysis revealed enrichment of DEGs in the TNF signaling pathway and the AGERAGE signaling pathway in diabetic complications (Fig. 2B; Table S5). GSEA results indicated enrichment in GO BP related to DNA repair, endoplasmic reticulum to Golgi vesicle-mediated transport, and protein targeting (Fig. 2C; Table S6), while nucleotide excision repair and peroxisome pathways were notably enriched in KEGG pathways (Fig. 2D; Table S6).

Screening and Verification of Diagnostic Markers

We selected genes with variances exceeding 20% (6222) for further analysis in the merged dataset. Setting the correlation coefficient threshold to 0.9, the soft thresholding power was adjusted to 7 (Figs. 3A-B). To consolidate similar modules in the cluster tree, the minimum gene count per module was set to 30, with a threshold of 0.3 (Fig. 3C). WGCNA analysis generated nine coexpression modules, some of which exhibited associations with NAFLD, as depicted by module-trait correlations (Fig. 3D). Notably, the pink module contained the most significant genes (221 in total). We identified 14 important biomarkers among DEGs by LASSO (Fig. 4A), while the SVM-RFE method recognized 22 genes from DEGs as important biomarkers (Fig. 4B). Additionally, the RF algorithm pinpointed four genes as crucial biomarkers (Fig. 4C-D). FMO1 emerged as an up-regulated gene in both the merged and validation datasets across all four algorithms (Fig. 5A-B; Table S7). ROC curves of FMO1 in the merged and validation datasets revealed their reliability as significant biomarkers, with AUCs of 0.744% and 0.833%, respectively (Fig. 5C-D), signifying high predictive accuracy. qRT-PCR analyses in NAFLD cell models validated the up-regulation of FMO1 (Fig. 5J), corroborating reproducibility and reliability was further supported by WB (Fig. 5H-I).

Results of Immune Cell Infiltration

Table S8 presents the results of CIBERSORT analysis on the merged data matrices. Statistically significant differences in immune cell infiltration were observed (Fig. 6A). Using the combined data matrix from GSE107231, GSE89632, and GSE666761 datasets, we investigated immune cell infiltration composition in NAFLD and normal samples (Fig. 6B). Results revealed a significantly higher proportion of naïve B lymphocytes, activated NK cells, eosinophils, and neutrophils in normal tissues. Conversely, the fraction of Macrophages M1, resting Dendritic cells (DCs), resting Mast cells, and active Mast cells was significantly higher in NAFLD tissues (Fig. 6C) (P < 0.05). Additionally, we assessed the relationship between 22 immune cells (Fig. 6D) and presented the associations among different immune cells in Table 2.

Table 2

Correlation results between immune cells.
Correlation	B cells naive	Mast cells activated	NK cells activated	Eosinophils	Neutrophils	Mast cells resting	Macrophages M1	Dendritic cells resting
B cells naive	1	0.17	0.16	−0.02	0.1	−0.32	−0.2	−0.07
Mast cells activated	0.17	1	0	0.3	0.38	−0.57	−0.02	−0.21
NK cells activated	0.16	0	1	−0.03	−0.08	−0.14	−0.03	−0.14
Eosinophils	−0.02	0.3	−0.03	1	0.27	−0.18	−0.14	−0.08
Neutrophils	0.1	0.38	−0.08	0.27	1	−0.41	−0.28	−0.31
Mast cells resting	−0.32	−0.57	−0.14	−0.18	−0.41	1	0.11	0.22
Macrophages M1	−0.2	−0.02	−0.03	−0.14	−0.28	0.11	1	0.2
Dendritic cells resting	−0.07	−0.21	−0.14	−0.08	−0.31	0.22	0.2	1

Correlation Analysis between FMO1 and Infiltrating Immune Cells

We investigated the relationship between immune infiltration results and FMO1. As shown in Fig. 8E, FMO1 was significantly associated with T cells gamma delta (P = 0.01), Macrophages M1 (P = 0.006), T cells CD8 (P = 0.002), DCs resting (P < 0.001), and Mast cells resting (P < 0.001), and negatively associated with Neutrophils (P < 0.001), Mast cells activated (P = 0.001), DCs activated (P = 0.001), and Monocyte (P = 0.002).

Two subtypes Identified by expression of the FAMRDEGs

The consensus clustering method was applied to identify two subgroups (cluster A and cluster B) based on the expression of 44 important FAMRDEGs (Fig. 7A-B). Cluster A had 33 instances, whereas cluster B comprised 39 instances (Table S9). The heat map and box plot were then applied to illustrate the differential expression levels of 44 important FAMRDEGs between the two clusters (Fig. 7C-D). PTGS2 and PPT1 were more highly expressed in cluster A than in cluster B, but GPD1, CYP2C19, IDH1, THRSP, AMACR, TP53INP2, ALAD, ACADL, ENO3, ERP29 and CYP4F2 demonstrated the contrary. The expression of 44 significant FAMRDEGs (Fig. 7E) and 13 subtype-specific FAMRDEGs (Fig. 7F) were able to properly distinguish the two clusters.

Functional enrichment analysis of differential genes between clusters

Between the two clusters, a total of 102 DEGs were uncovered (Fig. 8A; Table S10). GO and KEGG functional enrichment analysis show in Table S11. The majority of up-regulated genes in cluster A were enriched in the chemokinemediated signaling pathway, response to chemokine, collagen fibril organization, ECM − receptor interaction, Cytokine − cytokine receptor interaction, Chemokine signaling pathway and IL − 17 signaling pathway(Fig. 8B-C), while the majority of up-regulated genes in cluster B were enriched in the hormone metabolic process, alphaamino acid catabolic process, and cellular response to copper ion, Tryptophan metabolism, Arginine and proline metabolism and Biosynthesis of amino acids(Fig. 8D-E).

Comparing the clinical characteristics of two subgroups

The disease stage of NAFLD patients from the GSE49541 dataset was investigated to identify the clinical characteristics of two subgroups. In particular, cluster A revealed a more advanced illness stage than cluster B, indicating that cluster A individuals may be more severely damaged (Figure S3A). These results imply that the subgroups based on 44 FAMRDEGs may not only reflect the stages of NAFLD development, but also its fundamental biological characteristics. Comparing the clinical features of the two clusters revealed that cluster A had a more severe disease severity than cluster B.

Differential expression of 44 FAMRDEGs between clusters and correlation with immune cells

Using ssGSEA, we quantified immune cell quantities in NAFLD samples (Table S12) and compared differences between cluster A and B using a boxplot (Figure S3B). We explored the relationship between the 44 significant FAMRDEGs and immune cells, noting positive relationships with PTGS2, GPD1, THRSP, ALAD, ENO3, ERP29 and PPT1, particularly PTGS2 and PPT1 (Figure S3C). Conversely, CYP2C19, IDH1, AMACR, TP53INP2, ACADL and CYP4F2 exhibited negative associations with multiple immune cells. PTGS2 and PPT1, highly expressed in cluster A, positively correlated with immune cells. Conversely, lower-scoring groups displayed significant immune cell infiltration and negative correlations with specific immune cells. Additionally, the correlation between classification, clinical characteristics and immune cells indicates that immune cell infiltration differences contribute to the severity discrepancy between cluster A and B which illustrated in FigureS4-S6.

NAFLD is the most common chronic liver disease in the world and expected to be the leading cause of HCC worldwide by 2030(2). While liver biopsy is the diagnostic gold standard(14), it is invasive and risky highlighting the need for non-invasive biomarkers which could improve biomarker identification and therapy development. We analyzed publicly available human NAFLD microarray data, revealing FMO1 as a key gene in NAFLD pathophysiology.

We identified 2124 DEGs by comparing NAFLD and normal samples. GO pathways indicated regulation of lipopolysaccharide (LPS), autophagy, and chemokine synthesis. KEGG analysis revealed enrichment in the TNF-α and AGERAGE signaling pathways related to diabetes complications, while GSEA highlighted nucleotide excision repair and peroxisome pathways. Studies show that the TLR-4 signaling pathway, activated by intestinal endotoxin especially LPS, plays a role in chronic liver injury and the progression of NAFLD(15, 16). Oxidative stress from excessive triglycerides and free fatty acids in NAFLD inhibit autophagy while Autophagy-enhancing drugs can reduce hepatic steatosis and degrade Mallory-Denk bodies(17). Sara et al. found that severe histological profiles in pediatric NAFLD are associated with LPS-induced TNF-α overexpression in hepatic stellate cells (HSCs) which linked to elevated IL-1(18). In vitro study indicated that LPS-induced TNF-α is a key regulator of the pro-inflammatory response in HSCs, with p38MAPK inhibitors as potential therapies for hepatic inflammation in NASH. Zarei et al.(19)reported that PPARβ/δ, which regulates carbohydrate and lipid metabolism in the liver, may slow NAFLD progression when activated. The pathway enrichment analysis in this study supported these findings, confirming the reliability of the results.

This study is the first to conduct an in-depth analysis of fatty acid metabolism-related genes to elucidate their role in NAFLD. We employed various machine learning techniques to develop highly accurate predictive diagnostic models for screening significant NAFLD biomarkers. As results, FMO1 was identified as a key biomarker, validated through various ensemble approaches. The human FMO family includes FMO1 through FMO5, primarily responsible for oxidizing in drugs or foreign compounds with weak nucleophiles and catalyze the metabolism of drugs and dietary components using distinct catalytic processes and directly influence lipid metabolism(20–22). Mice lacking FMO1/FMO2/FMO4 are thinner with reduced triglycerides in white adipose tissue, and FMO1 protein is highly expressed in metabolic tissues(23). Shi et al.(24) found higher FMO1 levels in NASH, suggesting a crucial role in NAFLD pathogenesis. FMO1 and FMO3 primarily convert trimethylamine to trimethylamine-N-oxide (TMAO)(25), a compound linked to metabolic diseases like NAFLD and diabetes. Individuals with NAFLD have TMAO levels about three-and-a-half times higher than normal, and TMAO addition in mice leads to liver fat accumulation(26). TMAO also induces inflammation via NLRP3 inflammasome activation(27). Thus, FMO1 is crucial in regulating TMAO metabolism and related inflammation.

NAFLD is closely associated with various immune cells through CIBERSORT. In NAFLD tissues, there is a high expression of M1 macrophages, resting DCs, resting mast cells, and activated mast cells. Conversely, normal tissues exhibit high levels of naïve B cells, activated NK cells, eosinophils, and neutrophils. FMO1 expression is significantly elevated in NAFLD tissues and positively associated with gamma delta T cells, M1 macrophages, CD8 T cells, resting DCs, and resting mast cells. In contrast, neutrophils, activated mast cells, activated DCs, monocytes, eosinophils, and naïve B cells are negatively correlated with NAFLD. These findings suggest that FMO1 enhances immune cell infiltration and that both FMO1 and various inflammatory cells play a critical role in NAFLD progression, underscoring the importance of further research into the molecular mechanisms underlying NAFLD.

To further investigate the impact of key FAMRGs on NAFLD staging, we classified 72 NAFLD patients into clusters A and B based on FAMRDEG expression. Cluster A patients exhibited more advanced disease stages, indicating a more severe form of NAFLD compared to cluster B. These findings underscore the critical role of FAMRGs in NAFLD. By conducting an enrichment analysis of DEGs from the distinct clusters, we firstly delineated subtype differences based on FAMRGs and detailed the connections between specific pathways and NAFLD subtypes or clinical features. Functional enrichment analysis revealed that up-regulated genes in cluster A were primarily associated with inflammation and immune response, including the chemokine-mediated signaling pathway, response to chemokine, collagen fibril organization, ECM-receptor interaction, cytokine-cytokine receptor interaction, and the IL-17 signaling pathway.

Chemokines, which include the C, CC, CXC, and CX3C families, are small, highly conserved proteins that regulate pro-inflammatory and pro-oxidative processes in NAFLD pathogenesis(28). Previous research has shown that chemokines promote obesity by attracting pro-inflammatory monocytes into hypertrophic fat tissue, implicating them in energy metabolism, lipid metabolic diseases, and obesity(29). Cluster A was found to be primarily associated with chemokine-related pathways, further contributing to the inflammatory response in NAFLD. Thus, chemokine-related pathways may explain the differences in disease stage between clusters. ECM-receptor interaction, cytokine-cytokine receptor interaction, and the IL-17 signaling pathway are known to play significant roles in the immunological response and development of NAFLD. Our study identified a link between disease stage, collagen fibril organization, and severity in cluster A, suggesting that the collagen fibril organization pathway may have severe consequences in NAFLD patients. In contrast, up-regulated genes in cluster B are associated with hormone metabolism, alpha-amino acid catabolism, tryptophan metabolism, arginine and proline metabolism, and amino acid biosynthesis. Intestinal hormones influence food intake, body weight, and insulin resistance, thereby contributing to NAFLD(30, 31). These findings indicate that both hormonal and metabolic pathways significantly impact NAFLD, and the greater severity in cluster A is primarily due to immune cell infiltration and collagen fiber development.

Typing results showed NAFLD-related genes (PTGS2 and PPT1) were substantially expressed in cluster A. PTGS2 encodes Prostaglandin endoperoxide synthase 2 (PTGS2), also known as cyclooxygenase-2, is a key enzyme in prostaglandin biosynthesis and has been identified as a hallmark of ferroptosis, as have changes in NADPH levels and lipid peroxidation(32, 33). Studies have linked ferroptosis to various chronic liver diseases, including NASH and HCC. Unbalanced iron metabolism and reactive oxygen species-induced lipid peroxidation are proposed mechanisms of chronic liver injury in these conditions(34, 35). The role of PTGS2 in NAFLD is still under investigation. PPT1 encodes palmitoyl protein thioesterase 1 (PPT1), a thioesterase removes long-chain fatty acids from modified cysteine residues in proteins(36). Xu et al. reported elevated PPT1 expression in HCC tissues, and immunohistochemical analysis shown significantly higher PPT1 levels as well(37). Currently, there are no reports on PPT1 involvement in NAFLD, warranting further research.

We utilized ssGSEA to evaluate immune cell abundance in NAFLD samples and correlated 13 significant FAMRGs with immune cells at last. Cluster A exhibited expression of activated DCs, eosinophils, NKT cells, Th1s, and Th2s, indicating an inflammatory milieu. Liver DCs play a pivotal role in NAFLD by ingesting antigens and facilitating innate and adaptive immune responses(38). In NASH, DCs contribute to liver fibrosis and inflammation. Activation of cytotoxic T cells leads to increased release of proinflammatory cytokines such as IFN and TNF, and is crucial for macrophage recruitment and adipose tissue inflammation. Interaction between Tc cells and NKT cells enhances NASH development and hepatocellular carcinoma incidence(39). However, some studies suggest that NKT cells may have a limited role in NASH promotion(40, 41). NKT cells may either reduce liver inflammation and insulin resistance or potentially increase obesity and hepatic inflammation(42). Har et al. demonstrated that NAFLD progression is associated with increased eosinophilic type 2 liver inflammation, contrasting with type 1 inflammation and decreased eosinophils seen in expanding adipose tissue(43). These findings underscore the strong association of NK cells, eosinophils, T cells, and DCs with NAFLD progression, supporting the severity of cluster A. Furthermore, the favorable correlation between PTGS2, PPT1, and immune cells suggests their involvement in immune cell infiltration, highlighting the significance of FAMRGs in NAFLD development.

FMO1 was identified as a biomarker for NAFLD and associated with multiple immune cell types, suggesting their potential influence on NAFLD development. Further investigation into the role of these immune cells may unveil new immunotherapeutic targets. Elevated expression of PPT1 and PTGS2 in NAFLD correlates with disease severity, indicating their potential as promising therapeutic targets. The identification of NAFLD subgroups based on FAMGs has advanced our understanding of NAFLD etiology and facilitated the development of diagnostic, classification, and prognostic approaches for the disease.

FAMRGs, fatty acid metabolism-related genes; NAFLD, nonalcoholic fatty liver disease; DEGs, differentially expressed genes; FAMRDEGs, fatty acid metabolism-related DEGs; NASH, nonalcoholic steatohepatitis; HCC, hepatocellular carcinoma; LASSO, least absolute shrinkage and selection operator; SVM-RFE, support vector machine-recursive feature elimination; RF, random forests; WGCNA, weighted gene co-expression network analysis; GEO, Gene Expression Omnibus; PCA, Principal component analysis; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; GSEA, Gene Set Enrichment Analysis; AUC, area under the curve; WB, western blot; ssGSEA, Single sample GSEA; BP, biological process; DCs, dendritic cells; LPS, lipopolysaccharide; HSCs, hepatic stellate cells; TMAO, trimethylamine-N-oxide; PTGS2, Prostaglandin endoperoxide synthase 2; PPT1, palmitoyl protein thioesterase 1;

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

The datasets ( GSE107231, GSE89632, GSE66676, GSE48452 and GSE49541) for this study can be found in the GEO (https://www.ncbi.nlm.nih.gov/geo/).

Competing interests

The authors declare that they have no competing interests.

Funding

This work was supported by National Natural Science Foundation of China [No. 82302332 to HF] and The Department of Science and Technology of Yunnan Province Kunming Medical University joint project on basic research [NO.202401AY070001-122 to XL].

Author contributions

All authors made significant contributions to the conception and design of the study. XX and TX contributed equally to this work. XX, YZ, XL, WS, MG, XH were responsible for material preparation, data collection and analysis. The initial draft of the manuscript was written by XX and TX. HF and GW were responsible for formal analysis, investigation, conceptualization, writing review and editing. All authors contributed to manuscript revision, read, and approved the submitted version.

Acknowledgments

GEO facilitated in the completion of this work. We would like to express our gratitude to the GEO network for freely sharing huge amounts of data.

Kanwal F, Kramer JR, Li L, Dai J, Natarajan Y, Yu X, et al. Effect of Metabolic Traits on the Risk of Cirrhosis and Hepatocellular Cancer in Nonalcoholic Fatty Liver Disease. Hepatology. 2020;71(3):808–19.
Vitale A, Svegliati-Baroni G, Ortolani A, Cucco M, Dalla Riva GV, Giannini EG, et al. Epidemiological trends and trajectories of MAFLD-associated hepatocellular carcinoma 2002–2033: the ITA.LI.CA database. Gut. 2023;72(1):141–52.
Marusic M, Paic M, Knobloch M, Liberati Prso AM. NAFLD, Insulin Resistance, and Diabetes Mellitus Type 2. Can J Gastroenterol Hepatol. 2021;2021:6613827.
Seval GC, Kabacam G, Yakut M, Seven G, Savas B, Elhan A, et al. The natural course of non-alcoholic fatty liver disease. Hepatol Forum. 2020;1(1):20–4.
Martel C, Esposti DD, Bouchet A, Brenner C, Lemoine A. Non-alcoholic steatohepatitis: new insights from OMICS studies. Curr Pharm Biotechnol. 2012;13(5):726–35.
Kohjima M, Enjoji M, Higuchi N, Kato M, Kotoh K, Yoshimoto T, et al. Re-evaluation of fatty acid metabolism-related gene expression in nonalcoholic fatty liver disease. Int J Mol Med. 2007;20(3):351–8.
Kohjima M, Higuchi N, Kato M, Kotoh K, Yoshimoto T, Fujino T, et al. SREBP-1c, regulated by the insulin and AMPK signaling pathways, plays a role in nonalcoholic fatty liver disease. Int J Mol Med. 2008;21(4):507–11.
Wang T, Zheng X, Li R, Liu X, Wu J, Zhong X, et al. Integrated bioinformatic analysis reveals YWHAB as a novel diagnostic biomarker for idiopathic pulmonary arterial hypertension. J Cell Physiol. 2019;234(5):6449–62.
Guo J, Fang W, Sun L, Lu Y, Dou L, Huang X, et al. Ultraconserved element uc.372 drives hepatic lipid accumulation by suppressing miR-195/miR4668 maturation. Nat Commun. 2018;9(1):612.
Arendt BM, Comelli EM, Ma DW, Lou W, Teterina A, Kim T, et al. Altered hepatic gene expression in nonalcoholic fatty liver disease is associated with lower hepatic n-3 and n-6 polyunsaturated fatty acids. Hepatology. 2015;61(5):1565–78.
Xanthakos SA, Jenkins TM, Kleiner DE, Boyce TW, Mourya R, Karns R, et al. High Prevalence of Nonalcoholic Fatty Liver Disease in Adolescents Undergoing Bariatric Surgery. Gastroenterology. 2015;149(3):623–34 e8.
Ahrens M, Ammerpohl O, von Schonfels W, Kolarova J, Bens S, Itzel T, et al. DNA methylation analysis in nonalcoholic fatty liver disease suggests distinct disease-specific and remodeling signatures after bariatric surgery. Cell Metab. 2013;18(2):296–302.
Murphy SK, Yang H, Moylan CA, Pang H, Dellinger A, Abdelmalek MF, et al. Relationship between methylome and transcriptome in patients with nonalcoholic fatty liver disease. Gastroenterology. 2013;145(5):1076–87.
Meng Q, Li X, Xiong X. Identification of Hub Genes Associated With Non-alcoholic Steatohepatitis Using Integrated Bioinformatics Analysis. Front Genet. 2022;13:872518.
Soares JB, Pimentel-Nunes P, Roncon-Albuquerque R, Leite-Moreira A. The role of lipopolysaccharide/toll-like receptor 4 signaling in chronic liver diseases. Hepatol Int. 2010;4(4):659–72.
Miura K, Seki E, Ohnishi H, Brenner DA. Role of toll-like receptors and their downstream molecules in the development of nonalcoholic Fatty liver disease. Gastroenterol Res Pract. 2010;2010:362847.
Sun L, Zhang S, Yu C, Pan Z, Liu Y, Zhao J, et al. Hydrogen sulfide reduces serum triglyceride by activating liver autophagy via the AMPK-mTOR pathway. Am J Physiol Endocrinol Metab. 2015;309(11):E925-35.
Ceccarelli S, Panera N, Mina M, Gnani D, De Stefanis C, Crudele A, et al. LPS-induced TNF-alpha factor mediates pro-inflammatory and pro-fibrogenic pattern in non-alcoholic fatty liver disease. Oncotarget. 2015;6(39):41434–52.
Zarei M, Aguilar-Recarte D, Palomer X, Vazquez-Carrera M. Revealing the role of peroxisome proliferator-activated receptor beta/delta in nonalcoholic fatty liver disease. Metabolism. 2021;114:154342.
Krueger SK, Williams DE. Mammalian flavin-containing monooxygenases: structure/function, genetic polymorphisms and role in drug metabolism. Pharmacol Ther. 2005;106(3):357–87.
Varshavi D, Scott FH, Varshavi D, Veeravalli S, Phillips IR, Veselkov K, et al. Metabolic Biomarkers of Ageing in C57BL/6J Wild-Type and Flavin-Containing Monooxygenase 5 (FMO5)-Knockout Mice. Front Mol Biosci. 2018;5:28.
Zhang T, Yang P, Wei J, Li W, Zhong J, Chen H, et al. Overexpression of flavin-containing monooxygenase 5 predicts poor prognosis in patients with colorectal cancer. Oncol Lett. 2018;15(3):3923–7.
Veeravalli S, Omar BA, Houseman L, Hancock M, Gonzalez Malagon SG, Scott F, et al. The phenotype of a flavin-containing monooyxgenase knockout mouse implicates the drug-metabolizing enzyme FMO1 as a novel regulator of energy balance. Biochem Pharmacol. 2014;90(1):88–95.
Shi C, Pei M, Wang Y, Chen Q, Cao P, Zhang L, et al. Changes of flavin-containing monooxygenases and trimethylamine-N-oxide may be involved in the promotion of non-alcoholic fatty liver disease by intestinal microbiota metabolite trimethylamine. Biochem Biophys Res Commun. 2022;594:1–7.
Bennett BJ, de Aguiar Vallim TQ, Wang Z, Shih DM, Meng Y, Gregory J, et al. Trimethylamine-N-oxide, a metabolite associated with atherosclerosis, exhibits complex genetic and dietary regulation. Cell Metab. 2013;17(1):49–60.
Tan X, Liu Y, Long J, Chen S, Liao G, Wu S, et al. Trimethylamine N-Oxide Aggravates Liver Steatosis through Modulation of Bile Acid Metabolism and Inhibition of Farnesoid X Receptor Signaling in Nonalcoholic Fatty Liver Disease. Mol Nutr Food Res. 2019;63(17):e1900257.
Sun X, Jiao X, Ma Y, Liu Y, Zhang L, He Y, et al. Trimethylamine N-oxide induces inflammation and endothelial dysfunction in human umbilical vein endothelial cells via activating ROS-TXNIP-NLRP3 inflammasome. Biochem Biophys Res Commun. 2016;481(1–2):63–70.
Koutoukidis DA, Astbury NM, Tudor KE, Morris E, Henry JA, Noreik M, et al. Association of Weight Loss Interventions With Changes in Biomarkers of Nonalcoholic Fatty Liver Disease: A Systematic Review and Meta-analysis. JAMA Intern Med. 2019;179(9):1262–71.
Xu L, Kitade H, Ni Y, Ota T. Roles of Chemokines and Chemokine Receptors in Obesity-Associated Insulin Resistance and Nonalcoholic Fatty Liver Disease. Biomolecules. 2015;5(3):1563–79.
Gribble FM, Reimann F. Function and mechanisms of enteroendocrine cells and gut hormones in metabolism. Nat Rev Endocrinol. 2019;15(4):226–37.
Armstrong MJ, Gaunt P, Aithal GP, Barton D, Hull D, Parker R, et al. Liraglutide safety and efficacy in patients with non-alcoholic steatohepatitis (LEAN): a multicentre, double-blind, randomised, placebo-controlled phase 2 study. Lancet. 2016;387(10019):679–90.
Shimada K, Hayano M, Pagano NC, Stockwell BR. Cell-Line Selectivity Improves the Predictive Power of Pharmacogenomic Analyses and Helps Identify NADPH as Biomarker for Ferroptosis Sensitivity. Cell Chem Biol. 2016;23(2):225–35.
Wang H, An P, Xie E, Wu Q, Fang X, Gao H, et al. Characterization of ferroptosis in murine models of hemochromatosis. Hepatology. 2017;66(2):449–65.
Mehta KJ, Farnaud SJ, Sharp PA. Iron and liver fibrosis: Mechanistic and clinical aspects. World J Gastroenterol. 2019;25(5):521–38.
Li S, Tan HY, Wang N, Zhang ZJ, Lao L, Wong CW, et al. The Role of Oxidative Stress and Antioxidants in Liver Diseases. Int J Mol Sci. 2015;16(11):26087–124.
Segal-Salto M, Sapir T, Reiner O. Reversible Cysteine Acylation Regulates the Activity of Human Palmitoyl-Protein Thioesterase 1 (PPT1). PLoS One. 2016;11(1):e0146466.
Xu J, Su Z, Cheng X, Hu S, Wang W, Zou T, et al. High PPT1 expression predicts poor clinical outcome and PPT1 inhibitor DC661 enhances sorafenib sensitivity in hepatocellular carcinoma. Cancer Cell Int. 2022;22(1):115.
Schuster S, Cabrera D, Arrese M, Feldstein AE. Triggering and resolution of inflammation in NASH. Nat Rev Gastroenterol Hepatol. 2018;15(6):349–64.
Ghazarian M, Revelo XS, Nohr MK, Luck H, Zeng K, Lei H, et al. Type I Interferon Responses Drive Intrahepatic T cells to Promote Metabolic Syndrome. Sci Immunol. 2017;2(10).
Heymann F, Tacke F. Immunology in the liver–from homeostasis to disease. Nat Rev Gastroenterol Hepatol. 2016;13(2):88–110.
Bhattacharjee J, Kirby M, Softic S, Miles L, Salazar-Gonzalez RM, Shivakumar P, et al. Hepatic Natural Killer T-cell and CD8 + T-cell Signatures in Mice with Nonalcoholic Steatohepatitis. Hepatol Commun. 2017;1(4):299–310.
Hams E, Locksley RM, McKenzie AN, Fallon PG. Cutting edge: IL-25 elicits innate lymphoid type 2 and type II NKT cells that regulate obesity in mice. J Immunol. 2013;191(11):5349–53.
Hart KM, Fabre T, Sciurba JC, Gieseck RL, 3rd, Borthwick LA, Vannella KM, et al. Type 2 immunity is protective in metabolic disease but exacerbates NAFLD collaboratively with TGF-beta. Sci Transl Med. 2017;9(396).

No competing interests reported.

Identification of Fatty acid metabolism-Related Genes in the Progression from Non-Alcoholic Fatty Liver to Nonalcoholic fatty liver disease

Status:

Version 1

Abstract

Aim

Methods

Results

Conclusion

Figures

Introduction

Materials and methods

Data collection and processing data

Identification of Fatty acid metabolism-related Differentially Expressed Genes Between Normal and NAFLD samples

Functional Enrichment Analysis of DEGs

Screening and Verifification of Diagnostic Markers

Cell culture and treatment

Evaluation and Correlation Analysis of Infifiltration-Related Immune Cells

Construction of subgroups based on consensus clustering

Identification of DEGs Between Two Subgroups and Gene Ontology Functional Enrichment Analysis

Comparing the clinical characteristics of two subgroups

Subgroup-specific genes screening between subtypes

Estimation of the Subgroup-specific genes Signature

Estimation of Immune Cell Infiltration

Statistics analysis

Results

Data Preprocessing and DEG Screening

Functional Correlation Analysis

Screening and Verification of Diagnostic Markers

Results of Immune Cell Infiltration

Two subtypes Identified by expression of the FAMRDEGs

Functional enrichment analysis of differential genes between clusters

Comparing the clinical characteristics of two subgroups

Differential expression of 44 FAMRDEGs between clusters and correlation with immune cells

Discussion

Conclusions

Abbreviations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1