The workflow is shown in Figure 1 . Original bulk RNA sequencing (RNA-seq) data was collected from the gene expression profile GSE126848, based on the GPL18573 Illumina NextSeq 500 (Homo sapiens) downloaded from the Gene Expression Omnibus (GEO) database, a public depository database of gene expression data (13). Data of single-cell RNA sequencing (scRNA-seq) from GSE136103 were obtained from the GEO database (14). The GSE136103 dataset included single-cell transcriptome data from 11 healthy patient liver tissue samples, 9 cirrhotic patient liver tissue samples, 4 cirrhotic patient PBMC samples, and 2 mouse liver samples. Only 4 human NAFLD-associated cirrhotic samples were included, so all of them were included in our study as the NAFLD group and 4 human healthy liver tissue samples were randomly selected as the control group. Additional information regarding the aforementioned datasets is provided in Supplementary Table S1 .

We divided genes into different modules based on their similar expression in 57 liver samples using the WGCNA method (15). After evaluating the significance of genes in the module and analyzing the correlation between modules and subtypes, a module strongly associated with NAFLD was discovered. The genes within this module were then utilized for further investigation.

The RNA-seq data was analyzed for differentially expressed genes (DEGs) using the DESeq2 R package (version 1.40.2) (16). DEGs were classified as significant if their absolute log2 fold change (log2FC) was greater than 0.5 and their p-value was less than 0.05. Ensemble ID with no corresponding gene symbols were removed, and genes with more than one Ensemble ID were averaged.

With the “clusterProfiler” package (version 4.8.1), we conducted GSEA and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses (17, 18). The “c2.cp.v2023.2.Hs.symbols.gmt” was used as the annotation gene set (19). Significant terms were only considered if their p-value was less than 0.05. The outcomes of the enrichment analyses were graphically represented through the use of the “clusterProfiler” package (20), “enrichplot” package (version 1.20.0), and “GseaVis” package (version 0.0.9) in R 4.3.1 software.

In our research, we utilized three algorithms, CIBERSORT, xCell, and ssGSEA, to quantify immune cell infiltration in NAFLD liver tissues (21–23). We assessed the relative proportions of immune cells in each sample. Additionally, we employed Spearman correlation analysis to assess the correlation between the relative proportions of different cells in CIBERSORT and the mRNA levels of key DE-DRGs. Finally, we visualized the results of immune infiltration analysis using the R package “ggplot2” (version 3.4.4).

Liu’s study provided the information for 23 DRGs (SLC7A11, CD2AP, PDLIM1, ACTN4, MYH10, IQGAP1, FLNA, TLN1, MYL6, ACTB, CAPZB, GYS1, NDUFA11, NCKAP1, NDUFS1, RPN1, NUBPL, SLC3A2, INF2, MYH9, FLNB, DSTN, LRPPRC) (7). The “RCircos” package (version 1.2.2) was utilized to map the positioning of DRGs on chromosomes. Subsequently, the key DE-DRGs were identified from the DEGs, the module of WGCNA, and the DRGs using an interactive Venn diagram. Consensus clustering analysis was conducted using the R package “ConsensusClusterPlus” (version 1.64.0) to categorize patients into distinct clusters according to the expression levels of the key DE-DRGs. The range of the cluster category’s k value was between 2 and 9. Differences in expression of key DE-DRGs between clusters as reflected by heat map.

To identify reliable cell subpopulations from scRNA-seq data, we employed the “Seurat” package (version 4.3.0.1) for data processing (24). Using the Seurat R package, we applied normalization, clustering, dimensionality reduction, and visualization techniques to analyze the scRNA-seq data. Each gene must be expressed in a minimum of three cells. The PercentageFeatureSet function was used to calculate the proportions of mitochondria and ribosomal RNA, making sure that each cell had over 200 genes and less than 6,000 genes with less than 25% mitochondrial content and less than 50% ribosome content per cell. To eliminate the influence of confounding factors, the CellCycleSoring function was utilized. To correct batch effects, we employed the RunHarmony function from the harmony R package using the PCA reduction method to integrate samples (25). We used the FindAllMarkers function to identify DEGs for determining the cellular identity. We utilized a three-step workflow including automatic cell annotation, manual cell annotation and verification. The SingleR package (version 2.2.0) is used for automatic annotation. Marker genes used to define different cell types were selected based on the literature, and these cells were identified by their typical characteristic gene profiles (26–28). Finally, we validated the annotation results with the online tool CellMaker (29). During the examination of each cell type, any group of cells that exhibit the presence of multiple markers from various cell types will be eliminated. The condition-specific responses of cell subpopulations measured from two groups of patients were investigated using the “muscat” package (version 1.15.1) through differential state analysis (30). To demonstrate the functional characteristics of mononuclear phagocytes (MPs), we gathered gene collections from various sources and computed functional scores for each individual cell (28, 31–34). The gene list scores were calculated using the AddModuleScore function in Seurat. Online Supplementary Table S2 contained detailed information for each score. Finally, certain R packages were used for the visualization of single-cell data results, such as the “plot1cell”, “Nebulosa”, “scRNAtoolVis” and “ggplot2” packages (35–37).

This study involved 16 patients who underwent surgical procedures, with six diagnosed with NAFL, four with NASH, and six serving as controls ( Supplementary Table S3 ). Two experienced hepatopathologists, who were unaware of clinical data, centrally evaluated liver specimens from all of these patients. Liver samples were assessed for histologic characteristics using the NAFLD activity score (NAS) (38). This research received approval from the Ethics Committee of the Second Xiangya Hospital of Central South University (No. 2019–050). Participant consent was obtained in writing prior to participating in the study.

In the initial step, liver tissue samples were fixed using 4% paraformaldehyde. After fixation, the tissues were dehydrated in low to high concentrations of ethanol. Following dehydration, the tissues were transparent in xylene and then embedded in paraffin wax for fixation. Afterward, the tissue was sliced to a thickness of 3–5 um, baked for about 1 h, immersed in xylene to deparaffinize, and then hydrated in ethanol from high to low concentration. Once drained, the slices were placed in a vat of hematoxylin staining solution for a duration of 5 minutes, followed by a thorough rinse with running water to restore their blue color. The sections were then immersed in an eosin vat for 5 min and rinsed again under running water. Following a successful staining process, the sections underwent dehydration using varying concentrations of ethanol, followed by transparency achieved with xylene. Subsequently, the sections were sealed and finally examined and captured using a microscope.

The RNAmini kit (Qiagen, Germany) was used to extract total RNA from every liver tissue sample. The quality of RNA was assessed using gel electrophoresis and Qubit (Thermo, Waltham, MA, USA). The TruSeq RNA sample preparation kit (Illumina, San Diego, CA, USA) was used to construct strand-specific libraries. The Illumina Novaseq 6000 instrument was used for sequencing, which was performed by Genergy Biotechnology Co. Ltd. (Shanghai, China). Skewer was used to process the raw data and FastQC v0.11.2 was utilized to verify the quality of the data. The length of the read was 2×150 base pairs. STAR (version 2.5.3a) was employed to align the clean reads to the Human genome hg38. StringTie software (version 1.3.1c) was utilized to count the raw sequence counts of established genes in all samples.

Statistical analysis was performed using R version 4.3.1 and GraphPad Prism 9.3.0 software. The Wilcoxon test was utilized to complete the comparative analysis among the groups. Correlation analysis was conducted using Spearman analysis. The gene expression levels in this study were reported as log2(TPM+1) unless otherwise specified. The statistical difference between the sets was significant with a p-value less than 0.05.

In dataset GSE126848, we analyzed 57 liver tissue samples by applying median absolute deviation to screen the top 5000 genes for WGCNA analysis. A clustering tree for the dataset was established and no outliers were found ( Figure 2A ). A dendrogram was formed by clustering samples with clinical traits ( Figure 2B ). Afterwards, we assessed the fitting index for scale-free networks and the average connectivity for different soft threshold powers, using the scale-free R² as a basis ( Figure 2C ). Ultimately, we identified 13 modules through hierarchical clustering ( Figure 2D ). We utilized a correlation heatmap to investigate the correlation between each module and NAFLD. Our analysis revealed that the yellow module exhibited the strongest connection with NAFLD, displaying a correlation value of 0.85 and a p-value of 7E-17 ( Figure 2E ). Figure 2F showed that the yellow module genes exhibited a strong correlation not only with the yellow module itself but also with NAFLD. A heatmap was constructed to display the network of gene co-expression ( Figure 2G ). Figure 2H displayed the connections between the modules and the clinical trait weight. Therefore, we focused on the 564 genes in the yellow module in our follow-up study ( Supplementary Table S4 ).

We examined the DEGs in the GSE126848 dataset to investigate the underlying gene expression pattern and possible correlation between control and NAFLD. By removing the duplicate and irrelevant genes, we discovered a combined count of 2702 up-regulated DEGs and 2075 down-regulated DEGs ( Supplementary Table S5 ). Figure 3A displayed volcano maps of the DEGs. Subsequently, we identified 312 NAFLD-specific DEGs in DEGs and WGCNA, of which 293 were up-regulated and 19 were down-regulated ( Figure 3B ). We used a heat map to show the expression of NAFLD-specific DEGs in the two groups of samples ( Figure 3C ). In addition, KEGG pathway enrichment analysis was used in R to determine the pathways involved in NAFLD-specific DEGs ( Figure 3D ; Supplementary Table S6 ). The up-regulated NAFLD-specific DEGs were implicated in Ribosome, Oxidative phosphorylation, Non-alcoholic fatty liver disease, Diabetic cardiomyopathy, Motor proteins, Phagosome, and PPAR signaling pathways. Conversely, the down-regulated NAFLD-specific DEGs participated in signaling pathways related to Pentose and glucuronate interconversions, Arachidonic acid metabolism, Linoleic acid metabolism, and Retinol metabolism. Interestingly, we found that 17 of the screened NAFLD-specific DEGs were significantly enriched in the Non-alcoholic fatty liver disease signaling pathway ( Figure 3E ). This confirms the effectiveness of our screening method and also identifies novel genes related to NAFLD. Overall, the NAFLD-specific DEGs identified in this study might be involved in the occurrence and progression of NAFLD.

To investigate the immune microenvironment in patients with NAFLD, we applied to three algorithms, CIBERSORT, xCell, and ssGSEA, to revealed relative immune cell composition ( Figure 4 ; Supplementary Figure S1 ). The findings indicated that there is a variation in the immune microenvironment between patients with NAFLD and those without the condition ( Figure 4A ; Supplementary Table S7 ). In addition, there was a close interplay between immune cells, with a significant negative correlation between M1 macrophages and M2 macrophages, and a strong positive correlation between activated dendritic cells and follicular helper T cells ( Figure 4B ). In patients with NAFLD, the proportion of plasma cells, CD8 T cells, CD4 memory activated T cells, gamma delta T cells, M0 macrophages, M1 macrophages, and resting dendritic cells was found to be higher than those in control samples. However, the percentage of activated NK cells, M2 macrophages, resting mast cells, and neutrophils was lower in NAFLD samples ( Figures 4C, D ; Supplementary Figure S1 ). All of these results suggested that alterations in the immune microenvironment play a crucial role in the pathogenesis of NAFLD.

The location of 23 DRGs on chromosomes were demonstrated in Figure 5A . To explore the role of DRGs in NAFLD, we examined the expression patterns of 23 DRGs in the control and NAFLD groups. Figure 5B showed that 7 DRGs were differentially expressed in the two groups. After intersecting with the yellow module genes obtained from WGCNA, two key DE-DRGs associated with NAFLD, MYL6 and DSTN, were identified. We then grouped the samples using a consistent clustering algorithm based on the key DE-DRGs. Two disulfidptosis clusters were obtained by cluster analysis ( Figures 5C, D ). Compared to cluster 2, the key DE-DRGs were up-regulated in cluster 1 ( Figure 5E ). To clarify the characteristics of the two disulfidptosis clusters in the immune microenvironment, an immune infiltration analysis was performed ( Figure 5F ). It is worth noting that the analysis of correlation depicted in Figure 5G demonstrated a significant positive correlation between the main DE-DRGs and regulatory T cells, CD8 T cells, M0 macrophages, M1 macrophages, and resting dendritic cells. On the other hand, there was a significant negative correlation between the DE-DRGs and M2 macrophages, as well as resting mast cells. Subsequently, GSEA was applied to analyze the biological behavior of the two clusters. The GSEA enrichment analysis showed that the electron transport chain oxphos system in mitochondria, oxidative stress induced senescence, metabolic disorders of biological oxidation enzymes, triglyceride metabolism, diseases of programmed cell death and proinflammatory and profibrotic mediators were upregulated in cluster 1 ( Figure 5H ), all of which have been shown to be associated with the pathogenesis of NAFLD (39, 40). In conclusion, these results strongly indicated that the key DE-DRGs might be involved in the development and progression of NAFLD and were related to immune modulation.

We further explored the immunomodulatory role of the key DE-DRGs in NAFLD patients through scRNA-seq data. After rigorous quality control screening, a total of 20896 genes and 23929 high-quality cells were included in the subsequent analysis ( Figure 6A ). The FindAllMarkers function and the Wilcoxon rank sum test were used to identify specific gene signatures for each cluster, as shown in Figure 6B . We finally identified 10 major cell lineages based on marker genes from previous studies and SingleR annotation, including T/NK cells (13083 cells, 54.7%, marked with CD3D, CD3E, CD7, and KLRF1), B cells (930 cells, 3.9%, marked with CD79A, CD79B, and MS4A1), plasma cells (286 cells, 1.2%, marked with IGHA1, IGHG1, and MZB1), MPs (3496 cells, 14.6%, marked with CD68, CD163, CD14, and C1QA), dendritic cells (DCs) (479 cells, 2.0%, marked with CD1C and CLEC10A), mast cells (57 cells, 0.2%, marked with TPSAB1, TPSB2, and CPA3), endothelial cells (ECs) (3475 cells, 14.5%, marked with CLEC4G, TSPAN7, and VWF), hepatic stellate cells (HSCs) (483 cells, 2.0%, marked with ACTA2 and COL1A1), hepatocytes (122 cells, 0.5%, marked with ALB, APOA1, and APOC3), and cholangiocytes (1518 cells, 6.3%, marked with ANXA4, FXYD2, and KRT19) (26–28). These ten populations contained cells from both controls and NAFLD patients’ liver tissues ( Figure 6C ). The cells were grouped into clusters using an unsupervised graph-based clustering method, which was visualized in the uniform manifold approximation and projection (UMAP) plot ( Figure 6D ).

A weighted kernel density estimation of the key DE-DRGs’ gene expression was performed using the Nebulosa (v1.10.0) package with default parameters. Figure 7A showed that MYL6 was expressed in most intrahepatic cells. In addition, DSTN was highly expressed mainly in cholangiocytes and ECs ( Figure 7B ). We generated a multi-dimensional scaling (MDS) plot using aggregated signal to explore similarities among samples. It is shown that the first dimension (MDS1) separates control and NAFLD samples, while the second dimension (MDS2) separates cell clusters ( Figure 7C ). To further clarify which cell clusters have altered the expression of key DE-DRGs in NAFLD patients compared to controls, we performed differential state analysis ( Figures 7D, E ; Supplementary Table S8 ). Interestingly, we found that MYL6 was significantly differentially expressed only in the MPs cluster (log2FC=0.597, p value=0.0124), and its expression trend was the same as in the bulk RNA samples (log2FC=1.247, p value<0.0001). In addition, DSTN was differentially expressed in the HSCs cluster (log2FC= -1.61, p value=0.00257) and MPs (log2FC= -0.543, p value=0.0298), but both showed the opposite trend of expression to that in bulk RNA samples (log2FC=0.596, p value<0.0001). Consequently, we focused on MYL6 for our subsequent investigation.

We performed re-decimation clustering of MPs into 5 distinct clusters ( Figure 8A ; Supplementary Figure S2 ). It emerged 2 clusters for monocytes (with a suffix of Mono) and 3 clusters for macrophages (with a suffix of Macro). The phylogenetic tree relating to the “average” cell revealed that the 2 Mono clusters were closely related as a branch node, as were the 3 Macro clusters ( Figure 8A right panel). MYL6 was highly expressed mainly in C01, C02 and C04 clusters ( Figure 8B ). Heatmap showed that C02-CD5L-Macro and C05-FCN3-Macro expressed high levels of Kupffer markers such as CD5L, MARCO, VCAM1, and CD163 ( Figure 8C ; Supplementary Table S9 ). In addition, by examining signature genes defined previously (28, 31, 32), we observed distinct functional status for each MPs subset, with the high Kupffer mod-score for C02-CD5L-Macro and C05-FCN3-Macro, high phagocytosis mod-score in the Macro clusters, similar proinflammatory mod-score in all clusters, highest M1 mod-score for C02-CD5L-Macro, and highest NASH-associated macrophages (NAMs) mod-score for C03-RNASE1-Macro ( Figure 8D ). Subsequently, differential state analysis in each cluster showed that MYL6 was significantly differentially expressed mainly in C01-S100A8-Mono (log2FC=0.524, p value=0.00489) and C02-CD5L-Macro (log2FC= 0.544, p value=0.02170) ( Figures 8E, F ; Supplementary Table S10 ). The results of GSEA enrichment analysis for C02-CD5L-Macro indicated that, in comparison to the control group, NAFLD patients exhibited up-regulation in the pro-inflammatory signaling pathway, the signaling pathways associated with regulated cell death (e.g., regulated necrosis and ferroptosis), and the signaling pathway related to sulfide metabolism (e.g., sulfide oxidation leading to sulfate) ( Figure 8G ; Supplementary Table S11 ). These results suggested that MYL6 might play an important role through M1 phenotype Kupffer cells in NAFLD patients.

We collected liver tissues from 16 patients and conducted HE staining as well as RNA-seq analysis ( Figure 9A ; Supplementary Table S12 ). The results indicated a substantial up-regulation of the majority of DRGs in the NAFLD group when compared to the control group ( Figure 9B ). In addition, the heatmap showed that MYL6 was significantly positively correlated with the expression of Kupffer cell signature genes, M1 phenotype macrophage signature genes, and oxidative stress-related genes ( Figure 9C ). Subsequently, we divided 10 NAFLD patients into high and low-expression groups based on the median MYL6 expression and performed GSEA between the two groups. The group with high expression of MYL6 showed an upregulation of inflammation-related pathways such as the inflame pathway, cytokines, and inflammatory response. They also showed an increase in fatty acid metabolism pathways and oxidative stress-related pathways ( Figure 9D ; Supplementary Table S13 ). It is noteworthy that there was a positive correlation between MYL6 expression and disease severity in NAFLD ( Figures 9E, F ). In addition, we downloaded data related to liver hepatocellular carcinoma (LIHC) from the TCGA database and categorized patients in whom the history of hepatocellular carcinoma risk factors was NAFLD only into the NAFLD-LIHC group. Our analysis showed that MYL6 expression was significantly higher in patients with NAFLD-LIHC than in controls ( Figure 9G ). Moreover, to validate the effectiveness of MYL6 in the diagnosis of NAFLD, ROC logistic regression analysis was performed based on our data. It was found that the AUC value of MYL6 was 0.967 ( Figure 9H ), indicating good sensitivity and specificity of MYL6 in the diagnosis of NAFLD. Further Spearman correlation analysis revealed that there was a significant positive correlation between intrahepatic MYL6 mRNA levels and NAS, as well as serum concentrations of TC, TG, and LDL-C. On the other hand, researchers discovered a negative connection between MYL6 expression and serum concentrations of HDL-C ( Figures 9I–M ).